/*------------------------------------------------------------------------------------------------*/
/*                                        MECHANO FUNCTIONS                                       */
/*------------------------------------------------------------------------------------------------*/

void write_structure_to_file_Cb(double **loop_coordinates, int model_no, FILE *output, char sequence[], 
                             char chain[], int loop_length, double RMSD, int res_numbers[], 
                             int last_atom_no, int direction, char decoy_code[])
{
    int i;
    int j;
    int ins_count = 0;
    double x;
    double y;
    double z;
    char ins_codes[] = {' ','A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S',
                        'T','U','V','W','X','Y','Z'};
    char directions[] = {'N', 'C'};
    
    if (RMSD == 0) {
        fprintf(output, "MODEL %s_%d direction = %c\n", decoy_code, model_no, directions[direction]);
    }
    else {
        fprintf(output, "MODEL %s_%d RMSD = %f direction = %c\n", 
                decoy_code, model_no, RMSD, directions[direction]);
    }
    
    for (i=0; i<(loop_length)*5; i++) {
        
        /* Glycine has no beta carbon */
        if (sequence[i/5] == 'G' && i%5 == 4) {
            continue;
        }

        if (res_numbers[i/5] == res_numbers[i/5-1]) {
            if (i%5 == 0) {
                ins_count++;
            }
        }
        else {
            ins_count = 0;
        }
        
        fprintf(output, "ATOM%7d  %s%s %s%4d%c   %8.3f%8.3f%8.3f  1.00  0.00\n", 
                last_atom_no+i+1, atom_type_strings[i%5], to_three_letter[aa_indices[sequence[i/5]]],
                chain, res_numbers[i/5], ins_codes[ins_count], loop_coordinates[i][0], 
                loop_coordinates[i][1], loop_coordinates[i][2]);
    }
    
    fprintf(output, "ENDMDL\n");
}

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

void anchor_extractor(char pdb_file[], char chain[], int start, int end, double N_anchor[][3], 
     char N_seq[4], double C_anchor[][3], char C_seq[4], double **native_loop, int loop_length,
     int info_array[2])
{
    /*This function extracts the N and C anchor coordinates and sequences, as well as the 
      native loop coordinates*/
      
    /* Returns 1 if the whole native structure is present, 0 if not*/
    
    char record_name[7];
    char residue_no_char[5];
    char last_atom_no_char[10];
    char atom_type[5];
    char x_char[9],y_char[9],z_char[9];
    char res_type[4];
    int residue_no;
    int N_count=0,C_count=0,native_count=0;
    FILE *f;
    f = fopen(pdb_file,"r");
    
    if (f != NULL){
        char line[100];
        while (fgets(line, sizeof line, f) != NULL){
            
            //Choose the residues from the correct chain
            strncpy(record_name, line, 6);
            if (line[21] == chain[0] && strcmp(record_name,"ATOM  ")==0){
                residue_no = atoi(strncpy(residue_no_char,line+22,4));
                
                //N ANCHOR
                if (residue_no >= start - 3 && residue_no <= start - 1){
                    strncpy(atom_type,line+12,4);
                    atom_type[4] = '\0';
                    
                    info_array[1] = atoi(strncpy(last_atom_no_char,line+6,5));;
                    
                    //Extract only backbone atoms
                    if (strcmp(atom_type," N  ")==0 || strcmp(atom_type," CA ")==0 || 
                        strcmp(atom_type," C  ")==0 || strcmp(atom_type," O  ")==0){
                        N_anchor[N_count][0] = atof(strncpy(x_char,line+30,8));
                        N_anchor[N_count][1] = atof(strncpy(y_char,line+38,8));
                        N_anchor[N_count][2] = atof(strncpy(z_char,line+46,8));
                        N_count++;
                        
                        //Get sequence
                        if (strcmp(atom_type," CA ")==0){
                            strncpy(res_type,line+17,3);
                            res_type[3] = '\0';
                            if (strcmp(res_type,"ALA")==0) N_seq[N_count/4] = 'A';
                            else if (strcmp(res_type,"ARG")==0) N_seq[N_count/4] = 'R';
                            else if (strcmp(res_type,"ASN")==0) N_seq[N_count/4] = 'N';
                            else if (strcmp(res_type,"ASP")==0) N_seq[N_count/4] = 'D';
                            else if (strcmp(res_type,"CYS")==0) N_seq[N_count/4] = 'C';
                            else if (strcmp(res_type,"GLU")==0) N_seq[N_count/4] = 'E';
                            else if (strcmp(res_type,"GLN")==0) N_seq[N_count/4] = 'Q';
                            else if (strcmp(res_type,"GLY")==0) N_seq[N_count/4] = 'G';
                            else if (strcmp(res_type,"HIS")==0) N_seq[N_count/4] = 'H';
                            else if (strcmp(res_type,"ILE")==0) N_seq[N_count/4] = 'I';
                            else if (strcmp(res_type,"LEU")==0) N_seq[N_count/4] = 'L';
                            else if (strcmp(res_type,"LYS")==0) N_seq[N_count/4] = 'K';
                            else if (strcmp(res_type,"MET")==0) N_seq[N_count/4] = 'M';
                            else if (strcmp(res_type,"PHE")==0) N_seq[N_count/4] = 'F';
                            else if (strcmp(res_type,"PRO")==0) N_seq[N_count/4] = 'P';
                            else if (strcmp(res_type,"SER")==0) N_seq[N_count/4] = 'S';
                            else if (strcmp(res_type,"THR")==0) N_seq[N_count/4] = 'T';
                            else if (strcmp(res_type,"TRP")==0) N_seq[N_count/4] = 'W';
                            else if (strcmp(res_type,"TYR")==0) N_seq[N_count/4] = 'Y';
                            else if (strcmp(res_type,"VAL")==0) N_seq[N_count/4] = 'V';
                            else {
                                printf("Error: unrecognised residue type (%s)\n",res_type);
                            }
                        }     
                    }
                }
                              
                //C ANCHOR
                else if (residue_no > end && residue_no <= end + 3){
                    strncpy(atom_type,line+12,4);
                    atom_type[4] = '\0';
                    
                    //Extract only backbone atoms
                    if (strcmp(atom_type," N  ")==0 || strcmp(atom_type," CA ")==0 || 
                        strcmp(atom_type," C  ")==0 || strcmp(atom_type," O  ")==0){
                        C_anchor[C_count][0] = atof(strncpy(x_char,line+30,8));
                        C_anchor[C_count][1] = atof(strncpy(y_char,line+38,8));
                        C_anchor[C_count][2] = atof(strncpy(z_char,line+46,8));
                        C_count++;
                        
                        //Get sequence
                        if (strcmp(atom_type," CA ")==0){
                            strncpy(res_type,line+17,3);
                            res_type[3] = '\0';
                            if (strcmp(res_type,"ALA")==0) C_seq[C_count/4] = 'A';
                            else if (strcmp(res_type,"ARG")==0) C_seq[C_count/4] = 'R';
                            else if (strcmp(res_type,"ASN")==0) C_seq[C_count/4] = 'N';
                            else if (strcmp(res_type,"ASP")==0) C_seq[C_count/4] = 'D';
                            else if (strcmp(res_type,"CYS")==0) C_seq[C_count/4] = 'C';
                            else if (strcmp(res_type,"GLU")==0) C_seq[C_count/4] = 'E';
                            else if (strcmp(res_type,"GLN")==0) C_seq[C_count/4] = 'Q';
                            else if (strcmp(res_type,"GLY")==0) C_seq[C_count/4] = 'G';
                            else if (strcmp(res_type,"HIS")==0) C_seq[C_count/4] = 'H';
                            else if (strcmp(res_type,"ILE")==0) C_seq[C_count/4] = 'I';
                            else if (strcmp(res_type,"LEU")==0) C_seq[C_count/4] = 'L';
                            else if (strcmp(res_type,"LYS")==0) C_seq[C_count/4] = 'K';
                            else if (strcmp(res_type,"MET")==0) C_seq[C_count/4] = 'M';
                            else if (strcmp(res_type,"PHE")==0) C_seq[C_count/4] = 'F';
                            else if (strcmp(res_type,"PRO")==0) C_seq[C_count/4] = 'P';
                            else if (strcmp(res_type,"SER")==0) C_seq[C_count/4] = 'S';
                            else if (strcmp(res_type,"THR")==0) C_seq[C_count/4] = 'T';
                            else if (strcmp(res_type,"TRP")==0) C_seq[C_count/4] = 'W';
                            else if (strcmp(res_type,"TYR")==0) C_seq[C_count/4] = 'Y';
                            else if (strcmp(res_type,"VAL")==0) C_seq[C_count/4] = 'V';
                            else {
                                printf("Error: unrecognised residue type (%s)\n",res_type);
                            }
                        }
                    }
                }
                
                //NATIVE LOOP
                else if (residue_no >= start && residue_no <= end){
                    strncpy(atom_type,line+12,4);
                    atom_type[4] = '\0';
                    
                    //Extract only backbone atoms
                    if (strcmp(atom_type," N  ")==0 || strcmp(atom_type," CA ")==0 || 
                        strcmp(atom_type," C  ")==0 || strcmp(atom_type," O  ")==0){
                        native_loop[native_count][0] = atof(strncpy(x_char,line+30,8));
                        native_loop[native_count][1] = atof(strncpy(y_char,line+38,8));
                        native_loop[native_count][2] = atof(strncpy(z_char,line+46,8));
                        native_count++;
                    }
                }
            }
        }
        N_seq[3] = '\0';
        C_seq[3] = '\0';
    }
    else {
        printf("Error: File %s does not exist\n",pdb_file);
        exit(0);
    }
    
    if (native_count == loop_length * 4) {
        info_array[0] = 1;
    }
    else {
        info_array[0] = 0;
    }
    
    fclose(f);
}

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

