/*------------------------------------------------------------------------------------------------*/

void dihedral_angle_generator_C(char seq_inc_anchors[], double dihedral_angles[], 
                              double first_O_dihedral, double last_psi, double last_O_dihedral)
{
    /*Choose a random set of dihedral angles from the dihedral data stored in dihedrals.h*/
    /*Resulting array of values contains the angles needed to grow the loop from the N-anchor*/
    
    int i;
    int j;
    int r;
    int no_residues = strlen(seq_inc_anchors);  /*= no. of loop residues + 2 anchor residues*/
    int no_atoms = (no_residues - 1) * 4;       /*Number of atoms to be added to the structure*/
    int array_index;
    double phi;
    double psi;
    double omega;
    
    for (i=0; i<no_residues; i++) {
        
        if      (seq_inc_anchors[no_residues-(1+i)] == 'A') { array_index =  0; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'R') { array_index =  1; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'N') { array_index =  2; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'D') { array_index =  3; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'C') { array_index =  4; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'E') { array_index =  5; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'Q') { array_index =  6; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'G') { array_index =  7; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'H') { array_index =  8; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'I') { array_index =  9; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'L') { array_index = 10; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'K') { array_index = 11; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'M') { array_index = 12; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'F') { array_index = 13; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'P') { array_index = 14; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'S') { array_index = 15; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'T') { array_index = 16; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'W') { array_index = 17; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'Y') { array_index = 18; }
        else if (seq_inc_anchors[no_residues-(1+i)] == 'V') { array_index = 19; }
        else {
            printf("Error: unknown residue type: %c\n", seq_inc_anchors[no_residues-(1+i)]);
            exit(0);
        }
            
        /*The dihedral angle required for the addition of the first N atom is calculated from the
          N-Ca-C=O dihedral angle, since the C, O and N need to be in the same plane*/
        if (i == no_residues - 1) {
            if (first_O_dihedral < 0) {
                dihedral_angles[i*4-1] = 180 + first_O_dihedral;
            }
            else if (first_O_dihedral > 0) {
                dihedral_angles[i*4-1] = first_O_dihedral - 180;
            }
            else {
                dihedral_angles[i*4-1] = 180;
            }
        }
        
        /*The other dihedral angles are taken from the distributions in the dihedrals.h file.*/
        else if (i > 0 && i < no_residues - 1) {
            r = rand() % dihedral_array_lengths[array_index] * 2;
            phi = *(ptr_array[array_index] + r);
            psi = *(ptr_array[array_index] + r + 1);
            omega = omega_angle_generator();
            
            dihedral_angles[i * 4 - 1] = psi + (rand() + 1.0) / (RAND_MAX + 1.0);
            dihedral_angles[i * 4] = (phi + (rand() + 1.0) / (RAND_MAX + 1.0));
            dihedral_angles[i * 4 + 1] = 180 + omega;
            dihedral_angles[i * 4 + 2] = omega;
        }
        
        /*The dihedral angles for addition of the C-anchor residues are now added (depend on the
          psi angle of the C-anchor residue, which can be calculated)*/
        else if (i == 0) {
            phi = find_dihedral_match(array_index, (int)roundf(last_psi), 1);
            omega = omega_angle_generator();
            dihedral_angles[0] = (phi + (rand() + 1.0) / (RAND_MAX + 1.0));
            dihedral_angles[1] = 180 + omega;
            dihedral_angles[2] = omega;
        }
    }
    
    /*Convert to radians*/
    for (i=0;i<no_atoms;i++){
        dihedral_angles[i] = dihedral_angles[i] * pi / 180;
    } 
}

/*------------------------------------------------------------------------------------------------*/

void bond_length_generator_C(char seq_inc_anchors[], double bond_lengths[], double N_anchor[][3])
{
    int i;
    int j;
    int no_residues = strlen(seq_inc_anchors);
    double bond[3];

    for (i=no_residues-2; i>0; i--) {
        if (seq_inc_anchors[i] == 'G') {
            if (seq_inc_anchors[i+1] == 'P') {
                bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.341, 0.016);
            }
            else {
                bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.329, 0.014);
            }
            bond_lengths[(no_residues-2-i)*4+1] = random_gaussian(1.231, 0.020);
            bond_lengths[(no_residues-2-i)*4+2] = random_gaussian(1.516, 0.018);
            bond_lengths[(no_residues-2-i)*4+3] = random_gaussian(1.451, 0.016);
        }
        else if (seq_inc_anchors[i] == 'P') {
            if (seq_inc_anchors[i+1] == 'P') {
                bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.341, 0.016);
            }
            else {
                bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.329, 0.014);
            }
            bond_lengths[(no_residues-2-i)*4+1] = random_gaussian(1.231, 0.020);
            bond_lengths[(no_residues-2-i)*4+2] = random_gaussian(1.525, 0.021);
            bond_lengths[(no_residues-2-i)*4+3] = random_gaussian(1.466, 0.015);
        }
        else {
            if (seq_inc_anchors[i+1] == 'P') {
                bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.341, 0.016);
            }
            else {
                bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.329, 0.014);
            }
            bond_lengths[(no_residues-2-i)*4+1] = random_gaussian(1.231, 0.020);
            bond_lengths[(no_residues-2-i)*4+2] = random_gaussian(1.525, 0.021);
            bond_lengths[(no_residues-2-i)*4+3] = random_gaussian(1.458, 0.019);
        }
    }

    i = 0;
    if (seq_inc_anchors[i+1] == 'P') {
        bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.341, 0.016);
    }
    else {
        bond_lengths[(no_residues-2-i)*4]   = random_gaussian(1.329, 0.014);
    }

    /*Calculate bond lengths from N anchor structure*/
    for (j=1; j<4; j++) {
        vector_subtract(N_anchor[4-j], N_anchor[3-j], bond);
        bond_lengths[(no_residues-2)*4+j] = magnitude(bond);
    }  
}

/*------------------------------------------------------------------------------------------------*/