void count_atoms(char pdb_file[], char chain[], int start, int end, int atom_frequencies[4])
{
    char record_name[7];
    char residue_no_char[5];
    char atom_type[5];
    char x_char[9], y_char[9], z_char[9];
    int i;
    int residue_no;
    FILE *f;
    
    f = fopen(pdb_file,"r");
    
    for (i=0; i<4; i++) {
        atom_frequencies[i] = 0;
    }
    
    if (f != NULL){
        char line[100];
        while (fgets(line, sizeof line, f) != NULL){
            //Ignore the residues from the native loop structure and the anchor residue either side
            strncpy(record_name, line, 6);
            record_name[6] = '\0';
            if (line[21] == chain[0] && strcmp(record_name,"ATOM  ")==0) {
                residue_no = atoi(strncpy(residue_no_char,line+22,4));
                
                if (residue_no < start - 1 || residue_no > end + 1) {
                    strncpy(atom_type,line+12,4);
                    atom_type[4] = '\0';
                    
                    //Extract only backbone atoms
                    if (strcmp(atom_type," N  ")==0) { atom_frequencies[0]++; }
                    else if (strcmp(atom_type, " CA ")==0) { atom_frequencies[1]++; }
                    else if (strcmp(atom_type, " C  ")==0) { atom_frequencies[2]++; }
                    else if (strcmp(atom_type, " O  ")==0) { atom_frequencies[3]++; }
                }
            }
            else if (strcmp(record_name, "ATOM  ")==0) {
                strncpy(atom_type,line+12,4);
                atom_type[4] = '\0';
                
                //Extract only backbone atoms
                if (strcmp(atom_type," N  ")==0) { atom_frequencies[0]++; }
                else if (strcmp(atom_type, " CA ")==0) { atom_frequencies[1]++; }
                else if (strcmp(atom_type, " C  ")==0) { atom_frequencies[2]++; }
                else if (strcmp(atom_type, " O  ")==0) { atom_frequencies[3]++; }
            }
        }
    }
}
    
/*------------------------------------------------------------------------------------------------*/

void structure_extractor(char pdb_file[], char chain[], int start, int end, double **nitrogens, 
                         double **a_carbons, double **carbons, double **oxygens,
                         int atom_frequencies[4])
{
    char record_name[7];
    char residue_no_char[5];
    char atom_type[5];
    char x_char[9],y_char[9],z_char[9];
    int i;
    int residue_no;
    int N_count=0;
    int Ca_count=0;
    int C_count=0;
    int O_count=0;
    FILE *f;
    
    f = fopen(pdb_file, "r");
    
    if (f != NULL){
        char line[100];
        while (fgets(line, sizeof line, f) != NULL){
            //Ignore the residues from the native loop structure and the anchor residue either side
            strncpy(record_name, line, 6);
            record_name[6] = '\0';
            if (line[21] == chain[0] && strcmp(record_name,"ATOM  ")==0) {
                residue_no = atoi(strncpy(residue_no_char,line+22,4));
                
                if (residue_no < start - 1 || residue_no > end + 1) {
                    strncpy(atom_type,line+12,4);
                    atom_type[4] = '\0';
                    
                    //Extract only backbone atoms
                    if (strcmp(atom_type," N  ")==0) {
                        nitrogens[N_count][0] = atof(strncpy(x_char, line+30, 8));
                        nitrogens[N_count][1] = atof(strncpy(x_char, line+38, 8));
                        nitrogens[N_count][2] = atof(strncpy(x_char, line+46, 8));
                        N_count++;
                    }
                    else if (strcmp(atom_type, " CA ")==0) {
                        a_carbons[Ca_count][0] = atof(strncpy(x_char, line+30, 8));
                        a_carbons[Ca_count][1] = atof(strncpy(x_char, line+38, 8));
                        a_carbons[Ca_count][2] = atof(strncpy(x_char, line+46, 8));
                        Ca_count++;
                    }
                    else if (strcmp(atom_type, " C  ")==0) {
                        carbons[C_count][0] = atof(strncpy(x_char, line+30, 8));
                        carbons[C_count][1] = atof(strncpy(x_char, line+38, 8));
                        carbons[C_count][2] = atof(strncpy(x_char, line+46, 8));
                        C_count++;
                    }
                    else if (strcmp(atom_type, " O  ")==0) {
                        oxygens[O_count][0] = atof(strncpy(x_char, line+30, 8));
                        oxygens[O_count][1] = atof(strncpy(x_char, line+38, 8));
                        oxygens[O_count][2] = atof(strncpy(x_char, line+46, 8));
                        O_count++;
                    }
                }
            }
            else if (strcmp(record_name, "ATOM  ")==0) {
                strncpy(atom_type,line+12,4);
                atom_type[4] = '\0';
                
                //Extract only backbone atoms
                if (strcmp(atom_type," N  ")==0) {
                    nitrogens[N_count][0] = atof(strncpy(x_char, line+30, 8));
                    nitrogens[N_count][1] = atof(strncpy(x_char, line+38, 8));
                    nitrogens[N_count][2] = atof(strncpy(x_char, line+46, 8));
                    N_count++;
                }
                else if (strcmp(atom_type, " CA ")==0) {
                    a_carbons[Ca_count][0] = atof(strncpy(x_char, line+30, 8));
                    a_carbons[Ca_count][1] = atof(strncpy(x_char, line+38, 8));
                    a_carbons[Ca_count][2] = atof(strncpy(x_char, line+46, 8));
                    Ca_count++;
                }
                else if (strcmp(atom_type, " C  ")==0) {
                    carbons[C_count][0] = atof(strncpy(x_char, line+30, 8));
                    carbons[C_count][1] = atof(strncpy(x_char, line+38, 8));
                    carbons[C_count][2] = atof(strncpy(x_char, line+46, 8));
                    C_count++;
                }
                else if (strcmp(atom_type, " O  ")==0) {
                    oxygens[O_count][0] = atof(strncpy(x_char, line+30, 8));
                    oxygens[O_count][1] = atof(strncpy(x_char, line+38, 8));
                    oxygens[O_count][2] = atof(strncpy(x_char, line+46, 8));
                    O_count++;
                }
            }
        }
    }
}
                        
/*------------------------------------------------------------------------------------------------*/

int find_dihedral_match(int array_index, int angle, int angle_type)
{
    /*Find a dihedral angle that is possible given the other dihedral angle for that residue*/
    /*angle_type = 0 if the angle given is phi, 1 if psi*/
    /*extracts all possible angles, then chooses one at random*/
    /*if no match is present then the angle searched for is changed iteratively until one is found*/
    
    int extracted_dihedrals[dihedral_array_lengths[array_index]];
    int count = 0;
    int increased_angle;
    int decreased_angle;
    int i;
    int r;
    
    for (i = 0; i < dihedral_array_lengths[array_index] * 2; i+=2) {
        if (*(ptr_array[array_index] + i + angle_type) == angle) {
            extracted_dihedrals[count] = *(ptr_array[array_index] + i + (1-angle_type));
            count++;
        }
    }
    
    increased_angle = angle;
    decreased_angle = angle;
    
    while (count == 0) {
        increased_angle = increased_angle + 1;
        decreased_angle = decreased_angle - 1;
        
        for (i=0; i < dihedral_array_lengths[array_index] * 2; i+=2) {
            if (*(ptr_array[array_index] + i + angle_type) == increased_angle || 
                *(ptr_array[array_index] + i + angle_type) == decreased_angle) {
                     extracted_dihedrals[count] = *(ptr_array[array_index] + i + (1-angle_type));
                     count++;
            }
        }
    }
    
    r = rand() % count;
    
    return extracted_dihedrals[r];
}

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

double omega_angle_generator()
{
    int r;
    double omega;
    
    r = rand() % 10000;
    
    if (r < 59) {
        /*cis peptide bond*/
        omega = random_gaussian(1.220, 11.495);
    }
    else {
        /*trans peptide bond*/
        omega = random_gaussian(180.336, 5.930);
        
        if (omega > 180) {
            omega = omega - 360;
        }
    }
    
    return omega;
}

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

void dihedral_angle_generator_N(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[i] == 'A') { array_index =  0; }
        else if (seq_inc_anchors[i] == 'R') { array_index =  1; }
        else if (seq_inc_anchors[i] == 'N') { array_index =  2; }
        else if (seq_inc_anchors[i] == 'D') { array_index =  3; }
        else if (seq_inc_anchors[i] == 'C') { array_index =  4; }
        else if (seq_inc_anchors[i] == 'E') { array_index =  5; }
        else if (seq_inc_anchors[i] == 'Q') { array_index =  6; }
        else if (seq_inc_anchors[i] == 'G') { array_index =  7; }
        else if (seq_inc_anchors[i] == 'H') { array_index =  8; }
        else if (seq_inc_anchors[i] == 'I') { array_index =  9; }
        else if (seq_inc_anchors[i] == 'L') { array_index = 10; }
        else if (seq_inc_anchors[i] == 'K') { array_index = 11; }
        else if (seq_inc_anchors[i] == 'M') { array_index = 12; }
        else if (seq_inc_anchors[i] == 'F') { array_index = 13; }
        else if (seq_inc_anchors[i] == 'P') { array_index = 14; }
        else if (seq_inc_anchors[i] == 'S') { array_index = 15; }
        else if (seq_inc_anchors[i] == 'T') { array_index = 16; }
        else if (seq_inc_anchors[i] == 'W') { array_index = 17; }
        else if (seq_inc_anchors[i] == 'Y') { array_index = 18; }
        else if (seq_inc_anchors[i] == 'V') { array_index = 19; }
        else {
            printf("Error: unknown residue type: %c\n", seq_inc_anchors[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 == 0) {
            if (first_O_dihedral < 0) {
                dihedral_angles[0] = 180 + first_O_dihedral;
            }
            else if (first_O_dihedral > 0) {
                dihedral_angles[0] = first_O_dihedral - 180;
            }
            else {
                dihedral_angles[0] = 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 - 1) * 4 + 1] = omega;
            dihedral_angles[(i - 1) * 4 + 2] = (phi + (rand() + 1.0) / (RAND_MAX + 1.0));
            dihedral_angles[i * 4] = psi + (rand() + 1.0) / (RAND_MAX + 1.0);
            
            if (dihedral_angles[i * 4] < 0) {
                dihedral_angles[(i - 1) * 4 + 3] = 180 + dihedral_angles[i*4];
            }
            else if (dihedral_angles[i * 4] > 0) {
                dihedral_angles[(i - 1) * 4 + 3] = dihedral_angles[i*4] - 180;
            }
            else {
                dihedral_angles[(i - 1) * 4 + 3] = 180;
            }
        }
        
        /*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 == no_residues - 1) {
            phi = find_dihedral_match(array_index, (int)roundf(last_psi), 1);
            omega = omega_angle_generator();
            dihedral_angles[(i - 1) * 4 + 1] = omega;
            dihedral_angles[(i - 1) * 4 + 2] = (phi + (rand() + 1.0) / (RAND_MAX + 1.0));
            dihedral_angles[(i - 1) * 4 + 3] = last_O_dihedral;
        }
    }
    
    /*Convert to radians*/
    for (i=0;i<no_atoms;i++){
        dihedral_angles[i] = dihedral_angles[i] * pi / 180;
    }  
    
}

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

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

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

    i = no_residues-1;

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

    /*Calculate bond lengths from C anchor structure*/
    for (j=1; j<4; j++) {
        vector_subtract(C_anchor[j-1], C_anchor[j], bond);
        bond_lengths[(i-1)*4+j] = magnitude(bond);
    }  
}

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