void bond_angle_generator_C(char seq_inc_anchors[], double bond_angles[], double N_anchor[][3])
{
    /*Generates a set of random bond lengths from the distributions specified in the backbone_info
      header file*/
    
    double bond1[3],bond2[3];
    double cos_angle;
    int no_residues = strlen(seq_inc_anchors);
    int no_atoms = no_residues * 4;
    int i,j;
    
    for (i=no_residues-2; i>-1; i--) {
        if (seq_inc_anchors[i+1] == 'G') {
            bond_angles[(no_residues-2-i)*4] = random_gaussian(120.6, 1.7) * pi / 180;
            bond_angles[(no_residues-2-i)*4+1] = random_gaussian(123.0, 1.6) * pi / 180;
        }
        else if (seq_inc_anchors[i+1] == 'P') {
            bond_angles[(no_residues-2-i)*4] = random_gaussian(122.6, 5.0) * pi / 180;
            bond_angles[(no_residues-2-i)*4+1] = random_gaussian(122.0, 1.4) * pi / 180;
        }
        else {
            bond_angles[(no_residues-2-i)*4] = random_gaussian(121.7, 1.8) * pi / 180;
            bond_angles[(no_residues-2-i)*4+1] = random_gaussian(123.0, 1.6) * pi / 180;
        }

        if (seq_inc_anchors[i] == 'G') {
            if (seq_inc_anchors[i+1] == 'P') {
                bond_angles[(no_residues-2-i)*4+2] = random_gaussian(118.2, 2.1) * pi / 180;
            }
            else {
                bond_angles[(no_residues-2-i)*4+2] = random_gaussian(116.4, 2.1) * pi / 180;
            }
            bond_angles[(no_residues-2-i)*4+3] = random_gaussian(112.5, 2.9) * pi / 180;
        }
        else if (seq_inc_anchors[i] == 'P') {
            if (seq_inc_anchors[i+1] == 'P') {
                bond_angles[(no_residues-2-i)*4+2] = random_gaussian(116.9, 1.5) * pi / 180;
            }
            else {
                bond_angles[(no_residues-2-i)*4+2] = random_gaussian(116.2, 2.0) * pi / 180;
            }
            bond_angles[(no_residues-2-i)*4+3] = random_gaussian(111.8, 2.5) * pi / 180;
        }
        else {
            if (seq_inc_anchors[i+1] == 'P') {
                bond_angles[(no_residues-2-i)*4+2] = random_gaussian(116.9, 1.5) * pi / 180;
            }
            else {
                bond_angles[(no_residues-2-i)*4+2] = random_gaussian(116.2, 2.0) * pi / 180;
            }
            bond_angles[(no_residues-2-i)*4+3] = random_gaussian(111.2, 2.8) * pi / 180;
        }
    }

    /*Calculate N-Ca-C bond angle from N-anchor structure*/
    vector_subtract(N_anchor[1], N_anchor[0], bond1);
    vector_subtract(N_anchor[1], N_anchor[2], bond2);

    cos_angle = dot_product(bond1, bond2) / (magnitude(bond1) * magnitude(bond2));

    bond_angles[(no_residues-2)*4+3] = acos(cos_angle);
    
}    

/*------------------------------------------------------------------------------------------------*/

void build_loop_C(char seq_inc_anchors[], double N_anchor[][3], double C_anchor[][3], 
                double **loop_model, double first_O_dihedral, double last_psi, 
                double last_O_dihedral, double FREADbond_lengths[], double FREADbond_angles[], 
                double FREADdihedral_angles[], char bitstring[])
{
    /*Get a random set of dihedral angles, randomise bond lengths and angles, and use them to build
      a loop. Atoms added in the following order: C, O, Ca, N*/
      
    int no_residues = strlen(seq_inc_anchors)-1; /*No. of residues to be added, includes N anchor*/
    int no_atoms = no_residues * 4;
    int A_index,B_index,C_index;
    double dihedral_angles[no_atoms];
    double bond_lengths[no_atoms];
    double bond_angles[no_atoms];
    double length,angle,dihedral;
    double A[3], B[3], C[3], D[3], D2[3];
    double AB[3], BC[3], BC_unit[3], n[3], n_unit[3];
    double BC_magnitude, n_magnitude;
    double M[3][3];
    int i,j;
    int res;
    
    bond_length_generator_C(seq_inc_anchors, bond_lengths, N_anchor);
    
    bond_angle_generator_C(seq_inc_anchors, bond_angles, N_anchor);
    
    dihedral_angle_generator_C(seq_inc_anchors, dihedral_angles, first_O_dihedral, last_psi,
                             last_O_dihedral);
                                                          
    for (i=0; i<no_atoms; i++) {
        if (FREADbond_lengths[i] > 0.1) {
            bond_lengths[i] = FREADbond_lengths[i];
            dihedral_angles[i] = FREADdihedral_angles[i];
        }
        
        if (FREADbond_angles[i] > 0.1) {
            bond_angles[i] = FREADbond_angles[i];
        }
    }
                                 
    /*Fill in the first residue of the unclosed loop (the C anchor residue)*/
    for (i=0; i<3; i++) { 
        loop_model[0][i] = C_anchor[2][i];   /*C*/
    }
    
    for (i=0; i<3; i++) { 
        loop_model[1][i] = C_anchor[3][i];   /*O*/
    }
    
    for (i=0; i<3; i++) { 
        loop_model[2][i] = C_anchor[1][i];   /*Ca*/
    }
    
    for (i=0; i<3; i++) { 
        loop_model[3][i] = C_anchor[0][i];   /*N*/
    }

    /*Loop Building (see Parsons et al., 2005 DOI 10.1002/jcc.20237)*/
    
    for (i=0;i<no_atoms;i++) {

        length = bond_lengths[i];
        angle = pi - bond_angles[i];
        dihedral = dihedral_angles[i];
        
        /*Get required indices (of loop_model array)*/
        if (i % 4 == 0) {
            /*C*/
            A_index = i;
            B_index = i + 2;
            C_index = i + 3;
        }
        else if ((i - 1) % 4 == 0) {
            /*O*/
            A_index = i + 1;
            B_index = i + 2;
            C_index = i + 3;
        }
        else if ((i - 2) % 4 == 0) {
            /*Ca*/
            A_index = i;
            B_index = i + 1;
            C_index = i + 2;
        }
        else if ((i + 1) % 4 == 0) {
            /*N*/
            A_index = i;
            B_index = i + 1;
            C_index = i + 3;
        }
        
        /*Extract required atom coordinates from loop_model array*/
        for (j=0;j<3;j++) {
            A[j] = loop_model[A_index][j];
            B[j] = loop_model[B_index][j];
            C[j] = loop_model[C_index][j];
        }
               
        D2[0] = length * cos(angle);
        D2[1] = length * cos(dihedral) * sin(angle);
        D2[2] = length * sin(dihedral) * sin(angle);

        /*Calculate AB vector*/
        vector_subtract(A, B, AB);
  
        /*Calculate BC unit vector*/
        vector_subtract(B, C, BC);
     
        BC_magnitude = magnitude(BC);
    
        for (j=0; j<3; j++) {
            BC_unit[j] = BC[j] / BC_magnitude;
        }
     
        /*Calculate normal to the ABC plane*/
        cross_product(AB, BC_unit, n);

        n_magnitude = magnitude(n);
       
        for (j=0; j<3; j++) {
            n_unit[j] = n[j] / n_magnitude;
        }
       
        /*Insert values into M matrix*/
        for (j=0; j<3; j++) {
            M[j][0] = BC_unit[j];
            M[j][2] = n_unit[j];
        }
        
        M[0][1] = n_unit[1]*BC_unit[2] - n_unit[2]*BC_unit[1];
        M[1][1] = n_unit[2]*BC_unit[0] - n_unit[0]*BC_unit[2];
        M[2][1] = n_unit[0]*BC_unit[1] - n_unit[1]*BC_unit[0];
      
        /*Calculate atom position*/
        loop_model[i+4][0] = (M[0][0]*D2[0]) + (M[0][1]*D2[1]) + (M[0][2]*D2[2]) + C[0];
        loop_model[i+4][1] = (M[1][0]*D2[0]) + (M[1][1]*D2[1]) + (M[1][2]*D2[2]) + C[1];
        loop_model[i+4][2] = (M[2][0]*D2[0]) + (M[2][1]*D2[1]) + (M[2][2]*D2[2]) + C[2];
        
    }
}

/*------------------------------------------------------------------------------------------------*/

void CCD_C(double loop_model[][3], int i, double F1[3], double F2[3], double F3[3], int A_indices[],
         int loop_length)
{
    int j;
    double A[3], B[3];
    double M01[3], M02[3], M03[3];
    double O1[3], O2[3], O3[3];
    double f1[3], f2[3], f3[3];
    double r1[3], r2[3], r3[3];
    double r1_mag, r2_mag, r3_mag;
    double r1_unit[3], r2_unit[3], r3_unit[3];
    double theta[3];
    double theta_mag;
    double theta_unit[3];
    double s1_unit[3], s2_unit[3], s3_unit[3];
    double b, c, cos_angle, sin_angle, angle;
    double angle_limit = 3 * pi / 180;
    double R[3][3];
    double new_x, new_y, new_z;
    double d, e, f, u, v, w, x, y, z; 
    
    /*Find dihedral angle needed to minimise C anchor fixed/moving RMSD*/
    for (j=0; j<3; j++) {
        A[j] = loop_model[A_indices[i]][j];
        B[j] = loop_model[A_indices[i]+1][j];
        M01[j] = loop_model[(loop_length + 1) * 4 + 1][j];   /*C of moving N anchor*/
        M02[j] = loop_model[(loop_length + 1) * 4 + 2][j];   /*Ca of moving N anchor*/
        M03[j] = loop_model[(loop_length + 1) * 4 + 3][j];   /*N of moving N anchor*/
    }

    /*Calculate O position vectors*/
    find_foot_of_perpendicular(A, B, M01, O1);
    find_foot_of_perpendicular(A, B, M02, O2);
    find_foot_of_perpendicular(A, B, M03, O3);
    
    /*Calculate 'f' vectors*/       
    vector_subtract(O1, F1, f1);
    vector_subtract(O2, F2, f2);
    vector_subtract(O3, F3, f3);
           
    /*Calculate 'r' unit vectors*/
    vector_subtract(O1, M01, r1);
    vector_subtract(O2, M02, r2);
    vector_subtract(O3, M03, r3);
    
    r1_mag = magnitude(r1);
    r2_mag = magnitude(r2);
    r3_mag = magnitude(r3);
    
    for (j=0; j<3; j++) {
        r1_unit[j] = r1[j] / r1_mag;
        r2_unit[j] = r2[j] / r2_mag;
        r3_unit[j] = r3[j] / r3_mag;
    }
    
    /*Calculate unit vector along the rotation axis*/
    vector_subtract(O1, O2, theta);
    
    theta_mag = magnitude(theta);
    
    for (j=0; j<3; j++) {
        theta_unit[j] = theta[j] / theta_mag;
    }
    
    /*Calculate 's' vectors*/
    cross_product(r1_unit, theta_unit, s1_unit);
    cross_product(r2_unit, theta_unit, s2_unit);
    cross_product(r3_unit, theta_unit, s3_unit);
    
    /*Calculate cos(angle) and sin(angle) - see CCD paper*/
    b = 2 * r1_mag * dot_product(f1, r1_unit) + 2 * r2_mag * dot_product(f2, r2_unit) +
        2 * r3_mag * dot_product(f3, r3_unit);
        
    c = 2 * r1_mag * dot_product(f1, s1_unit) + 2 * r2_mag * dot_product(f2, s2_unit) +
        2 * r3_mag * dot_product(f3, s3_unit);
        
    cos_angle = b / sqrt(pow(b,2) + pow(c,2));
    sin_angle = c / sqrt(pow(b,2) + pow(c,2));
    
    /*Calculate angle to rotate by, limited to 3 degrees*/
    angle = atan2(sin_angle, cos_angle);
    
    /*if (angle < 0 && angle < -angle_limit) {
        angle = -angle_limit;
    }
    else if (angle > 0 && angle > angle_limit) {
        angle = angle_limit;
    }*/
            
    /*Rotate all atoms after the AB bond by this angle about the theta axis*/
    /*http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/, section 6.2*/
    d = O1[0];
    e = O1[1];
    f = O1[2];
    u = theta_unit[0];
    v = theta_unit[1];
    w = theta_unit[2];
    
    for (j=A_indices[i]+2; j<(loop_length + 2) * 4; j++) {
        x = loop_model[j][0];
        y = loop_model[j][1];
        z = loop_model[j][2];
        
        loop_model[j][0] = (d*(v*v + w*w) - u*(e*v + f*w - u*x - v*y - w*z)) * (1-cos(-angle)) 
                              + x*cos(-angle) + (-f*v + e*w - w*y + v*z) * sin(-angle);
                               
        loop_model[j][1] = (e*(u*u + w*w) - v*(d*u + f*w - u*x - v*y - w*z)) * (1-cos(-angle)) 
                              + y*cos(-angle) + (f*u - d*w + w*x - u*z) * sin(-angle);
                              
        loop_model[j][2] = (f*(u*u + v*v) - w*(d*u + e*v - u*x - v*y - w*z)) * (1-cos(-angle)) 
                              + z*cos(-angle) + (-e*u + d*v - v*x + u*y) * sin(-angle);
    }
}