void bond_angle_generator_N(char seq_inc_anchors[], double bond_angles[], double C_anchor[][3])
{
    double bond1[3], bond2[3];
    double cos_angle;
    int no_residues = strlen(seq_inc_anchors);
    int i, j;

    for (i=1; i<no_residues-1; i++) {

        if (seq_inc_anchors[i] == 'G') {
            if (seq_inc_anchors[i-1] == 'G') {
                bond_angles[(i-1)*4] = random_gaussian(116.4, 2.1) * pi / 180;
            }
            else {
                bond_angles[(i-1)*4] = random_gaussian(116.2, 2.0) * pi / 180;
            }
            bond_angles[(i-1)*4+1] = random_gaussian(120.6, 1.7) * pi / 180;
            bond_angles[(i-1)*4+2] = random_gaussian(112.5, 2.9) * pi / 180;
            bond_angles[(i-1)*4+3] = random_gaussian(120.8, 2.1) * pi / 180;
        }

        else if (seq_inc_anchors[i] == 'P') {
            if (seq_inc_anchors[i-1] == 'G') {
                bond_angles[(i-1)*4] = random_gaussian(118.2, 2.1) * pi / 180;
            }
            else {
                bond_angles[(i-1)*4] = random_gaussian(116.9, 2.0) * pi / 180;
            }
            bond_angles[(i-1)*4+1] = random_gaussian(116.9, 2.0) * pi / 180;
            bond_angles[(i-1)*4+2] = random_gaussian(111.8, 2.5) * pi / 180;
            bond_angles[(i-1)*4+3] = random_gaussian(120.8, 1.7) * pi / 180;
        }

        else {
            if (seq_inc_anchors[i-1] == 'G') {
                bond_angles[(i-1)*4] = random_gaussian(116.4, 2.1) * pi / 180;
            }
            else {
                bond_angles[(i-1)*4] = random_gaussian(116.2, 2.0) * pi / 180;
            }
            bond_angles[(i-1)*4+1] = random_gaussian(121.7, 1.8) * pi / 180;
            bond_angles[(i-1)*4+2] = random_gaussian(111.2, 2.8) * pi / 180;
            bond_angles[(i-1)*4+3] = random_gaussian(120.8, 1.7) * pi / 180;
        }
    }

    i = no_residues - 1;

    if (seq_inc_anchors[i] == 'G') {
        if (seq_inc_anchors[i-1] == 'G') {
            bond_angles[(i-1)*4] = random_gaussian(116.4, 2.1) * pi / 180;
        }
        else {
            bond_angles[(i-1)*4] = random_gaussian(116.2, 2.0) * pi / 180;
        }
        bond_angles[(i-1)*4+1] = random_gaussian(120.6, 1.7) * pi / 180;
    }

    else if (seq_inc_anchors[i] == 'P') {
        if (seq_inc_anchors[i-1] == 'G') {
            bond_angles[(i-1)*4] = random_gaussian(118.2, 2.1) * pi / 180;
        }
        else {
            bond_angles[(i-1)*4] = random_gaussian(116.9, 2.0) * pi / 180;
        }
        bond_angles[(i-1)*4+1] = random_gaussian(116.9, 2.0) * pi / 180;
    }

    else {
        if (seq_inc_anchors[i-1] == 'G') {
            bond_angles[(i-1)*4] = random_gaussian(116.4, 2.1) * pi / 180;
        }
        else {
            bond_angles[(i-1)*4] = random_gaussian(116.2, 2.0) * pi / 180;
        }
        bond_angles[(i-1)*4+1] = random_gaussian(121.7, 1.8) * pi / 180;
    }
    
    /*Calculate bond angles from C anchor structure*/
    for (i=1; i<3; i++) {
        vector_subtract(C_anchor[i], C_anchor[i-1], bond1);
        vector_subtract(C_anchor[i], C_anchor[i+1], bond2);

        cos_angle = dot_product(bond1, bond2) / (magnitude(bond1) * magnitude(bond2));
                                                              
        bond_angles[(no_residues-1)*4-(3-i)] = acos(cos_angle);
    }    
}
/*------------------------------------------------------------------------------------------------*/

void build_loop_N(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.*/
      
    int no_residues = strlen(seq_inc_anchors)-1; /*No. of residues to be added, includes C 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_N(seq_inc_anchors, bond_lengths, C_anchor);
 
    bond_angle_generator_N(seq_inc_anchors, bond_angles, C_anchor);

    dihedral_angle_generator_N(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];
            bond_angles[i] = FREADbond_angles[i];
        }
    }

    for (i=0; i<no_residues-1; i++) {
        if (bitstring[i] == '1') {
            dihedral_angles[i*4+2] = FREADdihedral_angles[i*4+2];
            dihedral_angles[i*4+3] = FREADdihedral_angles[i*4+3];
            dihedral_angles[i*4+4] = FREADdihedral_angles[i*4+4];
            dihedral_angles[i*4+5] = FREADdihedral_angles[i*4+5];
        }
    }
                             
    /*Fill in the first residue of the unclosed loop (the last N anchor residue)*/
    for (i=0;i<4;i++){
        loop_model[i][0] = N_anchor[i][0];
        loop_model[i][1] = N_anchor[i][1];
        loop_model[i][2] = N_anchor[i][2];
    }
    
    /*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) {
            /*N*/
            A_index = i;
            B_index = i + 1;
            C_index = i + 2;
        }
        else if ((i - 1) % 4 == 0) {
            /*C-alpha*/
            A_index = i;
            B_index = i + 1;
            C_index = i + 3;
        }
        else if ((i - 2) % 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;
        }
        
        /*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(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][j];       /*N of moving C anchor*/
        M02[j] = loop_model[(loop_length + 1) * 4 + 1][j];   /*Ca of moving C anchor*/
        M03[j] = loop_model[(loop_length + 1) * 4 + 2][j];   /*C of moving C 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);
    angle = atan2(sin_angle, cos_angle);   
    
    //printf("atan2 = %.3f, sin = %.3f, cos = %.3f\n",angle*180/pi, 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);
    }
} 

/*------------------------------------------------------------------------------------------------*/
/*Close loop using only angles not copied from database structure*/