/*------------------------------------------------------------------------------------------------*/

int close_loop_C_first(double **loop_model, double fixed_N_anchor[][3], int loop_length,
                 char seq_inc_anchors[], char bitstring[])
{
    /*CCD: Iterate through dihedrals, minimising the difference between the fixed and moving C 
      anchor residues (see "Cyclic coordinate descent: A robotics algorithm for protein loop
      closure", A. A. Canutescu and R. L. Dunbrack Jr., Protein Science, 12:963-972, 2003*/

    int A_indices[loop_length*2];
    int i, j, k, m;
    int iteration_count;
    int index;
    int P_old_index, P_new_index;
    double fixed_N_anchor_reordered[4][3];
    double moving_N_anchor[4][3];
    double N_anchor_RMSD;
    double RMSD_cutoff = 1.0;
    double F1[3], F2[3], F3[3];
    double initial_N_RMSD;
    double proposed_model[(loop_length+2)*4][3];
    double old_phi, old_psi;
    double new_phi, new_psi;
    double P_old, P_new;
    double P;
    double r;
    double probabilities[loop_length];
    double sum_probs;
    
        
    /*The central bond of the dihedral angle (i.e. the vector around which we rotate) is between 
      atoms A and B - get indices of the A atoms (this effectively creates a list containing which
      dihedral angles should be changed)*/
    
    for (i=0; i<loop_length * 2; i++) {
        if (i % 2 == 0) {
            A_indices[i] = (i / 2 * 4) + 5;
        }
        else {
            A_indices[i] = A_indices[i-1] + 1;
        }
    }
    
    /*Find initial RMSD between fixed and moving N anchor atoms*/
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) {
            moving_N_anchor[i][j] = loop_model[(loop_length + 1)*4 + i][j];
            F1[j] = fixed_N_anchor[2][j];    /*C of fixed N anchor*/
            F2[j] = fixed_N_anchor[1][j];    /*Ca of fixed N anchor*/
            F3[j] = fixed_N_anchor[0][j];    /*N of fixed N anchor*/
        }
    } 
    
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) { 
            fixed_N_anchor_reordered[i][j] = fixed_N_anchor[3-i][j];
        }
    }
    
    N_anchor_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);
    initial_N_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);
    
    i = 0;
    iteration_count = 0;
    
    for (j=0;j<(loop_length+2)*4;j++) {
        proposed_model[j][0] = loop_model[j][0];
        proposed_model[j][1] = loop_model[j][1];
        proposed_model[j][2] = loop_model[j][2];
    }
 
    /*CCD (number of iterations limited)*/
    while (N_anchor_RMSD > RMSD_cutoff) {  
    
        //Allow residues next to MECHANO-modelled residues to change--------------------------------------------------------
        
        //If this is the first residue:
        /*if ((loop_length-1) - (A_indices[i]-5)/4 < 0) {
            if (bitstring[(loop_length-1) - (A_indices[i]-5)/4 + 1] == '1' && bitstring[(loop_length-1) - (A_indices[i]-5)/4] == '1') {
                
                if (i == loop_length * 2 - 2) {
                    i = 0;
                    iteration_count += 1;
                }
                else {
                    i += 2;
                }
                
                continue;
            }
        }
        
        //If this is the last residue:
        if ((loop_length-1) - (A_indices[i]-5)/4 > loop_length-1) {
            if (bitstring[(loop_length-1) - (A_indices[i]-5)/4 -1] == '1' && bitstring[(loop_length-1) - (A_indices[i]-5)/4] == '1') {
                
                if (i == loop_length * 2 - 2) {
                    i = 0;
                    iteration_count += 1;
                }
                else {
                    i += 2;
                }
                
                continue;
            }
        }
        
        //If this is other residues:
        if (bitstring[(loop_length-1) - (A_indices[i]-5)/4 + 1] == '1' && bitstring[(loop_length-1) - (A_indices[i]-5)/4] == '1' && bitstring[(loop_length-1) - (A_indices[i]-5)/4 - 1] == '1') {
            
            if (i == loop_length * 2 - 2) {
                i = 0;
                iteration_count += 1;
            }
            else {
                i += 2;
            }
            
            continue;
        }*/
        //------------------------------------------------------------------------------------------------------------------
        
        //Ignore residues modelled using FREAD------------------------------------------------------------------------------
        
        if (bitstring[(loop_length-1) - (A_indices[i]-5)/4] == '1') {
            /*Set up next iteration*/
            if (i == loop_length * 2 - 2) {
                i = 0;
                iteration_count += 1;
            }
            else {
                i += 2;
            }
            continue;
        }
        
        //------------------------------------------------------------------------------------------------------------------
                                    
        CCD_C(proposed_model, i, F1, F2, F3, A_indices, loop_length);
        
        CCD_C(proposed_model, i+1, F1, F2, F3, A_indices, loop_length);
                    
        /*Calculate previous dihedrals*/
        index = A_indices[i];
        old_psi = calculate_dihedral(loop_model[index-2], loop_model[index],
                                 loop_model[index+1], loop_model[index+2]);
                                 
        old_phi = calculate_dihedral(loop_model[index], loop_model[index+1],
                                 loop_model[index+2], loop_model[index+4]);
                                 
        /*Calculate new dihedrals*/
        new_psi = calculate_dihedral(proposed_model[index-2], proposed_model[index],
                                     proposed_model[index+1], proposed_model[index+2]);
                                 
        new_phi = calculate_dihedral(proposed_model[index], proposed_model[index+1],
                                     proposed_model[index+2], proposed_model[index+4]);
                                 
        /*Decide whether to accept/reject*/
        P_old_index = (int)((180+old_phi)/5)*72 + (int)((180+old_psi)/5);
        P_new_index = (int)((180+new_phi)/5)*72 + (int)((180+new_psi)/5);
        
        if      (seq_inc_anchors[loop_length+1-(index/4)] == 'A') { P_old = ALA_CCD_dist[P_old_index];
                                                    P_new = ALA_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'R') { P_old = ARG_CCD_dist[P_old_index];
                                                    P_new = ARG_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'N') { P_old = ASN_CCD_dist[P_old_index];
                                                    P_new = ASN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'D') { P_old = ASP_CCD_dist[P_old_index];
                                                    P_new = ASP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'C') { P_old = CYS_CCD_dist[P_old_index];
                                                    P_new = CYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'E') { P_old = GLU_CCD_dist[P_old_index];
                                                    P_new = GLU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'Q') { P_old = GLN_CCD_dist[P_old_index];
                                                    P_new = GLN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'G') { P_old = GLY_CCD_dist[P_old_index];
                                                    P_new = GLY_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'H') { P_old = HIS_CCD_dist[P_old_index];
                                                    P_new = HIS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'I') { P_old = ILE_CCD_dist[P_old_index];
                                                    P_new = ILE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'L') { P_old = LEU_CCD_dist[P_old_index];
                                                    P_new = LEU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'K') { P_old = LYS_CCD_dist[P_old_index];
                                                    P_new = LYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'M') { P_old = MET_CCD_dist[P_old_index];
                                                    P_new = MET_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'F') { P_old = PHE_CCD_dist[P_old_index];
                                                    P_new = PHE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'P') { P_old = PRO_CCD_dist[P_old_index];
                                                    P_new = PRO_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'S') { P_old = SER_CCD_dist[P_old_index];
                                                    P_new = SER_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'T') { P_old = THR_CCD_dist[P_old_index];
                                                    P_new = THR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'W') { P_old = TRP_CCD_dist[P_old_index];
                                                    P_new = TRP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'Y') { P_old = TYR_CCD_dist[P_old_index];
                                                    P_new = TYR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'V') { P_old = VAL_CCD_dist[P_old_index];
                                                    P_new = VAL_CCD_dist[P_new_index]; }
        else {
            printf("Error: unknown residue type: %c\n", seq_inc_anchors[loop_length+1-(index/4)]);
            exit(0);
        }
        
        probabilities[i/2] = P_new;
        
        /*Decide whether to accept the new angles or not*/
        if (P_new >= P_old) {  
            for (j=0;j<(loop_length+2)*4;j++) {
                loop_model[j][0] = proposed_model[j][0];
                loop_model[j][1] = proposed_model[j][1];
                loop_model[j][2] = proposed_model[j][2];
            }   
        }             
        else {
            P = P_new/P_old;
            r = (double)rand()/(double)RAND_MAX;
            if (r <= P) {
                for (j=0;j<(loop_length+2)*4;j++) {
                    loop_model[j][0] = proposed_model[j][0];
                    loop_model[j][1] = proposed_model[j][1];
                    loop_model[j][2] = proposed_model[j][2];
                }
            }
            else {
                for (j=0;j<(loop_length+2)*4;j++) {
                    proposed_model[j][0] = loop_model[j][0];
                    proposed_model[j][1] = loop_model[j][1];
                    proposed_model[j][2] = loop_model[j][2];
                }
            }                
        }       
                    
        for (j=0; j<4; j++) {
            for (k=0; k<3; k++) {
                moving_N_anchor[j][k] = loop_model[(loop_length + 1)*4 + j][k];
            }
        } 
    
        N_anchor_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);

        if (iteration_count >= 100) {
            return 0;
        } 
        
        sum_probs = 0;
        if (iteration_count > 0) {
            for (j=0; j<loop_length;j++) {
                sum_probs = sum_probs + probabilities[j];
            }
            if (sum_probs == 0) {
                return 0; 
            }
        }
           
        /*Set up next iteration*/
        if (i == loop_length * 2 - 2) {
            i = 0;
            iteration_count += 1;
        }
        else {
            i += 2;
        }
    }
        
    /*Insert the true N-anchor coordinates*/
    /*for (j=0; j<4; j++) {
        loop_model[(loop_length + 1) * 4 + j][0] = fixed_N_anchor_reordered[j][0];
        loop_model[(loop_length + 1) * 4 + j][1] = fixed_N_anchor_reordered[j][1];
        loop_model[(loop_length + 1) * 4 + j][2] = fixed_N_anchor_reordered[j][2];
    }*/

    return 1;
}