int close_loop_N_first(double **loop_model, double fixed_C_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 moving_C_anchor[4][3];
    double N_anchor[4][3];
    double C_anchor_RMSD;
    double anchor_separation;
    double RMSD_cutoff = 1.0;
    double F1[3], F2[3], F3[3];
    double initial_C_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 + 1) * 4;
        }
        else {
            A_indices[i] = A_indices[i - 1] + 1;
        }
    }
    
    /*Find initial RMSD between fixed and moving C anchor atoms*/
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) {
            moving_C_anchor[i][j] = loop_model[(loop_length + 1)*4 + i][j];
            N_anchor[i][j] = loop_model[i][j];
            F1[j] = fixed_C_anchor[0][j];    /*N of fixed C anchor*/
            F2[j] = fixed_C_anchor[1][j];    /*Ca of fixed C anchor*/
            F3[j] = fixed_C_anchor[2][j];    /*C of fixed C anchor*/
        }
    } 
    
    C_anchor_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);
    initial_C_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);   
    
    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 (C_anchor_RMSD > RMSD_cutoff) {  
    
        //Allow residues next to MECHANO-modelled residues to change--------------------------------------------------------
        
        //If this is the first residue:
       /* if (A_indices[i]/4-2 < 0) {
            if (bitstring[A_indices[i]/4-1] == '1' && bitstring[A_indices[i]/4] == '1') {
                
                if (i == loop_length * 2 - 2) {
                    i = 0;
                    iteration_count += 1;
                }
                else {
                    i += 2;
                }
                
                continue;
            }
        }*/
        
        //If this is the last residue:
       /* if (A_indices[i]/4 > loop_length-1) {
            if (bitstring[A_indices[i]/4] == '1' && bitstring[A_indices[i]/4-1] == '1') {
                
                if (i == loop_length * 2 - 2) {
                    i = 0;
                    iteration_count += 1;
                }
                else {
                    i += 2;
                }
                
                continue;
            }
        }*/
        
        //If this is other residues:
        /*if (bitstring[A_indices[i]/4] == '1' && bitstring[A_indices[i]/4-1] == '1' && bitstring[A_indices[i]/4-2] == '1') {
            
            if (i == loop_length * 2 - 2) {
                i = 0;
                iteration_count += 1;
            }
            else {
                i += 2;
            }
            
            continue;
        }*/
        //------------------------------------------------------------------------------------------------------------------
        
        //Ignore residues modelled using FREAD------------------------------------------------------------------------------
        
        if (bitstring[A_indices[i]/4 - 1] == '1') {
            /*Set up next iteration*/
            if (i == loop_length * 2 - 2) {
                i = 0;
                iteration_count += 1;
            }
            else {
                i += 2;
            }
            continue;
        }
        
        //------------------------------------------------------------------------------------------------------------------
            
                
        CCD(proposed_model, i, F1, F2, F3, A_indices, loop_length);
        
        CCD(proposed_model, i+1, F1, F2, F3, A_indices, loop_length);
        
        /*Calculate previous dihedrals*/
        index = A_indices[i];
        old_phi = calculate_dihedral(loop_model[index-2], loop_model[index],
                                 loop_model[index+1], loop_model[index+2]);
                                 
        old_psi = calculate_dihedral(loop_model[index], loop_model[index+1],
                                 loop_model[index+2], loop_model[index+4]);
                                 
        /*Calculate new dihedrals*/
        new_phi = calculate_dihedral(proposed_model[index-2], proposed_model[index],
                                     proposed_model[index+1], proposed_model[index+2]);
                                 
        new_psi = 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[index/4] == 'A') { P_old = ALA_CCD_dist[P_old_index];
                                                    P_new = ALA_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'R') { P_old = ARG_CCD_dist[P_old_index];
                                                    P_new = ARG_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'N') { P_old = ASN_CCD_dist[P_old_index];
                                                    P_new = ASN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'D') { P_old = ASP_CCD_dist[P_old_index];
                                                    P_new = ASP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'C') { P_old = CYS_CCD_dist[P_old_index];
                                                    P_new = CYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'E') { P_old = GLU_CCD_dist[P_old_index];
                                                    P_new = GLU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'Q') { P_old = GLN_CCD_dist[P_old_index];
                                                    P_new = GLN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'G') { P_old = GLY_CCD_dist[P_old_index];
                                                    P_new = GLY_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'H') { P_old = HIS_CCD_dist[P_old_index];
                                                    P_new = HIS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'I') { P_old = ILE_CCD_dist[P_old_index];
                                                    P_new = ILE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'L') { P_old = LEU_CCD_dist[P_old_index];
                                                    P_new = LEU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'K') { P_old = LYS_CCD_dist[P_old_index];
                                                    P_new = LYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'M') { P_old = MET_CCD_dist[P_old_index];
                                                    P_new = MET_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'F') { P_old = PHE_CCD_dist[P_old_index];
                                                    P_new = PHE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'P') { P_old = PRO_CCD_dist[P_old_index];
                                                    P_new = PRO_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'S') { P_old = SER_CCD_dist[P_old_index];
                                                    P_new = SER_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'T') { P_old = THR_CCD_dist[P_old_index];
                                                    P_new = THR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'W') { P_old = TRP_CCD_dist[P_old_index];
                                                    P_new = TRP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'Y') { P_old = TYR_CCD_dist[P_old_index];
                                                    P_new = TYR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[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[index/4]);
            exit(0);
        }
        
        /*if (seq_inc_anchors[index/4] == 'P') {
            P_old = PRO_CCD_dist[P_old_index];
            P_new = PRO_CCD_dist[P_new_index];
        }
        else if (seq_inc_anchors[index/4] == 'G') { 
            P_old = GLY_CCD_dist[P_old_index];
            P_new = GLY_CCD_dist[P_new_index];
        }
        else { 
            P_old = OTR_CCD_dist[P_old_index];
            P_new = OTR_CCD_dist[P_new_index]; 
        }*/
        
        probabilities[i/2] = P_new;
        
        /*Decide whether to accept the new angles or not*/
        //if (i==0) {printf("old = %.15f, new = %.15f, phi = %.3f, psi = %.3f\n",P_old,P_new,new_phi,new_psi);}  
        if (P_new >= P_old) {          
            //if (i==0) {printf("    accepted automatically\n");}  
            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;
            //printf("%.10f\n",P);
            if (r <= P) {
                //if (i==0) {printf("    r = %f, accepted\n",r);}
                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 {
                //if (i==0) {printf("    r = %f, p = %f, rejected\n",r,P_new);}
                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_C_anchor[j][k] = loop_model[(loop_length + 1)*4 + j][k];
            }
        }
   
        C_anchor_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);
        
        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];
            }
            //printf("%.15f\n",sum_probs);
            if (sum_probs == 0) {
                //printf("         exiting loop closure...\n");
                return 0; 
            }
        }
           
        /*Set up next iteration*/
        if (i == loop_length * 2 - 2) {
            i = 0;
            iteration_count += 1;
        }
        else {
            i += 2;
        }
    }
        
    /*Insert the true C-anchor coordinates*/
    /*for (j=0; j<4; j++) {
        loop_model[(loop_length + 1) * 4 + j][0] = fixed_C_anchor[j][0];
        loop_model[(loop_length + 1) * 4 + j][1] = fixed_C_anchor[j][1];
        loop_model[(loop_length + 1) * 4 + j][2] = fixed_C_anchor[j][2];
    }*/
    
    return 1;
}

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