/*------------------------------------------------------------------------------------------------*/

int close_loop_C(double **loop_model, double fixed_N_anchor[][3], int loop_length,
                 char seq_inc_anchors[])
{
    /*CCD: Iterate through dihedrals, minimising the difference between the fixed and moving C 
      anchor residues (see "Cyclic coordinate descent: A robotics algorithm for protein loop
      closure", A. A. Canutescu and R. L. Dunbrack Jr., Protein Science, 12:963-972, 2003*/

    int A_indices[loop_length*2];
    int i, j, k, m;
    int iteration_count;
    int index;
    int P_old_index, P_new_index;
    double fixed_N_anchor_reordered[4][3];
    double moving_N_anchor[4][3];
    double N_anchor_RMSD;
    double RMSD_cutoff = 0.2;
    double F1[3], F2[3], F3[3];
    double initial_N_RMSD;
    double proposed_model[(loop_length+2)*4][3];
    double old_phi, old_psi;
    double new_phi, new_psi;
    double P_old, P_new;
    double P;
    double r;
    double probabilities[loop_length];
    double sum_probs;
    
        
    /*The central bond of the dihedral angle (i.e. the vector around which we rotate) is between 
      atoms A and B - get indices of the A atoms (this effectively creates a list containing which
      dihedral angles should be changed)*/
    
    for (i=0; i<loop_length * 2; i++) {
        if (i % 2 == 0) {
            A_indices[i] = (i / 2 * 4) + 5;
        }
        else {
            A_indices[i] = A_indices[i-1] + 1;
        }
    }
    
    /*Find initial RMSD between fixed and moving N anchor atoms*/
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) {
            moving_N_anchor[i][j] = loop_model[(loop_length + 1)*4 + i][j];
            F1[j] = fixed_N_anchor[2][j];    /*C of fixed N anchor*/
            F2[j] = fixed_N_anchor[1][j];    /*Ca of fixed N anchor*/
            F3[j] = fixed_N_anchor[0][j];    /*N of fixed N anchor*/
        }
    } 
    
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) { 
            fixed_N_anchor_reordered[i][j] = fixed_N_anchor[3-i][j];
        }
    }
    
    N_anchor_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);
    initial_N_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);
    
    i = 0;
    iteration_count = 0;
    
    for (j=0;j<(loop_length+2)*4;j++) {
        proposed_model[j][0] = loop_model[j][0];
        proposed_model[j][1] = loop_model[j][1];
        proposed_model[j][2] = loop_model[j][2];
    }
 
    /*CCD (number of iterations limited)*/
    while (N_anchor_RMSD > RMSD_cutoff) {                                
                                
        CCD_C(proposed_model, i, F1, F2, F3, A_indices, loop_length);
        
        CCD_C(proposed_model, i+1, F1, F2, F3, A_indices, loop_length);
                    
        /*Calculate previous dihedrals*/
        index = A_indices[i];
        old_psi = calculate_dihedral(loop_model[index-2], loop_model[index],
                                 loop_model[index+1], loop_model[index+2]);
                                 
        old_phi = calculate_dihedral(loop_model[index], loop_model[index+1],
                                 loop_model[index+2], loop_model[index+4]);
                                 
        /*Calculate new dihedrals*/
        new_psi = calculate_dihedral(proposed_model[index-2], proposed_model[index],
                                     proposed_model[index+1], proposed_model[index+2]);
                                 
        new_phi = calculate_dihedral(proposed_model[index], proposed_model[index+1],
                                     proposed_model[index+2], proposed_model[index+4]);
                                 
        /*Decide whether to accept/reject*/
        P_old_index = (int)((180+old_phi)/5)*72 + (int)((180+old_psi)/5);
        P_new_index = (int)((180+new_phi)/5)*72 + (int)((180+new_psi)/5);
        
        if      (seq_inc_anchors[loop_length+1-(index/4)] == 'A') { P_old = ALA_CCD_dist[P_old_index];
                                                    P_new = ALA_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'R') { P_old = ARG_CCD_dist[P_old_index];
                                                    P_new = ARG_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'N') { P_old = ASN_CCD_dist[P_old_index];
                                                    P_new = ASN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'D') { P_old = ASP_CCD_dist[P_old_index];
                                                    P_new = ASP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'C') { P_old = CYS_CCD_dist[P_old_index];
                                                    P_new = CYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'E') { P_old = GLU_CCD_dist[P_old_index];
                                                    P_new = GLU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'Q') { P_old = GLN_CCD_dist[P_old_index];
                                                    P_new = GLN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'G') { P_old = GLY_CCD_dist[P_old_index];
                                                    P_new = GLY_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'H') { P_old = HIS_CCD_dist[P_old_index];
                                                    P_new = HIS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'I') { P_old = ILE_CCD_dist[P_old_index];
                                                    P_new = ILE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'L') { P_old = LEU_CCD_dist[P_old_index];
                                                    P_new = LEU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'K') { P_old = LYS_CCD_dist[P_old_index];
                                                    P_new = LYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'M') { P_old = MET_CCD_dist[P_old_index];
                                                    P_new = MET_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'F') { P_old = PHE_CCD_dist[P_old_index];
                                                    P_new = PHE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'P') { P_old = PRO_CCD_dist[P_old_index];
                                                    P_new = PRO_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'S') { P_old = SER_CCD_dist[P_old_index];
                                                    P_new = SER_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'T') { P_old = THR_CCD_dist[P_old_index];
                                                    P_new = THR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'W') { P_old = TRP_CCD_dist[P_old_index];
                                                    P_new = TRP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'Y') { P_old = TYR_CCD_dist[P_old_index];
                                                    P_new = TYR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'V') { P_old = VAL_CCD_dist[P_old_index];
                                                    P_new = VAL_CCD_dist[P_new_index]; }
        else {
            printf("Error: unknown residue type: %c\n", seq_inc_anchors[loop_length+1-(index/4)]);
            exit(0);
        }
        
        probabilities[i/2] = P_new;
        
        /*Decide whether to accept the new angles or not*/
        if (P_new >= P_old) {  
            for (j=0;j<(loop_length+2)*4;j++) {
                loop_model[j][0] = proposed_model[j][0];
                loop_model[j][1] = proposed_model[j][1];
                loop_model[j][2] = proposed_model[j][2];
            }   
        }             
        else {
            P = P_new/P_old;
            r = (double)rand()/(double)RAND_MAX;
            if (r <= P) {
                for (j=0;j<(loop_length+2)*4;j++) {
                    loop_model[j][0] = proposed_model[j][0];
                    loop_model[j][1] = proposed_model[j][1];
                    loop_model[j][2] = proposed_model[j][2];
                }
            }
            else {
                for (j=0;j<(loop_length+2)*4;j++) {
                    proposed_model[j][0] = loop_model[j][0];
                    proposed_model[j][1] = loop_model[j][1];
                    proposed_model[j][2] = loop_model[j][2];
                }
            }                
        }       
                    
        for (j=0; j<4; j++) {
            for (k=0; k<3; k++) {
                moving_N_anchor[j][k] = loop_model[(loop_length + 1)*4 + j][k];
            }
        } 
    
        N_anchor_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);

        if (iteration_count >= 100) {
            return 0;
        } 
        
        sum_probs = 0;
        if (iteration_count > 0) {
            for (j=0; j<loop_length;j++) {
                sum_probs = sum_probs + probabilities[j];
            }
            if (sum_probs == 0) {
                return 0; 
            }
        }
           
        /*Set up next iteration*/
        if (i == loop_length * 2 - 2) {
            i = 0;
            iteration_count += 1;
        }
        else {
            i += 2;
        }
    }
        
    /*Insert the true N-anchor coordinates*/
    for (j=0; j<4; j++) {
        loop_model[(loop_length + 1) * 4 + j][0] = fixed_N_anchor_reordered[j][0];
        loop_model[(loop_length + 1) * 4 + j][1] = fixed_N_anchor_reordered[j][1];
        loop_model[(loop_length + 1) * 4 + j][2] = fixed_N_anchor_reordered[j][2];
    }

    return 1;
}

/*------------------------------------------------------------------------------------------------*/

int close_loop_C_limited(double **loop_model, double fixed_N_anchor[][3], int loop_length,
                 char seq_inc_anchors[], char bitstring[], double FREADdihedralsN[])
{
    /*CCD: Iterate through dihedrals, minimising the difference between the fixed and moving C 
      anchor residues (see "Cyclic coordinate descent: A robotics algorithm for protein loop
      closure", A. A. Canutescu and R. L. Dunbrack Jr., Protein Science, 12:963-972, 2003*/

    int A_indices[loop_length*2];
    int i, j, k, m;
    int iteration_count;
    int index;
    int P_old_index, P_new_index;
    int residue_no;
    int FREAD_modelled;
    double fixed_N_anchor_reordered[4][3];
    double moving_N_anchor[4][3];
    double N_anchor_RMSD;
    double RMSD_cutoff = 1.0;
    double F1[3], F2[3], F3[3];
    double initial_N_RMSD;
    double proposed_model[(loop_length+2)*4][3];
    double old_phi, old_psi;
    double new_phi, new_psi;
    double initial_phi, initial_psi;
    double P_old, P_new;
    double P;
    double r;
    double probabilities[loop_length];
    double sum_probs;
    
        
    /*The central bond of the dihedral angle (i.e. the vector around which we rotate) is between 
      atoms A and B - get indices of the A atoms (this effectively creates a list containing which
      dihedral angles should be changed)*/
    
    for (i=0; i<loop_length * 2; i++) {
        if (i % 2 == 0) {
            A_indices[i] = (i / 2 * 4) + 5;
        }
        else {
            A_indices[i] = A_indices[i-1] + 1;
        }
    }
    
    /*Find initial RMSD between fixed and moving N anchor atoms*/
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) {
            moving_N_anchor[i][j] = loop_model[(loop_length + 1)*4 + i][j];
            F1[j] = fixed_N_anchor[2][j];    /*C of fixed N anchor*/
            F2[j] = fixed_N_anchor[1][j];    /*Ca of fixed N anchor*/
            F3[j] = fixed_N_anchor[0][j];    /*N of fixed N anchor*/
        }
    } 
    
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) { 
            fixed_N_anchor_reordered[i][j] = fixed_N_anchor[3-i][j];
        }
    }
    
    N_anchor_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);
    initial_N_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);
    
    i = 0;
    iteration_count = 0;
    
    for (j=0;j<(loop_length+2)*4;j++) {
        proposed_model[j][0] = loop_model[j][0];
        proposed_model[j][1] = loop_model[j][1];
        proposed_model[j][2] = loop_model[j][2];
    }
 
    /*CCD (number of iterations limited)*/
    while (N_anchor_RMSD > RMSD_cutoff) {  
            
        //Was this residue modelled using FREAD?------------------------------------------------------------------------------
        
        residue_no = (loop_length-1) - (A_indices[i]-5)/4;
        if (bitstring[residue_no] == '1') {
                FREAD_modelled = 1;
        }
        else {
            FREAD_modelled = 0;
        }
                                                 
        CCD_C(proposed_model, i, F1, F2, F3, A_indices, loop_length);
        
        CCD_C(proposed_model, i+1, F1, F2, F3, A_indices, loop_length);
                    
        /*Calculate previous dihedrals*/
        index = A_indices[i];
        old_psi = calculate_dihedral(loop_model[index-2], loop_model[index],
                                 loop_model[index+1], loop_model[index+2]);
                                 
        old_phi = calculate_dihedral(loop_model[index], loop_model[index+1],
                                 loop_model[index+2], loop_model[index+4]);
                                 
        /*Calculate new dihedrals*/
        new_psi = calculate_dihedral(proposed_model[index-2], proposed_model[index],
                                     proposed_model[index+1], proposed_model[index+2]);
                                 
        new_phi = calculate_dihedral(proposed_model[index], proposed_model[index+1],
                                     proposed_model[index+2], proposed_model[index+4]);
                                     
        /*If this residue was modelled using FREAD, how far has it moved from where it started?*/
        if (FREAD_modelled == 1) {
            initial_phi = FREADdihedralsN[(residue_no*4)+2];
            initial_psi = FREADdihedralsN[(residue_no+1)*4];
            
            if (fabs(initial_phi-new_phi) > 0.5236 || fabs(initial_psi-new_psi) > 0.5236) {
                for (j=0;j<(loop_length+2)*4;j++) {
                    proposed_model[j][0] = loop_model[j][0];
                    proposed_model[j][1] = loop_model[j][1];
                    proposed_model[j][2] = loop_model[j][2];
                }
                
                if (i == loop_length * 2 - 2) {
                    i = 0;
                    iteration_count += 1;
                }
                else {
                    i += 2;
                }
                
                continue;
            }               
        }
                                 
        /*Decide whether to accept/reject*/
        P_old_index = (int)((180+old_phi)/5)*72 + (int)((180+old_psi)/5);
        P_new_index = (int)((180+new_phi)/5)*72 + (int)((180+new_psi)/5);
        
        if      (seq_inc_anchors[loop_length+1-(index/4)] == 'A') { P_old = ALA_CCD_dist[P_old_index];
                                                    P_new = ALA_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'R') { P_old = ARG_CCD_dist[P_old_index];
                                                    P_new = ARG_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'N') { P_old = ASN_CCD_dist[P_old_index];
                                                    P_new = ASN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'D') { P_old = ASP_CCD_dist[P_old_index];
                                                    P_new = ASP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'C') { P_old = CYS_CCD_dist[P_old_index];
                                                    P_new = CYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'E') { P_old = GLU_CCD_dist[P_old_index];
                                                    P_new = GLU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'Q') { P_old = GLN_CCD_dist[P_old_index];
                                                    P_new = GLN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'G') { P_old = GLY_CCD_dist[P_old_index];
                                                    P_new = GLY_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'H') { P_old = HIS_CCD_dist[P_old_index];
                                                    P_new = HIS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'I') { P_old = ILE_CCD_dist[P_old_index];
                                                    P_new = ILE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'L') { P_old = LEU_CCD_dist[P_old_index];
                                                    P_new = LEU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'K') { P_old = LYS_CCD_dist[P_old_index];
                                                    P_new = LYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'M') { P_old = MET_CCD_dist[P_old_index];
                                                    P_new = MET_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'F') { P_old = PHE_CCD_dist[P_old_index];
                                                    P_new = PHE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'P') { P_old = PRO_CCD_dist[P_old_index];
                                                    P_new = PRO_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'S') { P_old = SER_CCD_dist[P_old_index];
                                                    P_new = SER_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'T') { P_old = THR_CCD_dist[P_old_index];
                                                    P_new = THR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'W') { P_old = TRP_CCD_dist[P_old_index];
                                                    P_new = TRP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'Y') { P_old = TYR_CCD_dist[P_old_index];
                                                    P_new = TYR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[loop_length+1-(index/4)] == 'V') { P_old = VAL_CCD_dist[P_old_index];
                                                    P_new = VAL_CCD_dist[P_new_index]; }
        else {
            printf("Error: unknown residue type: %c\n", seq_inc_anchors[loop_length+1-(index/4)]);
            exit(0);
        }
        
        probabilities[i/2] = P_new;
        
        /*Decide whether to accept the new angles or not*/
        if (P_new >= P_old) {  
            for (j=0;j<(loop_length+2)*4;j++) {
                loop_model[j][0] = proposed_model[j][0];
                loop_model[j][1] = proposed_model[j][1];
                loop_model[j][2] = proposed_model[j][2];
            }   
        }             
        else {
            P = P_new/P_old;
            r = (double)rand()/(double)RAND_MAX;
            if (r <= P) {
                for (j=0;j<(loop_length+2)*4;j++) {
                    loop_model[j][0] = proposed_model[j][0];
                    loop_model[j][1] = proposed_model[j][1];
                    loop_model[j][2] = proposed_model[j][2];
                }
            }
            else {
                for (j=0;j<(loop_length+2)*4;j++) {
                    proposed_model[j][0] = loop_model[j][0];
                    proposed_model[j][1] = loop_model[j][1];
                    proposed_model[j][2] = loop_model[j][2];
                }
            }                
        }       
                    
        for (j=0; j<4; j++) {
            for (k=0; k<3; k++) {
                moving_N_anchor[j][k] = loop_model[(loop_length + 1)*4 + j][k];
            }
        } 
    
        N_anchor_RMSD = calculate_RMSD2(moving_N_anchor, fixed_N_anchor_reordered, 1, 4);

        if (iteration_count >= 100) {
            return 0;
        } 
        
        sum_probs = 0;
        if (iteration_count > 0) {
            for (j=0; j<loop_length;j++) {
                sum_probs = sum_probs + probabilities[j];
            }
            if (sum_probs == 0) {
                return 0; 
            }
        }
           
        /*Set up next iteration*/
        if (i == loop_length * 2 - 2) {
            i = 0;
            iteration_count += 1;
        }
        else {
            i += 2;
        }
    }
        
    /*Insert the true N-anchor coordinates*/
    /*for (j=0; j<4; j++) {
        loop_model[(loop_length + 1) * 4 + j][0] = fixed_N_anchor_reordered[j][0];
        loop_model[(loop_length + 1) * 4 + j][1] = fixed_N_anchor_reordered[j][1];
        loop_model[(loop_length + 1) * 4 + j][2] = fixed_N_anchor_reordered[j][2];
    }*/

    return 1;
}