int close_loop_N(double **loop_model, double fixed_C_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 moving_C_anchor[4][3];
    double N_anchor[4][3];
    double C_anchor_RMSD;
    double anchor_separation;
    double RMSD_cutoff = 0.2;
    double F1[3], F2[3], F3[3];
    double initial_C_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 + 1) * 4;
        }
        else {
            A_indices[i] = A_indices[i - 1] + 1;
        }
    }
    
    /*Find initial RMSD between fixed and moving C anchor atoms*/
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) {
            moving_C_anchor[i][j] = loop_model[(loop_length + 1)*4 + i][j];
            N_anchor[i][j] = loop_model[i][j];
            F1[j] = fixed_C_anchor[0][j];    /*N of fixed C anchor*/
            F2[j] = fixed_C_anchor[1][j];    /*Ca of fixed C anchor*/
            F3[j] = fixed_C_anchor[2][j];    /*C of fixed C anchor*/
        }
    } 
    
    C_anchor_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);
    initial_C_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);   
    
    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 (C_anchor_RMSD > RMSD_cutoff) {  
    
        CCD(proposed_model, i, F1, F2, F3, A_indices, loop_length);
        
        CCD(proposed_model, i+1, F1, F2, F3, A_indices, loop_length);
        
        /*Calculate previous dihedrals*/
        index = A_indices[i];
        old_phi = calculate_dihedral(loop_model[index-2], loop_model[index],
                                 loop_model[index+1], loop_model[index+2]);
                                 
        old_psi = calculate_dihedral(loop_model[index], loop_model[index+1],
                                 loop_model[index+2], loop_model[index+4]);
                                 
        /*Calculate new dihedrals*/
        new_phi = calculate_dihedral(proposed_model[index-2], proposed_model[index],
                                     proposed_model[index+1], proposed_model[index+2]);
                                 
        new_psi = 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[index/4] == 'A') { P_old = ALA_CCD_dist[P_old_index];
                                                    P_new = ALA_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'R') { P_old = ARG_CCD_dist[P_old_index];
                                                    P_new = ARG_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'N') { P_old = ASN_CCD_dist[P_old_index];
                                                    P_new = ASN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'D') { P_old = ASP_CCD_dist[P_old_index];
                                                    P_new = ASP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'C') { P_old = CYS_CCD_dist[P_old_index];
                                                    P_new = CYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'E') { P_old = GLU_CCD_dist[P_old_index];
                                                    P_new = GLU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'Q') { P_old = GLN_CCD_dist[P_old_index];
                                                    P_new = GLN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'G') { P_old = GLY_CCD_dist[P_old_index];
                                                    P_new = GLY_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'H') { P_old = HIS_CCD_dist[P_old_index];
                                                    P_new = HIS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'I') { P_old = ILE_CCD_dist[P_old_index];
                                                    P_new = ILE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'L') { P_old = LEU_CCD_dist[P_old_index];
                                                    P_new = LEU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'K') { P_old = LYS_CCD_dist[P_old_index];
                                                    P_new = LYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'M') { P_old = MET_CCD_dist[P_old_index];
                                                    P_new = MET_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'F') { P_old = PHE_CCD_dist[P_old_index];
                                                    P_new = PHE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'P') { P_old = PRO_CCD_dist[P_old_index];
                                                    P_new = PRO_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'S') { P_old = SER_CCD_dist[P_old_index];
                                                    P_new = SER_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'T') { P_old = THR_CCD_dist[P_old_index];
                                                    P_new = THR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'W') { P_old = TRP_CCD_dist[P_old_index];
                                                    P_new = TRP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'Y') { P_old = TYR_CCD_dist[P_old_index];
                                                    P_new = TYR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[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[index/4]);
            exit(0);
        }
        
        /*if (seq_inc_anchors[index/4] == 'P') {
            P_old = PRO_CCD_dist[P_old_index];
            P_new = PRO_CCD_dist[P_new_index];
        }
        else if (seq_inc_anchors[index/4] == 'G') { 
            P_old = GLY_CCD_dist[P_old_index];
            P_new = GLY_CCD_dist[P_new_index];
        }
        else { 
            P_old = OTR_CCD_dist[P_old_index];
            P_new = OTR_CCD_dist[P_new_index]; 
        }*/
        
        probabilities[i/2] = P_new;
        
        /*Decide whether to accept the new angles or not*/
        //if (i==0) {printf("old = %.15f, new = %.15f, phi = %.3f, psi = %.3f\n",P_old,P_new,new_phi,new_psi);}  
        if (P_new >= P_old) {          
            //if (i==0) {printf("    accepted automatically\n");}  
            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;
            //printf("%.10f\n",P);
            if (r <= P) {
                //if (i==0) {printf("    r = %f, accepted\n",r);}
                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 {
                //if (i==0) {printf("    r = %f, p = %f, rejected\n",r,P_new);}
                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_C_anchor[j][k] = loop_model[(loop_length + 1)*4 + j][k];
            }
        }
   
        C_anchor_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);
        
        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];
            }
            //printf("%.15f\n",sum_probs);
            if (sum_probs == 0) {
                //printf("         exiting loop closure...\n");
                return 0; 
            }
        }
           
        /*Set up next iteration*/
        if (i == loop_length * 2 - 2) {
            i = 0;
            iteration_count += 1;
        }
        else {
            i += 2;
        }
    }
        
    /*Insert the true C-anchor coordinates*/
    for (j=0; j<4; j++) {
        loop_model[(loop_length + 1) * 4 + j][0] = fixed_C_anchor[j][0];
        loop_model[(loop_length + 1) * 4 + j][1] = fixed_C_anchor[j][1];
        loop_model[(loop_length + 1) * 4 + j][2] = fixed_C_anchor[j][2];
    }
    
    return 1;
}

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

/*Close loop using only angles not copied from database structure*/

int close_loop_N_limited(double **loop_model, double fixed_C_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 FREAD_modelled; //1 = yes, 0 = no
    int residue_no;
    double moving_C_anchor[4][3];
    double N_anchor[4][3];
    double C_anchor_RMSD;
    double anchor_separation;
    double RMSD_cutoff = 1.0;
    double F1[3], F2[3], F3[3];
    double initial_C_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 + 1) * 4;
        }
        else {
            A_indices[i] = A_indices[i - 1] + 1;
        }
    }
    
    /*Find initial RMSD between fixed and moving C anchor atoms*/
    for (i=0; i<4; i++) {
        for (j=0; j<3; j++) {
            moving_C_anchor[i][j] = loop_model[(loop_length + 1)*4 + i][j];
            N_anchor[i][j] = loop_model[i][j];
            F1[j] = fixed_C_anchor[0][j];    /*N of fixed C anchor*/
            F2[j] = fixed_C_anchor[1][j];    /*Ca of fixed C anchor*/
            F3[j] = fixed_C_anchor[2][j];    /*C of fixed C anchor*/
        }
    } 
    
    C_anchor_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);
    initial_C_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);   
    
    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 (C_anchor_RMSD > RMSD_cutoff) {  
        
        //Was residue modelled using FREAD?-----------------------------------------------------------------------------
        residue_no = (A_indices[i]/4 - 1);        
        if (bitstring[residue_no] == '1') {
            FREAD_modelled = 1;
        }
        else {
            FREAD_modelled = 0;
        }
        
        //------------------------------------------------------------------------------------------------------------------
            
                
        CCD(proposed_model, i, F1, F2, F3, A_indices, loop_length);
        
        CCD(proposed_model, i+1, F1, F2, F3, A_indices, loop_length);
        
        /*Calculate previous dihedrals*/
        index = A_indices[i];
        old_phi = calculate_dihedral(loop_model[index-2], loop_model[index],
                                 loop_model[index+1], loop_model[index+2]);
                                 
        old_psi = calculate_dihedral(loop_model[index], loop_model[index+1],
                                 loop_model[index+2], loop_model[index+4]);
                                 
        /*Calculate new dihedrals*/
        new_phi = calculate_dihedral(proposed_model[index-2], proposed_model[index],
                                     proposed_model[index+1], proposed_model[index+2]);
                                 
        new_psi = 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[index/4] == 'A') { P_old = ALA_CCD_dist[P_old_index];
                                                    P_new = ALA_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'R') { P_old = ARG_CCD_dist[P_old_index];
                                                    P_new = ARG_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'N') { P_old = ASN_CCD_dist[P_old_index];
                                                    P_new = ASN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'D') { P_old = ASP_CCD_dist[P_old_index];
                                                    P_new = ASP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'C') { P_old = CYS_CCD_dist[P_old_index];
                                                    P_new = CYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'E') { P_old = GLU_CCD_dist[P_old_index];
                                                    P_new = GLU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'Q') { P_old = GLN_CCD_dist[P_old_index];
                                                    P_new = GLN_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'G') { P_old = GLY_CCD_dist[P_old_index];
                                                    P_new = GLY_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'H') { P_old = HIS_CCD_dist[P_old_index];
                                                    P_new = HIS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'I') { P_old = ILE_CCD_dist[P_old_index];
                                                    P_new = ILE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'L') { P_old = LEU_CCD_dist[P_old_index];
                                                    P_new = LEU_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'K') { P_old = LYS_CCD_dist[P_old_index];
                                                    P_new = LYS_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'M') { P_old = MET_CCD_dist[P_old_index];
                                                    P_new = MET_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'F') { P_old = PHE_CCD_dist[P_old_index];
                                                    P_new = PHE_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'P') { P_old = PRO_CCD_dist[P_old_index];
                                                    P_new = PRO_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'S') { P_old = SER_CCD_dist[P_old_index];
                                                    P_new = SER_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'T') { P_old = THR_CCD_dist[P_old_index];
                                                    P_new = THR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'W') { P_old = TRP_CCD_dist[P_old_index];
                                                    P_new = TRP_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[index/4] == 'Y') { P_old = TYR_CCD_dist[P_old_index];
                                                    P_new = TYR_CCD_dist[P_new_index]; }
        else if (seq_inc_anchors[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[index/4]);
            exit(0);
        }
        
        probabilities[i/2] = P_new;
        
        /*Decide whether to accept the new angles or not*/
        //if (i==0) {printf("old = %.15f, new = %.15f, phi = %.3f, psi = %.3f\n",P_old,P_new,new_phi,new_psi);}  
        if (P_new >= P_old) {          
            //if (i==0) {printf("    accepted automatically\n");}  
            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;
            //printf("%.10f\n",P);
            if (r <= P) {
                //if (i==0) {printf("    r = %f, accepted\n",r);}
                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 {
                //if (i==0) {printf("    r = %f, p = %f, rejected\n",r,P_new);}
                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_C_anchor[j][k] = loop_model[(loop_length + 1)*4 + j][k];
            }
        }
   
        C_anchor_RMSD = calculate_RMSD2(moving_C_anchor, fixed_C_anchor, 0, 3);
        
        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];
            }
            //printf("%.15f\n",sum_probs);
            if (sum_probs == 0) {
                //printf("         exiting loop closure...\n");
                return 0; 
            }
        }
           
        /*Set up next iteration*/
        if (i == loop_length * 2 - 2) {
            i = 0;
            iteration_count += 1;
        }
        else {
            i += 2;
        }
    }
        
    /*Insert the true C-anchor coordinates*/
    for (j=0; j<4; j++) {
        loop_model[(loop_length + 1) * 4 + j][0] = fixed_C_anchor[j][0];
        loop_model[(loop_length + 1) * 4 + j][1] = fixed_C_anchor[j][1];
        loop_model[(loop_length + 1) * 4 + j][2] = fixed_C_anchor[j][2];
    }
    
    return 1;
}

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

int self_clash_check(double **loop_model, int loop_length)
{
    int i;
    int j;
    int k;
    double atom1[3];
    double atom2[3];
    double vector12[3];
    double distance;
    double allowed_distance;
    double radii[4] = {1.23, 1.32, 1.21, 1.07};
  
    for (i=0; i<(loop_length+1)*4; i++) {
        for (j=i+4; j<(loop_length+2)*4; j++) {
            
            for (k=0; k<3; k++) {
                atom1[k] = loop_model[i][k];
                atom2[k] = loop_model[j][k];
            }
            
            vector_subtract(atom1, atom2, vector12);
            
            distance = magnitude(vector12);
            
            allowed_distance = (radii[i%4] + radii[j%4]) / 2;
            
            if (distance < allowed_distance) {
                return 0;
            }
        }
    }
    
    return 1;
}

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

int clash_check(double **loop_model, int loop_length, double nitrogens[][3], double a_carbons[][3], 
                double carbons[][3], double oxygens[][3], int atom_frequencies[4])
{
    int i;
    int j;
    double sq_dist;
    double allowed_dist;
    double radii[4] = {1.23, 1.32, 1.21, 1.07};
    
    for (i=4; i<(loop_length+1)*4; i++) {
        /*Check Backbone Nitrogens*/
        for (j=0; j<atom_frequencies[0]; j++) {
            sq_dist = pow((loop_model[i][0] - nitrogens[j][0]),2) + 
                      pow((loop_model[i][1] - nitrogens[j][1]),2) + 
                      pow((loop_model[i][2] - nitrogens[j][2]),2);
                      
            allowed_dist = pow((radii[i%4] + 1.23),2);
            
            if (sq_dist < allowed_dist) {
                return 0;
            }
        }
        
        /*Check Alpha Carbons*/    
        for (j=0; j<atom_frequencies[1]; j++) {
            sq_dist = pow((loop_model[i][0] - a_carbons[j][0]),2) + 
                      pow((loop_model[i][1] - a_carbons[j][1]),2) + 
                      pow((loop_model[i][2] - a_carbons[j][2]),2);
                      
            allowed_dist = pow((radii[i%4] + 1.32),2);
            
            if (sq_dist < allowed_dist) {
                return 0;
            }
        }
        
        /*Check Carbonyl Carbons*/
        for (j=0; j<atom_frequencies[2]; j++) {
            sq_dist = pow((loop_model[i][0] - carbons[j][0]),2) + 
                      pow((loop_model[i][1] - carbons[j][1]),2) + 
                      pow((loop_model[i][2] - carbons[j][2]),2);
                      
            allowed_dist = pow((radii[i%4] + 1.21),2);
            
            if (sq_dist < allowed_dist) {
                return 0;
            }
        }
        
        /*Check Oxygens*/
        for (j=0; j<atom_frequencies[3]; j++) {
            sq_dist = pow((loop_model[i][0] - oxygens[j][0]),2) + 
                      pow((loop_model[i][1] - oxygens[j][1]),2) + 
                      pow((loop_model[i][2] - oxygens[j][2]),2);
                      
            allowed_dist = pow((radii[i%4] + 1.07),2);
            
            if (sq_dist < allowed_dist) {
                return 0;
            }
        }
    }
    
    return 1;
}

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

void build_Cb_atoms(double **final_model, double **final_model_Cb, int loop_length, char sequence[])
{
    int residue_no;
    int i;
    double N[3];
    double Ca[3];
    double C[3];
    double Cb[3];
    double CaN[3];
    double CaC[3];
    double norm1[3];
    double norm2[3];
    double CaN_mag;
    double CaC_mag;
    double norm1_mag;
    double norm2_mag;
    double CaN_unit[3];
    double CaC_unit[3];
    double norm1_unit[3];
    double norm2_unit[3];
    double n_factor;
    double norm1_factor;
    double norm2_factor;
    double cb_dist;
    double Cb_mag;
    
    for (residue_no=0; residue_no<loop_length; residue_no++) {
        if (sequence[residue_no] == 'G') {
            for (i=0; i<3; i++) {
                /*Copy across atoms*/
                final_model_Cb[residue_no*5][i]   = final_model[residue_no*4][i];    //N
                final_model_Cb[residue_no*5+1][i] = final_model[residue_no*4+1][i];  //Ca
                final_model_Cb[residue_no*5+2][i] = final_model[residue_no*4+2][i];  //C
                final_model_Cb[residue_no*5+3][i] = final_model[residue_no*4+3][i];  //O
                final_model_Cb[residue_no*5+4][i] = 0;  //Cb, doesn't exist for glycine
            }
        }
        else {
            for (i=0; i<3; i++) {
                N[i]  = final_model[residue_no*4][i];
                Ca[i] = final_model[residue_no*4+1][i];
                C[i]  = final_model[residue_no*4+2][i];
            }
            
            vector_subtract(Ca, N, CaN);
            CaN_mag = magnitude(CaN);
            for (i=0; i<3; i++) {
                CaN_unit[i] = CaN[i] / CaN_mag;
            }
            
            vector_subtract(Ca, C, CaC);
            CaC_mag = magnitude(CaC);
            for (i=0; i<3; i++) {
                CaC_unit[i] = CaC[i] / CaC_mag;
            }
            
            cross_product(CaN_unit, CaC_unit, norm1);
            norm1_mag = magnitude(norm1);
            for (i=0; i<3; i++) {
                norm1_unit[i] = norm1[i] / norm1_mag;
            }
            
            cross_product(CaN_unit, norm1_unit, norm2);
            norm2_mag = magnitude(norm2);
            for (i=0; i<3; i++) {
                norm2_unit[i] = norm2[i] / norm2_mag;
            }
            
            if (sequence[residue_no] == 'A')      { n_factor     = -0.34777316024597776;
                                                    norm1_factor =  0.7887568908818895;
                                                    norm2_factor =  0.5056639074261424;
                                                    cb_dist      =  1.5253623347586374;}
            else if (sequence[residue_no] == 'R') { n_factor     = -0.3507360024977345;
                                                    norm1_factor =  0.7883874912065479;
                                                    norm2_factor =  0.5030476890981802;
                                                    cb_dist      =  1.5309172528392483;}
            else if (sequence[residue_no] == 'N') { n_factor     = -0.34954774900564267;
                                                    norm1_factor =  0.7821300410199126;
                                                    norm2_factor =  0.5126836623595018;
                                                    cb_dist      =  1.5321433083090938;}
            else if (sequence[residue_no] == 'D') { n_factor     = -0.3499688809989311;
                                                    norm1_factor =  0.786800521994171;
                                                    norm2_factor =  0.5054582115824542;
                                                    cb_dist      =  1.5318055409680786;}
            else if (sequence[residue_no] == 'C') { n_factor     = -0.3494607343363032;
                                                    norm1_factor =  0.7915552932428591;
                                                    norm2_factor =  0.49818968247044454;
                                                    cb_dist      =  1.5296357057135819;}
            else if (sequence[residue_no] == 'E') { n_factor     = -0.35008028783676887;
                                                    norm1_factor =  0.7871416517612738;
                                                    norm2_factor =  0.5055332413409879;
                                                    cb_dist      =  1.5316362178823235;}
            else if (sequence[residue_no] == 'Q') { n_factor     = -0.3504418980727456;
                                                    norm1_factor =  0.7878986800263459;
                                                    norm2_factor =  0.5039157935682486;
                                                    cb_dist      =  1.5307608842551144;}
            else if (sequence[residue_no] == 'H') { n_factor     = -0.3478745720232258;
                                                    norm1_factor =  0.7887577012128891;
                                                    norm2_factor =  0.5038179368915092;
                                                    cb_dist      =  1.53204409902269;}
            else if (sequence[residue_no] == 'I') { n_factor     = -0.3631729792461987;
                                                    norm1_factor =  0.7820277368466476;
                                                    norm2_factor =  0.5039364531320552;
                                                    cb_dist      =  1.5451861881417324;}
            else if (sequence[residue_no] == 'L') { n_factor     = -0.3492731366881571;
                                                    norm1_factor =  0.7911706392471316;
                                                    norm2_factor =  0.49990908813589296;
                                                    cb_dist      =  1.5313231292876908;}
            else if (sequence[residue_no] == 'K') { n_factor     = -0.35014881263058767;
                                                    norm1_factor =  0.7886286378593781;
                                                    norm2_factor =  0.5033890381484109;
                                                    cb_dist      =  1.5311461087301732;}
            else if (sequence[residue_no] == 'M') { n_factor     = -0.3507077306077249;
                                                    norm1_factor =  0.7891140701626578;
                                                    norm2_factor =  0.501435278701145;
                                                    cb_dist      =  1.530863273877535;}
            else if (sequence[residue_no] == 'F') { n_factor     = -0.3504095675472226;
                                                    norm1_factor =  0.7891263129068525;
                                                    norm2_factor =  0.5016236698213021;
                                                    cb_dist      =  1.5339592272837324;}
            else if (sequence[residue_no] == 'P') { n_factor     = -0.22910526303159487;
                                                    norm1_factor =  0.8441480630546168;
                                                    norm2_factor =  0.4831096753232061;
                                                    cb_dist      =  1.5324482100664505;}
            else if (sequence[residue_no] == 'S') { n_factor     = -0.3490800523345143;
                                                    norm1_factor =  0.790925255695327;
                                                    norm2_factor =  0.5000881622552047;
                                                    cb_dist      =  1.5298749212857365;}
            else if (sequence[residue_no] == 'T') { n_factor     = -0.3608824783469435;
                                                    norm1_factor =  0.787447497345662;
                                                    norm2_factor =  0.4970810557771811;
                                                    cb_dist      =  1.5401087763522436;}
            else if (sequence[residue_no] == 'W') { n_factor     = -0.3496482864883515;
                                                    norm1_factor =  0.7890753847881438;
                                                    norm2_factor =  0.5023111623195795;
                                                    cb_dist      =  1.5332689020009584;}
            else if (sequence[residue_no] == 'Y') { n_factor     = -0.3497958790092702;
                                                    norm1_factor =  0.7898336280013343;
                                                    norm2_factor =  0.5011455140152965;
                                                    cb_dist      =  1.5342036619521204;}
            else if (sequence[residue_no] == 'V') { n_factor     = -0.3648163340818305;
                                                    norm1_factor =  0.7822145583392209;
                                                    norm2_factor =  0.5026404428436114;
                                                    cb_dist      =  1.5453102405120656;}
                                                    
            for (i=0; i<3; i++) {
                Cb[i] = (CaN_unit[i]*n_factor) + (norm1_unit[i]*norm1_factor) + (norm2_unit[i]*norm2_factor);
            }
            Cb_mag = magnitude(Cb);
            for (i=0; i<3; i++) {
                Cb[i] = Cb[i] / Cb_mag * cb_dist;
            }
            
            for (i=0; i<3; i++) {
                /*Copy across atoms*/
                final_model_Cb[residue_no*5][i]   = final_model[residue_no*4][i];    //N
                final_model_Cb[residue_no*5+1][i] = final_model[residue_no*4+1][i];  //Ca
                final_model_Cb[residue_no*5+2][i] = final_model[residue_no*4+2][i];  //C
                final_model_Cb[residue_no*5+3][i] = final_model[residue_no*4+3][i];  //O
                final_model_Cb[residue_no*5+4][i] = Cb[i] + final_model[residue_no*4+1][i];  
            }
        }
    }        
} 
