/* EGAD: read_forcefield.cpp Navin Pokala and Tracy Handel Dept. of Molecular and Cell Biology University of California, Berkeley Copyright (C) 2003 Regents of the University of California GNU Public License Aug 12 2003 Absolutely no warranties are made or are implied with the use of this program or its parts. This file contains functions for reading forcefield parameter, residue topology, and rotamer library files. Some energy function invariants are also precalculated. */ #include "read_forcefield.h" ATOMRESPARAM PSEUDO_SIDECHAIN_ATOMRESPARAM, BKBN_N_ATOMRESPARAM,BKBN_CA_ATOMRESPARAM, BKBN_C_ATOMRESPARAM,BKBN_O_ATOMRESPARAM, BKBN_HA_ATOMRESPARAM; ATOMRESPARAM *CB_ATOMRESPARAM, *H_ATOMRESPARAM, *HA_ATOMRESPARAM, *C_ATOMRESPARAM, *O_ATOMRESPARAM, *N_ATOMRESPARAM, *CA_ATOMRESPARAM, *U_ATOMRESPARAM; ATOMRESPARAM *BHY_ATOMRESPARAM=NULL; ATOMRESPARAM PSEB; /* pseudoatom sidechain for Born radii precalculations */ RESPARAM *MEA; /* pseudoresidue for identifying core residues */ RESPARAM *PRO_RESPARAM, *GLY_RESPARAM; RESPARAM *N_TERM_RESPARAM, *C_TERM_RESPARAM; double NONBOND_FACTOR_1_4; double **VDW_WELLDEPTH=NULL; double **VDW_SIGMA6_14=NULL; double **VDW_WELLDEPTH_14=NULL; double **VDW_SIGMA6_SCALE=NULL; double **VDW_LINEAR_SLOPE=NULL; double **VDW_LINEAR_INTERCEPT=NULL; double **VDW_LINEAR_SWITCHPOINT=NULL; double **VDW_SIGMA=NULL; double **CHARGE_PRODUCT=NULL; double **VDW_CONTACT_DISTANCE=NULL; double WATER_DIELECTRIC = 80.0; double ONE_OVER_WATER_DIELECTRIC = 0.0125; double HBOND_DIELECTRIC = 0.0; double INTERNAL_DIELECTRIC = 0.0; double ONE_OVER_INTERNAL_DIELECTRIC = 0.125; double BORN_TAU = -0.1125; double COULOMB_CONST_OVER_INTERNAL_DIELECTRIC = 41.50795; double COUL_CONST_BORN_TAU = -37.357155; double HALF_COUL_CONST_BORN_TAU = -18.6785775; double COULOMB_CONST_HBOND_INTERNAL_DIELECTRIC; // used for hbonding; kc*(1/hbond_diel - 1/protein_diel) double CORE_INTERFACIAL_BORN_RADIUS_CUTOFF = 3.0; double INTERFACIAL_SURFACE_SASA_CUTOFF= 25; /* defaults based on survey of non-redundant pdb subset; avg/median of gaussian */ double FRACTION_HYDROPHOBIC_SASA_CUTOFF = 0.61; /* 1.0 std dev out */ double TRANSFER_FREE_ENERGY_DENSITY_CUTOFF = -0.000318; /* 1.0 std dev out */ int MAX_RESIMERS; int MAX_ROTAMERS; int MAX_RES_SIZE; double VAR_FIX_BORN_SCALE; double BORN_P1; double BORN_P2; double BORN_P3; double BORN_P4; double BORN_P5; double ONE_OVER_BORN_P5; double PI_TIMES_BORN_P5; double BORN_LAMBDA; double ONE_OVER_BORN_LAMBDA; double VDW_RADII_SCALE=0.0; double **UNITED_RADII_SUM_SQUARED=NULL; BACKBONE DUMMY_BKBN, PROLINE_BKBN; double HBOND_PROBE_BONDLENGTH = 0.1; double HBOND_WELLDEPTH = 236957.413356; double HBOND_REPULSIVE = 39.2; double HBOND_ATTRACTIVE = 6; int BHY_CLASS; double DONOR_UNSAT_HBOND_REFSTATE = ENDFLAG; double H_UNSAT_HBOND_REFSTATE = ENDFLAG; double ACCEPTOR_UNSAT_HBOND_REFSTATE = ENDFLAG; double ENERGY_ERROR = ENDFLAG; /* generate and calculate CHARGE_PRODUCT matrix for all possible atompairs called by read_forcefield.cpp: read_forcefield */ void calculate_charge_products(RESPARAM *resparam) { int i,j,q,n,k,atom_classes; extern double PH, **CHARGE_PRODUCT; double P_prot,P_unprot; extern int CHARGES_PH_INDEPENDENT_FLAG; atom_classes=0; i=1; while(strcmp(resparam[i].residuetype,"zzz")!=0) { for(j=4;j<=(resparam[i].numberofAtoms+3);++j) { ++atom_classes; resparam[i].atom[j].charge_class = atom_classes; } ++i; } if(CHARGE_PRODUCT==NULL) { CHARGE_PRODUCT = (double **)calloc(atom_classes+5, sizeof(double *)); for(i=0;i<=atom_classes;++i) CHARGE_PRODUCT[i] = (double *)calloc((atom_classes+5), sizeof(double)); CHARGE_PRODUCT[0][0] = 0; } if(CHARGES_PH_INDEPENDENT_FLAG==0 && PH!=ENDFLAG) { i=1; while(strcmp(resparam[i].residuetype,"zzz")!=0) { if(resparam[i].pKa != 20 && abs(resparam[i].residuecode)<100 && resparam[i].protonation_flag==1) { /* i is the protonated form; j is deprotonated */ j=1; while(resparam[i].residuecode != -resparam[j].residuecode) ++j; P_prot = exp(-(PH - resparam[i].pKa)*LN10); P_prot = P_prot/(1.0 + P_prot); /* probability of the protonated form */ P_unprot = 1.0 - P_prot; /* probability of the unprotonated form */ for(q=1;qvolume = (4.0/3.0)*PI*pow(PSEB.atom_ptr->united_radius, 3); PSEB.BORN_P4_volume = BORN_P4*PSEB.atom_ptr->volume; n=1; while(strcmp(inputresparam[n].residuetype, "zzz")!=0) { resparam = &inputresparam[n]; for(j=4;j<=(resparam->numberofAtoms+3);++j) resparam->atom[j].volume_bonded = resparam->atom[j].atom_ptr->volume; chi = (double *)calloc(resparam->rotamerlib_ptr->numChi+2, sizeof(double)); side_pdb = (pdbATOM *)calloc(MAX_RES_SIZE, sizeof(pdbATOM)); all_pdb = (mini_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(mini_pdbATOM)); side.atom = (micro_pdbATOM *)calloc(MAX_RES_SIZE, sizeof(micro_pdbATOM)); flag=0; if(resparam->residuetype[1]=='T' && resparam->residuetype[2]=='E' /* strstr(resparam->residuetype, "TE")!=0 */ ) /* terminii */ flag=1; bkbn = DUMMY_BKBN; if(strcmp(resparam->residuetype,"PRO")==0) bkbn = PROLINE_BKBN; chi[0]=0; for(q=1;q<=resparam->rotamerlib_ptr->numChi;++q) chi[q] = 180.0; chi[resparam->rotamerlib_ptr->numChi+1]=0; if(resparam->ligand_flag==1) { if(resparam->rotamerlib_ptr->numChi==1) { chi[1]=1; } else for(q=1;q<=12;++q) chi[q] = STATIONARY_LIGAND_TRANSFORM_MATRIX[q]; } q=1; if(flag==0) { build_a_SideChain(&side, &bkbn,resparam,&q,chi); make_side_pdbATOM(side_pdb,&side); } else strcpy(side_pdb[1].residuetype,"END"); q=4; all_pdb[q].coord = bkbn.N; all_pdb[q].atom_ptr = &resparam->atom[q]; all_pdb[q].sasa=0; ++q; all_pdb[q].sasa=0; all_pdb[q].coord = bkbn.H; all_pdb[q].atom_ptr = &resparam->atom[q]; ++q; all_pdb[q].sasa=0; all_pdb[q].coord = bkbn.CA; all_pdb[q].atom_ptr = &resparam->atom[q]; ++q; all_pdb[q].sasa=0; all_pdb[q].coord = bkbn.HA; all_pdb[q].atom_ptr = &resparam->atom[q]; ++q; all_pdb[q].sasa=0; all_pdb[q].coord = bkbn.CB; all_pdb[q].atom_ptr = &resparam->atom[q]; ++q; i=1; while(strcmp(side_pdb[i].residuetype,"END")!=0) { all_pdb[q].coord = side_pdb[i].coord; all_pdb[q].atom_ptr = side_pdb[i].atom_ptr; all_pdb[q].sasa = 0; ++q; ++i; } all_pdb[q].sasa=0; all_pdb[q].coord = bkbn.C; all_pdb[q].atom_ptr = &resparam->atom[q]; ++q; all_pdb[q].sasa=0; all_pdb[q].coord = bkbn.O; all_pdb[q].atom_ptr = &resparam->atom[q]; ++q; all_pdb[q].sasa=ENDFLAG; all_pdb[q].seq_position=ENDFLAG; r=0; bondstate=0; for(i=4;i<=(resparam->numberofAtoms+3);++i) if(resparam->atom[i].atom_ptr->united_radius>0) /* only non-zero radius atoms */ for(j=(i+1);j<=(resparam->numberofAtoms+3);++j) if(resparam->atom[j].atom_ptr->united_radius>0) /* only non-zero radius atoms */ { /* consider here only atoms that are within <=2 bonds or within the same aromatic group */ if(resparam->atom[i].other_info>0 && resparam->atom[j].other_info>0 && resparam->atom[i].other_info == resparam->atom[j].other_info) /* for fixed rings in a sidechain */ { r = Distance(all_pdb[i].coord, all_pdb[j].coord); if(resparam->atom[i].contactAtomNumber==j || resparam->atom[j].contactAtomNumber==i || resparam->atom[i].ringContactAtom==j || resparam->atom[j].ringContactAtom==i) bondstate=1; } else if(strcmp(resparam->residuetype, "NTE")==0 && /* special case N-term HN-N */ ((strcmp(resparam->atom[i].atomname, "N")==0 && strcmp(resparam->atom[j].atomtype, "HN")==0) || (strcmp(resparam->atom[j].atomname, "N")==0 && strcmp(resparam->atom[i].atomtype, "HN")==0))) { r = 1.01; bondstate=1; } else if(strcmp(resparam->residuetype, "CTE")==0 && /* special case C-term C-O */ ((strcmp(resparam->atom[i].atomname, "C")==0 && strncmp(resparam->atom[j].atomname, "O", 1)==0) || (strcmp(resparam->atom[j].atomname, "C")==0 && strncmp(resparam->atom[i].atomname, "O", 1)==0))) { r = 1.229; bondstate=1; } else if(strcmp(resparam->residuetype, "NTE")==0 && /* special case for the N-terminal group HN-CA interaction */ ((strcmp(resparam->atom[i].atomname, "CA")==0 && strcmp(resparam->atom[j].atomtype, "HN")==0) || (strcmp(resparam->atom[j].atomname, "CA")==0 && strcmp(resparam->atom[i].atomtype, "HN")==0))) { bondstate=2; } else /* special case for the N-terminal HN-HN interaction */ if(strcmp(resparam->residuetype, "NTE")==0 && ((strcmp(resparam->atom[i].atomtype, "HN")==0 && strcmp(resparam->atom[j].atomtype, "HN")==0) && (resparam->atom[i].atomnumber!=resparam->atom[j].atomnumber))) { bondstate=2; } else if(strcmp(resparam->residuetype, "CTE")==0 && /* C-term CA-O */ ((strcmp(resparam->atom[i].atomname, "CA")==0 || strncmp(resparam->atom[j].atomname, "O", 1)==0) || (strncmp(resparam->atom[i].atomtype, "O", 1)==0 || strcmp(resparam->atom[j].atomname, "CA")==0))) { bondstate=2; } else if(strcmp(resparam->residuetype, "CTE")==0 && /* C-term O-O */ ((strncmp(resparam->atom[i].atomname, "O", 1)==0 && strncmp(resparam->atom[j].atomtype, "O", 1)==0) && (resparam->atom[i].atomnumber!=resparam->atom[j].atomnumber))) { bondstate=2; } else if(resparam->atom[i].contactAtomNumber==j) { r=resparam->atom[i].bondlength; bondstate=1; } else if(resparam->atom[j].contactAtomNumber==i) { r=resparam->atom[j].bondlength; bondstate=1; } else if(resparam->atom[i].ringContactAtom==j || resparam->atom[j].ringContactAtom==i || resparam->atom[i].ringTwoBondAtom==j || resparam->atom[j].ringTwoBondAtom==i) { r = Distance(all_pdb[i].coord, all_pdb[j].coord); if(resparam->atom[i].ringContactAtom==j || resparam->atom[j].ringContactAtom==i) bondstate=1; } if(r!=0 && bondstate==1) { if(resparam->atom[i].charge!=0 && strcmp(resparam->atom[j].atomname, "HU")!=0) { h = volume_h(resparam->atom[i].atom_ptr->united_radius, resparam->atom[j].atom_ptr->united_radius, r); resparam->atom[i].volume_bonded -= PI_over3*h*h*(3*resparam->atom[i].atom_ptr->united_radius - h); } if(resparam->atom[j].charge!=0 && strcmp(resparam->atom[i].atomname, "HU")!=0) { h = volume_h(resparam->atom[j].atom_ptr->united_radius, resparam->atom[i].atom_ptr->united_radius, r); resparam->atom[j].volume_bonded -= PI_over3*h*h*(3*resparam->atom[j].atom_ptr->united_radius - h); } r=0; bondstate=0; } } free_memory(chi); free_memory(side_pdb); free_memory(all_pdb); free_memory(side.atom); ++n; } i=1; while(strcmp(inputresparam[i].residuetype, "zzz")!=0) { for(j=4;j<=(inputresparam[i].numberofAtoms+3);++j) { if(inputresparam[i].atom[j].volume_bonded <= 0) inputresparam[i].atom[j].volume_bonded = 0; } ++i; } } /* This function reads in a rotamer library from rotamerfilename into an allocated array of ROTAMERLIB Also reads in ligand "rotamers" called by read_forcefield.cpp: read_forcefield and ligand_stuff.cpp: read_ligand_data.cpp */ void read_rotamers(char *rotamerfilename, ROTAMERLIB *rotamerlib) { FILE *input, *ligand_coords_file_ptr; char *line, doh[6], *ligand_coord_filename, **word, *dummystring; int i,q,numChi,r,first, n, m, p, d, num_library_rotamers, library_spread_ratio; int nn, dd, version_flag,read_coordinates_flag, coord_file_flag, ligand_flag; double duh[10], denom; ROTAMER dummyrotamer; extern int LOTS_OF_ROTAMERLETS_FLAG; long fp; line = (char *)calloc(MAXLINE, sizeof(char)); ligand_coord_filename = (char *)calloc(MAXLINE, sizeof(char)); dummystring = (char *)calloc(MAXLINE, sizeof(char)); first=1; input = fopen_file(rotamerfilename,"r"); /* decide if the file is the old or new format */ version_flag=1; fgets(line,MAXLINE,input); while(strncmp(line,"START",5)!=0) { if(strstr(line,"VERSION 2")!=0) version_flag=2; fgets(line,MAXLINE,input); } if(strstr(line,"VERSION 2")!=0) version_flag=2; fgets(line,MAXLINE,input); r=1; while(strncmp(line,"STOP",4)!=0) { read_coordinates_flag=0; ligand_flag=0; rotamerlib[r].ligand_base_coords=NULL; rotamerlib[r].generate_ligamers_flag = 0; if(version_flag==1) { sscanf(line,"%s %d %d %d",rotamerlib[r].residuetype,&num_library_rotamers, &rotamerlib[r].number_of_rotamers, &numChi); library_spread_ratio = rotamerlib[r].number_of_rotamers/num_library_rotamers; } else /* version 2 rotamer library */ { num_library_rotamers=0; library_spread_ratio=0; if(word_count(line)==3 && strstr(line,"COORD")==0 && strstr(line,"coord")==0 && strstr(line,"GENER")==0 && strstr(line,"gener")==0) { sscanf(line,"%s %d %d",rotamerlib[r].residuetype, &numChi, &library_spread_ratio); } else { sscanf(line,"%s %s",rotamerlib[r].residuetype, dummystring); if(dummystring[0] > 58) /* ie : not a number, but a letter */ { convert_string_to_all_caps(dummystring); /* coords defined */ if(strstr(dummystring,"COORD")!=0) { read_coordinates_flag=1; coord_file_flag=0; ligand_flag=1; if(strstr(dummystring,"FILE")!=0) /* listed in seperate file */ { sscanf(line,"%s %s %s",rotamerlib[r].residuetype, dummystring, ligand_coord_filename); ligand_coords_file_ptr = fopen_file(ligand_coord_filename,"r"); coord_file_flag=1; } else /* for coord files: */ { /* pdb lines 1*/ ligand_coords_file_ptr = input; /* END */ } /* pdb lines 2*/ /* END */ /* loop to count "rotamers" */ /* STOP */ fp=ftell(ligand_coords_file_ptr); fseek(ligand_coords_file_ptr, fp, 0); fgets(line,MAXLINE,ligand_coords_file_ptr); /* line after header line */ sscanf(line,"%s",dummystring); while(strcmp(dummystring,"STOP")!=0) { while(strncmp(line,"END",3)!=0) /* move past coords for this "rotamer" */ { fgets(line,MAXLINE,ligand_coords_file_ptr); } ++num_library_rotamers; fgets(line,MAXLINE,ligand_coords_file_ptr); sscanf(line,"%s",dummystring); } fseek(ligand_coords_file_ptr, fp, 0); /* reset file pointer */ /* allocate memory and set some variables etc */ rotamerlib[r].number_of_rotamers=num_library_rotamers; rotamerlib[r].numChi=1; numChi=1; rotamerlib[r].rotamer = (ROTAMER *)calloc((rotamerlib[r].number_of_rotamers + 5), sizeof(ROTAMER)); rotamerlib[r].freq_array = (double *)calloc((rotamerlib[r].number_of_rotamers + 5), sizeof(double)); MAX_RESIMERS += rotamerlib[r].number_of_rotamers; if(rotamerlib[r].number_of_rotamers > MAX_ROTAMERS) MAX_ROTAMERS = rotamerlib[r].number_of_rotamers + 5; /* loop to read each set of coordinates */ for(i=1;i<=rotamerlib[r].number_of_rotamers+3;++i) /* make room for more 'rotamers' from other files */ { rotamerlib[r].rotamer[i].chi = (double *)calloc(rotamerlib[r].numChi+2,sizeof(double)); rotamerlib[r].rotamer[i].std_dev = (double *)calloc(rotamerlib[r].numChi+2,sizeof(double)); rotamerlib[r].rotamer[i].chi[1]=i; rotamerlib[r].rotamer[i].chi[2]=0; rotamerlib[r].rotamer[i].std_dev[1]=0; rotamerlib[r].rotamer[i].chi1_class=1; rotamerlib[r].rotamer[i].freq=1.0; rotamerlib[r].rotamer[i].ligand_coords = NULL; rotamerlib[r].rotamer[i].seqpos_text=NULL; if(i<=rotamerlib[r].number_of_rotamers) { fgets(line,MAXLINE,ligand_coords_file_ptr); rotamerlib[r].rotamer[i].ligand_coords = (pdbATOM *)calloc(MAX_RES_SIZE,sizeof(pdbATOM)); read_ligand_coords(ligand_coords_file_ptr, rotamerlib[r].rotamer[i].ligand_coords); } } if(coord_file_flag==0) /* set input file ptr */ input = ligand_coords_file_ptr; else fclose(ligand_coords_file_ptr); } else if(strstr(dummystring,"GENER")!=0 || strstr(dummystring,"gener")!=0) // generate ligamers internally { ligand_flag=1; numChi = 12; rotamerlib[r].generate_ligamers_flag = 1; // returns input file ptr to generated ligamer file generate_ligamers(input, &rotamerlib[r], line); } else { fprintf(stderr,"ERROR Must define number of chi or 'COORDINATES' or 'COORDINATE_FILE' for %s rotamerlib\n", rotamerlib[r].residuetype); exit(1); } } else { sscanf(line,"%s %d",rotamerlib[r].residuetype, &numChi); } } if(library_spread_ratio > 2) { fprintf(stderr,"ERROR Don't spread beyond chi2 for %s rotamerlib entry\n",rotamerlib[r].residuetype); exit(1); } library_spread_ratio = (int)pow(3,library_spread_ratio); if(num_library_rotamers==0) /* ie: dihedrals or transform matrix given, so count the number of rotamers */ { fp = ftell(input); fseek(input, fp, 0); fgets(line,MAXLINE,input); while(word_count(line)>3) { ++num_library_rotamers; fgets(line,MAXLINE,input); } fseek(input, fp, 0); /* reset file pointer */ rotamerlib[r].number_of_rotamers = num_library_rotamers*library_spread_ratio; } } if(read_coordinates_flag==0) { if(numChi>5 && numChi<12) { fprintf(stderr,"ERROR Can't deal with more than 5 dihedrals for residue %s\n",rotamerlib[r].residuetype); exit(1); } rotamerlib[r].numChi=numChi; rotamerlib[r].rotamer = (ROTAMER *)calloc((rotamerlib[r].number_of_rotamers + 3), sizeof(ROTAMER)); rotamerlib[r].freq_array = (double *)calloc((rotamerlib[r].number_of_rotamers + 3), sizeof(double)); MAX_RESIMERS += rotamerlib[r].number_of_rotamers; if(rotamerlib[r].number_of_rotamers > MAX_ROTAMERS) MAX_ROTAMERS = rotamerlib[r].number_of_rotamers + 3; for(q=1;q<=rotamerlib[r].number_of_rotamers;++q) { rotamerlib[r].freq_array[q] = 0; rotamerlib[r].rotamer[q].rotamerlet=NULL; rotamerlib[r].rotamer[q].born_factor=NULL; rotamerlib[r].rotamer[q].first_coulombic=NULL; rotamerlib[r].rotamer[q].coulombic=NULL; rotamerlib[r].rotamer[q].seqpos_text=NULL; } n=1; for(i=1;i<=num_library_rotamers;++i) { fgets(line,MAXLINE,input); dummyrotamer.chi = (double *)calloc(numChi+2,sizeof(double)); dummyrotamer.std_dev = (double *)calloc(numChi+2,sizeof(double)); dummyrotamer.ligand_coords = NULL; dummyrotamer.seqpos_text=NULL; for(q=0;q<=numChi+1;++q) { dummyrotamer.chi[q]=0; dummyrotamer.std_dev[q]=0; } if(numChi==1) { sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1]); // position-specific rotamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %s", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], dummyrotamer.seqpos_text); } } else if(numChi==2) { sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2]); // position-specific rotamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %s", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], dummyrotamer.seqpos_text); } } else if(numChi==3) { sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %lf%lf", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], &dummyrotamer.chi[3], &dummyrotamer.std_dev[3]); // position-specific rotamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %lf%lf %s", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], &dummyrotamer.chi[3], &dummyrotamer.std_dev[3], dummyrotamer.seqpos_text); } } else if(numChi==4) { sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %lf%lf %lf%lf", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], &dummyrotamer.chi[3], &dummyrotamer.std_dev[3], &dummyrotamer.chi[4], &dummyrotamer.std_dev[4]); // position-specific rotamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %lf%lf %lf%lf %s", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], &dummyrotamer.chi[3], &dummyrotamer.std_dev[3], &dummyrotamer.chi[4], &dummyrotamer.std_dev[4], dummyrotamer.seqpos_text); } } else if(numChi==5) { sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %lf%lf %lf%lf %lf%lf", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], &dummyrotamer.chi[3], &dummyrotamer.std_dev[3], &dummyrotamer.chi[4], &dummyrotamer.std_dev[4], &dummyrotamer.chi[5], &dummyrotamer.std_dev[5]); // position-specific rotamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(line,"%s %d%lf%lf%lf %lf%lf %lf %lf%lf%lf %lf%lf %lf%lf %lf%lf %lf%lf %lf%lf %s", doh, &dummyrotamer.chi1_class,&duh[2],&duh[3],&duh[4], &duh[5],&duh[6], &dummyrotamer.freq,&duh[7], &dummyrotamer.chi1_freq,&duh[9], &dummyrotamer.chi[1], &dummyrotamer.std_dev[1], &dummyrotamer.chi[2], &dummyrotamer.std_dev[2], &dummyrotamer.chi[3], &dummyrotamer.std_dev[3], &dummyrotamer.chi[4], &dummyrotamer.std_dev[4], &dummyrotamer.chi[5], &dummyrotamer.std_dev[5], dummyrotamer.seqpos_text); } } else if(numChi==12) /* transformation matrix for rigid ligand */ { sscanf(line,"%lf%lf%lf%lf %lf%lf%lf%lf %lf%lf%lf%lf", &dummyrotamer.chi[1], &dummyrotamer.chi[2], &dummyrotamer.chi[3], &dummyrotamer.chi[4], &dummyrotamer.chi[5], &dummyrotamer.chi[6], &dummyrotamer.chi[7], &dummyrotamer.chi[8], &dummyrotamer.chi[9], &dummyrotamer.chi[10], &dummyrotamer.chi[11], &dummyrotamer.chi[12]); check_transform_matrix(dummyrotamer.chi); dummyrotamer.chi1_class=1; dummyrotamer.freq=1.0; dummyrotamer.ligand_coords = NULL; dummyrotamer.seqpos_text=NULL; ligand_flag=1; // position-specific ligamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(line,"%lf%lf%lf%lf %lf%lf%lf%lf %lf%lf%lf%lf %s", &dummyrotamer.chi[1], &dummyrotamer.chi[2], &dummyrotamer.chi[3], &dummyrotamer.chi[4], &dummyrotamer.chi[5], &dummyrotamer.chi[6], &dummyrotamer.chi[7], &dummyrotamer.chi[8], &dummyrotamer.chi[9], &dummyrotamer.chi[10], &dummyrotamer.chi[11], &dummyrotamer.chi[12], dummyrotamer.seqpos_text); } } else /* combination of transforms and rotatable dihedral bonds */ { /* first 12 values are transform matrix elements, the remainder are dihedrals */ /* in principle, this should allow one to move any flexible ligand both internally and externally */ /* in addition to small molecules, one may even represent a protein or peptide in this manner */ ligand_flag=1; word = (char **)calloc(numChi +5,sizeof(char *)); for(p=1;p<=numChi+2;++p) word[p] = (char *)calloc(MAXLINE,sizeof(char)); extract_words(line, word); for(p=1;p<=numChi;++p) sscanf(word[p],"%lf",&dummyrotamer.chi[p]); check_transform_matrix(dummyrotamer.chi); dummyrotamer.chi1_class=1; dummyrotamer.freq=1.0; dummyrotamer.ligand_coords = NULL; dummyrotamer.seqpos_text=NULL; // position-specific ligamer if(word_count(line) == numChi+1) { dummyrotamer.seqpos_text = (char *)calloc(10,sizeof(char)); sscanf(word[numChi+1],"%s",dummyrotamer.seqpos_text); } for(p=1;p<=numChi+2;++p) free_memory(word[p]); free_memory(word); } dummyrotamer.freq = (dummyrotamer.freq/100); /* convert % -> fraction */ dummyrotamer.chi1_freq = (dummyrotamer.chi1_freq/100); /* spread rotamers if specified */ if(library_spread_ratio != 1) { for(m=-1;m<=1;++m) { if(library_spread_ratio == 9) for(d=-1;d<=1;++d) { rotamerlib[r].rotamer[n] = dummyrotamer; rotamerlib[r].rotamer[n].freq += d*EPS + m*EPS; /* just so the rotamers have unique freqs */ rotamerlib[r].rotamer[n].chi = (double *)calloc(numChi+2,sizeof(double)); rotamerlib[r].rotamer[n].std_dev = (double *)calloc(numChi+2,sizeof(double)); rotamerlib[r].rotamer[n].half_std_dev = (double *)calloc(numChi+2,sizeof(double)); for(p=0;p<=rotamerlib[r].numChi+1;++p) { rotamerlib[r].rotamer[n].chi[p]=dummyrotamer.chi[p]; rotamerlib[r].rotamer[n].std_dev[p]=dummyrotamer.std_dev[p]; rotamerlib[r].rotamer[n].half_std_dev[p]=0.49*dummyrotamer.std_dev[p]; } rotamerlib[r].rotamer[n].chi[1] = rotamerlib[r].rotamer[n].chi[1] + m*rotamerlib[r].rotamer[n].std_dev[1]; rotamerlib[r].rotamer[n].chi[2] = rotamerlib[r].rotamer[n].chi[2] + d*rotamerlib[r].rotamer[n].std_dev[2]; ++n; } if(library_spread_ratio == 3) { rotamerlib[r].rotamer[n] = dummyrotamer; rotamerlib[r].rotamer[n].freq += m*EPS; /* just so the rotamers have unique freqs */ rotamerlib[r].rotamer[n].chi = (double *)calloc(numChi+2,sizeof(double)); rotamerlib[r].rotamer[n].std_dev = (double *)calloc(numChi+2,sizeof(double)); rotamerlib[r].rotamer[n].half_std_dev = (double *)calloc(numChi+2,sizeof(double)); for(p=0;p<=rotamerlib[r].numChi+1;++p) { rotamerlib[r].rotamer[n].chi[p]=dummyrotamer.chi[p]; rotamerlib[r].rotamer[n].std_dev[p]=dummyrotamer.std_dev[p]; rotamerlib[r].rotamer[n].half_std_dev[p]=0.49*dummyrotamer.std_dev[p]; } rotamerlib[r].rotamer[n].chi[1] = rotamerlib[r].rotamer[n].chi[1] + m*rotamerlib[r].rotamer[n].std_dev[1]; ++n; } } } else { rotamerlib[r].rotamer[n] = dummyrotamer; rotamerlib[r].rotamer[n].chi = (double *)calloc(numChi+2,sizeof(double)); rotamerlib[r].rotamer[n].std_dev = (double *)calloc(numChi+2,sizeof(double)); rotamerlib[r].rotamer[n].half_std_dev = (double *)calloc(numChi+2,sizeof(double)); for(p=0;p<=rotamerlib[r].numChi+1;++p) { rotamerlib[r].rotamer[n].chi[p]=dummyrotamer.chi[p]; rotamerlib[r].rotamer[n].std_dev[p]=dummyrotamer.std_dev[p]; rotamerlib[r].rotamer[n].half_std_dev[p]=0.49*dummyrotamer.std_dev[p]; } ++n; } } } denom = 0; for(i=1;i<=rotamerlib[r].number_of_rotamers;++i) denom += rotamerlib[r].rotamer[i].freq; for(i=1;i<=rotamerlib[r].number_of_rotamers;++i) { rotamerlib[r].rotamer[i].freq = rotamerlib[r].rotamer[i].freq/denom; rotamerlib[r].rotamer[i].native_rotamer_flag = 0; } /* sort rotamers by database frequency */ if(ligand_flag==0) sort_ROTAMER(&first,&rotamerlib[r].number_of_rotamers,rotamerlib[r].rotamer); for(i=1;i<=rotamerlib[r].number_of_rotamers;++i) rotamerlib[r].rotamer[i].index = i; rotamerlib[r].freq_array[0] = 0; rotamerlib[r].freq_array[1] = 0; rotamerlib[r].freq_array[2] = 0; /* for local minimization....generate rotamerlets */ rotamerlib[r].num_rotamerlets = 0; if(rotamerlib[r].number_of_rotamers>=3) { if(numChi >= 2) { rotamerlib[r].num_rotamerlets = 9; if(LOTS_OF_ROTAMERLETS_FLAG == 1) if(numChi >= 3) rotamerlib[r].num_rotamerlets = 27; } else rotamerlib[r].num_rotamerlets = 3; } if(rotamerlib[r].number_of_rotamers==1) rotamerlib[r].num_rotamerlets = 1; for(i=1;i<=rotamerlib[r].number_of_rotamers;++i) { rotamerlib[r].rotamer[i].rank = i; rotamerlib[r].freq_array[i] = rotamerlib[r].rotamer[i].freq + rotamerlib[r].freq_array[i-1]; if(rotamerlib[r].num_rotamerlets!=0) { rotamerlib[r].rotamer[i].rotamerlet = (double **)calloc(rotamerlib[r].num_rotamerlets+2, sizeof(double *)); for(q=1;q<=rotamerlib[r].num_rotamerlets;++q) { rotamerlib[r].rotamer[i].rotamerlet[q] = (double *)calloc(rotamerlib[r].numChi+2, sizeof(double)); for(n=0;n<=rotamerlib[r].numChi+1;++n) rotamerlib[r].rotamer[i].rotamerlet[q][n] = rotamerlib[r].rotamer[i].chi[n]; } if(rotamerlib[r].num_rotamerlets >= 3) { q=1; for(m=-1;m<=1;++m) { rotamerlib[r].rotamer[i].rotamerlet[q][1] = rotamerlib[r].rotamer[i].chi[1] + ROTAMER_SPREAD_FACTOR*m*rotamerlib[r].rotamer[i].std_dev[1]; if(rotamerlib[r].num_rotamerlets>=9) { n=q; for(d=-1;d<=1;++d) { rotamerlib[r].rotamer[i].rotamerlet[q][1] = rotamerlib[r].rotamer[i].rotamerlet[n][1]; rotamerlib[r].rotamer[i].rotamerlet[q][2] = rotamerlib[r].rotamer[i].rotamerlet[q][2] + ROTAMER_SPREAD_FACTOR*d*rotamerlib[r].rotamer[i].std_dev[2]; if(rotamerlib[r].num_rotamerlets>=27) { nn = q; for(dd = -1;dd<=1;++dd) { rotamerlib[r].rotamer[i].rotamerlet[q][1] = rotamerlib[r].rotamer[i].rotamerlet[nn][1]; rotamerlib[r].rotamer[i].rotamerlet[q][2] = rotamerlib[r].rotamer[i].rotamerlet[nn][2]; rotamerlib[r].rotamer[i].rotamerlet[q][3] = rotamerlib[r].rotamer[i].rotamerlet[q][3] + ROTAMER_SPREAD_FACTOR*dd*rotamerlib[r].rotamer[i].std_dev[3]; ++q; } } else ++q; } } else ++q; } } } } rotamerlib[r].freq_array[rotamerlib[r].number_of_rotamers+1] = ENDFLAG; fgets(line,MAXLINE,input); ++r; } while(rtorsion_class[1] = atomparam[i].torsion_class; for(k=2;k<=4;++k) { n=1; while(strcmp(atomparam[n].torsiontype, dummytorsion[j].atom[k])!=0) ++n; atomparam[i].torsion->torsion_class[k] = atomparam[n].torsion_class; } atomparam[i].torsion->next = (TORSION *)malloc(sizeof(TORSION)); atomparam[i].torsion = atomparam[i].torsion->next; atomparam[i].torsion->next = NULL; } } ++i; } free_memory(dummytorsion); free_memory(line); } /* reads in residue topologies from resdatafilename into an allocated array of RESPARAM. if count_max_res_size_flag=1, merely counts the size of residues, adjusts MAX_RES_SIZE, and returns. called by read_forcefield.cpp: read_forcefield and ligand_stuff.cpp: read_ligand_data.cpp */ void read_residuedata(char *resdatafilename, RESPARAM *resparam, ATOMPARAM *atomparam, int count_max_res_size_flag) { FILE *input; char *line, *dummystring; int j, flag, k,i, q, n, m, methyl_flag, methyl_ctr, MeH_ctr; int max_methyl_actual, version_flag; static int warning_flag=0; ATOMRESPARAM temp_atom; extern RESPARAM *N_TERM_RESPARAM, *C_TERM_RESPARAM; extern char *ROTAMERFILE; double P_prot, P_unprot; double m1, m2, m3; double E_rand, pi, pj, ala_tripep, refstate_sasa; extern double TEMPERATURE, PH, RT; extern int SPECIFICITY_FLAG; line = (char *)calloc(MAXLINE, sizeof(char)); if(count_max_res_size_flag==1) { input = fopen_file(resdatafilename,"r"); fgets(line,MAXLINE,input); while(strncmp(line,"START",5)!=0) fgets(line,MAXLINE,input); while(strncmp(line,"STOP",4)!=0) { if(strncmp(line,"START",5)==0 || strncmp(line,"done",4)==0) { fgets(line,MAXLINE,input); n=1; while(strncmp(line,"done",4)!=0 && strncmp(line,"STOP",4)!=0) { ++n; fgets(line,MAXLINE,input); } if(n+10>MAX_RES_SIZE) MAX_RES_SIZE=n+10; } } fclose(input); free_memory(line); return; } dummystring = (char *)calloc(MAXLINE, sizeof(char)); max_methyl_actual = 0; strcpy(resparam[0].residuetype,""); flag=0;q=0; input = fopen_file(resdatafilename,"r"); version_flag=1; fgets(line,MAXLINE,input); while(strncmp(line,"START",5)!=0) { /* parse ROTAMERFILE if defined in this file and not defined in user input file */ if(strstr(line,"ROTAMER")!=0) { if(ROTAMERFILE == NULL) /* if ROTAMERFILE has been defined in the input file, ignore this line */ { ROTAMERFILE = (char *)calloc(MAXLINE,sizeof(char)); sscanf(line,"%s %s",dummystring,ROTAMERFILE); } } /* determine version number */ if(strstr(line,"VERSION 2")!=0) version_flag=2; fgets(line,MAXLINE,input); } if(strstr(line,"VERSION 2")!=0) version_flag=2; free_memory(dummystring); while(strncmp(line,"STOP",4)!=0) { if(strncmp(line,"START",5)==0 || strncmp(line,"done",4)==0) { ++q; flag=1; j=4; /* this way, param.atom.[j].atomnumber = index */ resparam[q].atom = (ATOMRESPARAM *)calloc(MAX_RES_SIZE, sizeof(ATOMRESPARAM)); resparam[q].sorted_atom = (ATOMRESPARAM *)calloc(MAX_RES_SIZE, sizeof(ATOMRESPARAM)); resparam[q].atom_ptr=NULL; resparam[q].non_MeH_atom_ptr=NULL; resparam[q].MeH_atom_ptr=NULL; resparam[q].numberofAtoms = 0; resparam[q].num_non_H = 0; resparam[q].num_charged = 0; resparam[q].overall_charge = 0; resparam[q].charged_CB_flag = 0; resparam[q].num_charge_H = 0; resparam[q].num_rotatable_bonds = 0; resparam[q].num_methyl_groups=0; methyl_flag=0; resparam[q].num_methyl_confs=0; resparam[q].numberofSideAtoms = 0; resparam[q].E_unfolded=0; resparam[q].ligand_flag=0; resparam[q].num_hbonding_sidechain_atoms = 0; resparam[q].numhbondsideatoms_minus1=0; resparam[q].max_solvation = 1000; resparam[q].max_vdw = 10.0; resparam[q].sidechain_hbond_atom_ptr = NULL; fgets(line,MAXLINE,input); resparam[q].protonation_flag = 0; refstate_sasa = 0; if(version_flag==1) { sscanf(line,"%s %s %d %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf",resparam[q].residuetype, resparam[q].one_letter_code, &resparam[q].residuecode, &resparam[q].pKa, &resparam[q].max_vdw, &m1, &m2, &m3, &resparam[q].S, &resparam[q].E_random, &resparam[q].E_specificity, &resparam[q].E_specificity_core, &resparam[q].E_specificity_interfacial, &resparam[q].E_specificity_surface, &refstate_sasa ); if(m2!=0) resparam[q].E_tripeptide = m1 + m2*TEMPERATURE + m3*TEMPERATURE*TEMPERATURE + HYDROPHOBIC_ASP*refstate_sasa; else { if(m1!=0) resparam[q].E_tripeptide = m1 - TEMPERATURE*resparam[q].S; else resparam[q].E_tripeptide = 0.0; } } else { if(word_count(line)==10) /* global specificity energies */ { sscanf(line,"%s %s %d %lf %lf %lf %lf %lf %lf %lf", resparam[q].residuetype, resparam[q].one_letter_code, &resparam[q].residuecode, &resparam[q].pKa, &resparam[q].max_vdw, &resparam[q].E_tripeptide, &resparam[q].S, &resparam[q].E_random, &resparam[q].E_specificity, &resparam[q].E_unfolded ); if(resparam[q].E_tripeptide==0) resparam[q].S=0.0; resparam[q].E_specificity_surface = resparam[q].E_specificity; resparam[q].E_specificity_core = resparam[q].E_specificity; resparam[q].E_specificity_interfacial = resparam[q].E_specificity; } else if(word_count(line)==3) /* probably for a ligand */ { sscanf(line,"%s %s %d", resparam[q].residuetype, resparam[q].one_letter_code, &resparam[q].residuecode); resparam[q].pKa=20.0; resparam[q].max_vdw=10.0; resparam[q].E_tripeptide=0; resparam[q].S=0; resparam[q].E_random=0; resparam[q].E_unfolded=0; resparam[q].E_specificity=0; resparam[q].E_specificity_surface = resparam[q].E_specificity; resparam[q].E_specificity_core = resparam[q].E_specificity; resparam[q].E_specificity_interfacial = resparam[q].E_specificity; } else if(word_count(line)==4) /* probably for a ligand, but with max_vdw specified */ { sscanf(line,"%s %s %d %lf", resparam[q].residuetype, resparam[q].one_letter_code, &resparam[q].residuecode, &resparam[q].max_vdw); resparam[q].pKa=20.0; resparam[q].E_tripeptide=0; resparam[q].S=0; resparam[q].E_random=0; resparam[q].E_specificity=0; resparam[q].E_unfolded=0; resparam[q].E_specificity_surface = resparam[q].E_specificity; resparam[q].E_specificity_core = resparam[q].E_specificity; resparam[q].E_specificity_interfacial = resparam[q].E_specificity; } else if(word_count(line)==5) /* probably for a ligand, but with max_vdw and max_solvation specified */ { sscanf(line,"%s %s %d %lf %lf", resparam[q].residuetype, resparam[q].one_letter_code, &resparam[q].residuecode, &resparam[q].max_vdw, &resparam[q].max_solvation); resparam[q].pKa=20.0; resparam[q].E_tripeptide=0; resparam[q].S=0; resparam[q].E_random=0; resparam[q].E_specificity=0; resparam[q].E_unfolded=0; resparam[q].E_specificity_surface = resparam[q].E_specificity; resparam[q].E_specificity_core = resparam[q].E_specificity; resparam[q].E_specificity_interfacial = resparam[q].E_specificity; } else if(word_count(line) >= 11) /* location-dependent specificity energies */ { sscanf(line,"%s %s %d %lf %lf %lf %lf %lf %lf %lf %lf", resparam[q].residuetype, resparam[q].one_letter_code, &resparam[q].residuecode, &resparam[q].pKa, &resparam[q].max_vdw, &resparam[q].E_tripeptide, &resparam[q].S, &resparam[q].E_random, &resparam[q].E_specificity_core, &resparam[q].E_specificity_interfacial, &resparam[q].E_specificity_surface ); resparam[q].E_unfolded = 0; resparam[q].E_specificity=0; if(resparam[q].E_tripeptide==0) resparam[q].S=0.0; } if(resparam[q].E_tripeptide==0) resparam[q].S=0.0; if(resparam[q].E_random==0) resparam[q].E_random = resparam[q].E_tripeptide; if(resparam[q].E_unfolded==0) resparam[q].E_unfolded = resparam[q].E_random; if(resparam[q].E_specificity_core == 0) resparam[q].E_specificity_core = resparam[q].E_specificity; if(resparam[q].E_specificity_interfacial == 0) resparam[q].E_specificity_interfacial = resparam[q].E_specificity; if(resparam[q].E_specificity_surface == 0) resparam[q].E_specificity_surface = resparam[q].E_specificity; resparam[q].E_unfolded = resparam[q].E_unfolded - TEMPERATURE*resparam[q].S; resparam[q].E_tripeptide = resparam[q].E_tripeptide - TEMPERATURE*resparam[q].S; resparam[q].E_reference = resparam[q].E_random - TEMPERATURE*resparam[q].S; } if(SPECIFICITY_FLAG==0) { resparam[q].E_specificity_core=0; resparam[q].E_specificity_surface=0; resparam[q].E_specificity_interfacial=0; resparam[q].E_specificity=0; } fgets(line,MAXLINE,input); for(i=0;i<4;++i) { strcpy(resparam[q].atom[i].atomname,""); strcpy(resparam[q].atom[i].atomtype,""); } for(i=4;i0) { if(temp_atom.donorflag[0] != ' ') ++resparam[q].num_hbonding_sidechain_atoms; } else { if(temp_atom.donorflag[0] != ' ') { if(strcmp(temp_atom.atomname, "N")==0) temp_atom.E_hbond_specificity = DONOR_UNSAT_HBOND_REFSTATE; if(strcmp(temp_atom.atomname, "H")==0) temp_atom.E_hbond_specificity = H_UNSAT_HBOND_REFSTATE; if(strcmp(temp_atom.atomname, "O")==0) temp_atom.E_hbond_specificity = ACCEPTOR_UNSAT_HBOND_REFSTATE; } } /* find and attach pointer to the appropriate ATOMPARAM entry */ k=0; while(k0 || (strcmp(resparam[q].residuetype, "PRO")==0 && resparam[q].atom[j].atomnumber>8 && resparam[q].atom[j].atomnumber<17)) { ++resparam[q].numberofSideAtoms; if(strncmp(resparam[q].atom[j].atomtype, "H", 1)!=0) /* count-non-backbone-non-hydrogen atoms */ ++resparam[q].num_non_H; if(resparam[q].atom[j].charge!=0) /* non-backbone charged atoms */ { ++resparam[q].num_charged; resparam[q].overall_charge = resparam[q].overall_charge + resparam[q].atom[j].charge; if(strncmp(resparam[q].atom[j].atomtype, "H", 1)==0) /* count charged H */ ++resparam[q].num_charge_H; } } if(resparam[q].atom[j].other_info>0 || strcmp(resparam[q].atom[j].atomname, "CB")==0) { if(strcmp(resparam[q].atom[j].atomtype, "CH3")==0) /* methyl C */ { ++resparam[q].num_methyl_groups; methyl_flag = resparam[q].atom[j].atomnumber; } if(resparam[q].num_methyl_groups!=0) if(resparam[q].atom[j].contactAtomNumber == methyl_flag) resparam[q].atom[j].methyl_flag = methyl_flag; } if(resparam[q].atom[j].other_info<0) if(strcmp(resparam[q].atom[j].atomname, "CB")==0) /* charged CB flag */ if(resparam[q].atom[j].charge!=0) { resparam[q].overall_charge = resparam[q].overall_charge + resparam[q].atom[j].charge; resparam[q].charged_CB_flag = 1; } ++j; } fgets(line,MAXLINE,input); } resparam[q].atom[j].atomnumber = ENDFLAG; if(strcmp(resparam[q].residuetype, "NTE")==0) { resparam[q].overall_charge = 1; resparam[q].num_charged = 7; resparam[q].num_charge_H = 3; } else if(strcmp(resparam[q].residuetype, "NTN")==0) { resparam[q].overall_charge = 0; resparam[q].num_charged = 7; resparam[q].num_charge_H = 3; } else if(strcmp(resparam[q].residuetype, "CTE")==0) { resparam[q].overall_charge = -1; resparam[q].num_charged = 7; resparam[q].num_charge_H = 2; } else if(strcmp(resparam[q].residuetype, "CTN")==0) { resparam[q].overall_charge = 0; resparam[q].num_charged = 7; resparam[q].num_charge_H = 2; } if(resparam[q].num_charged + resparam[q].charged_CB_flag > 0) resparam[q].charged_flag = 1; else resparam[q].charged_flag = 0; resparam[q].overall_charge_pH_avg = resparam[q].overall_charge; resparam[q].atom_ptr = NULL; resparam[q].non_MeH_atom_ptr = NULL; resparam[q].MeH_atom_ptr = NULL; if(resparam[q].num_hbonding_sidechain_atoms!=0) { resparam[q].sidechain_hbond_atom_ptr = (ATOMRESPARAM **)calloc((resparam[q].num_hbonding_sidechain_atoms+3), sizeof(ATOMRESPARAM *)); j=8; n=1; resparam[q].sidechain_hbond_atom_ptr[0] = NULL; while(resparam[q].atom[j].atomnumber!=ENDFLAG) { if(resparam[q].atom[j].other_info > 0) if(resparam[q].atom[j].donorflag[0] != ' ') { resparam[q].sidechain_hbond_atom_ptr[n] = &resparam[q].atom[j]; ++n; } ++j; } resparam[q].sidechain_hbond_atom_ptr[n] = NULL; } if(resparam[q].numberofSideAtoms!=0) { resparam[q].atom_ptr = (ATOMRESPARAM **)calloc((resparam[q].numberofSideAtoms+3), sizeof(ATOMRESPARAM *)); if(resparam[q].num_methyl_groups != 0) { resparam[q].non_MeH_atom_ptr = (ATOMRESPARAM **)calloc((resparam[q].numberofSideAtoms+3 - 3*resparam[q].num_methyl_groups), sizeof(ATOMRESPARAM *)); resparam[q].MeH_atom_ptr = (ATOMRESPARAM ***)calloc(resparam[q].num_methyl_groups+1, sizeof(ATOMRESPARAM **)); } else { resparam[q].non_MeH_atom_ptr = NULL; resparam[q].MeH_atom_ptr = NULL; } n=1; j=8; m=0; methyl_ctr=0; MeH_ctr=0; resparam[q].atom_ptr[0] = NULL; while(n<=(resparam[q].numberofSideAtoms +1)) { /* array of pointers to atomresparam entries for sidechain atoms (including CB) used for calculating var_var energies in the lookup table */ resparam[q].atom_ptr[n] = &resparam[q].atom[j]; /* methyl groups for local minimization */ if(resparam[q].num_methyl_groups != 0) { if(resparam[q].atom[j].methyl_flag!=0) { if(MeH_ctr%3 == 0) { ++methyl_ctr; resparam[q].MeH_atom_ptr[methyl_ctr] = (ATOMRESPARAM **)calloc(5,sizeof(ATOMRESPARAM *)); resparam[q].MeH_atom_ptr[methyl_ctr][4] = NULL; MeH_ctr = 0; } ++MeH_ctr; resparam[q].MeH_atom_ptr[methyl_ctr][MeH_ctr] = &resparam[q].atom[j]; } else { ++m; resparam[q].non_MeH_atom_ptr[m] = &resparam[q].atom[j]; } } ++n; ++j; } resparam[q].atom_ptr[n] = NULL; if(resparam[q].num_methyl_groups != 0) { resparam[q].non_MeH_atom_ptr[m+1] = NULL; resparam[q].num_methyl_confs = (int)pow(3, resparam[q].num_methyl_groups); if(resparam[q].num_methyl_groups > max_methyl_actual) max_methyl_actual = resparam[q].num_methyl_groups; } } } if(resparam[q].num_hbonding_sidechain_atoms>0) resparam[q].numhbondsideatoms_minus1 = resparam[q].num_hbonding_sidechain_atoms - 1; if(strcmp(resparam[q].residuetype, "NTE")==0) N_TERM_RESPARAM = &resparam[q]; else if(strcmp(resparam[q].residuetype, "CTE")==0) C_TERM_RESPARAM = &resparam[q]; else if(strcmp(resparam[q].residuetype, "PRO")==0) PRO_RESPARAM= &resparam[q]; else if(strcmp(resparam[q].residuetype, "GLY")==0) GLY_RESPARAM= &resparam[q]; } /* pseudo-sidechain for born/sasa calcs */ if(PSEUDO_SIDECHAIN_ATOMRESPARAM.atomnumber==ENDFLAG) { fgets(line,MAXLINE,input); sscanf(line,"%d %s %s %d %d %d %d %d %d %d %lf %lf %lf %lf", &PSEUDO_SIDECHAIN_ATOMRESPARAM.atomnumber, PSEUDO_SIDECHAIN_ATOMRESPARAM.atomname, PSEUDO_SIDECHAIN_ATOMRESPARAM.atomtype, &PSEUDO_SIDECHAIN_ATOMRESPARAM.other_info, &PSEUDO_SIDECHAIN_ATOMRESPARAM.contactAtomNumber, &PSEUDO_SIDECHAIN_ATOMRESPARAM.twobondAtomNumber, &PSEUDO_SIDECHAIN_ATOMRESPARAM.dihedAtomNumber, &PSEUDO_SIDECHAIN_ATOMRESPARAM.ringContactAtom, &PSEUDO_SIDECHAIN_ATOMRESPARAM.ringTwoBondAtom, &PSEUDO_SIDECHAIN_ATOMRESPARAM.ringDihedAtom, &PSEUDO_SIDECHAIN_ATOMRESPARAM.bondlength, &PSEUDO_SIDECHAIN_ATOMRESPARAM.bondangle, &PSEUDO_SIDECHAIN_ATOMRESPARAM.dihedral, &PSEUDO_SIDECHAIN_ATOMRESPARAM.charge); PSEUDO_SIDECHAIN_ATOMRESPARAM.donorflag[0] = ' '; } fclose(input); q=q+1; while(q resparam[j].overall_charge) resparam[i].protonation_flag=1; else resparam[j].protonation_flag=1; if(PH > 0) { P_prot = exp(-(PH - resparam[i].pKa)*LN10); P_prot = P_prot/(1.0 + P_prot); /* probability of the protonated form */ P_unprot = 1.0 - P_prot; /* probability of the unprotonated form */ if(resparam[i].protonation_flag==1) { pi = P_prot; pj = P_unprot; } else { pj = P_prot; pi = P_unprot; } /* boltz avg energy of protonated and deprot random and specificity states */ E_rand = pi*resparam[i].E_random + pj*resparam[j].E_random; resparam[i].E_random = E_rand; resparam[j].E_random = E_rand; /* penalize the charged form for core positions, since charging it is unfavorable */ if(fabs(resparam[i].overall_charge) > EPS) resparam[i].E_specificity_core -= RT*fabs(resparam[i].pKa - PH)*LN10; else resparam[j].E_specificity_core -= RT*fabs(resparam[j].pKa - PH)*LN10; E_rand = pi*resparam[i].E_specificity_core + pj*resparam[j].E_specificity_core; resparam[i].E_specificity_core = E_rand; resparam[j].E_specificity_core = E_rand; E_rand = pi*resparam[i].E_specificity_interfacial + pj*resparam[j].E_specificity_interfacial; resparam[i].E_specificity_interfacial = E_rand; resparam[j].E_specificity_interfacial = E_rand; E_rand = pi*resparam[i].E_specificity_surface + pj*resparam[j].E_specificity_surface; resparam[i].E_specificity_surface = E_rand; resparam[j].E_specificity_surface = E_rand; E_rand = pi*resparam[i].overall_charge + pj*resparam[j].overall_charge; resparam[i].overall_charge_pH_avg = E_rand; resparam[j].overall_charge_pH_avg = E_rand; resparam[i].E_reference = pi*resparam[i].E_reference + pj*resparam[j].E_reference; resparam[j].E_reference = resparam[i].E_reference; if(version_flag==1) /* E_unf and E_ref are same for old version */ { resparam[i].E_unfolded = resparam[i].E_reference; resparam[j].E_unfolded = resparam[j].E_reference; } else /* new version, E_unfolded is different */ { E_rand = pi*resparam[i].E_unfolded + pj*resparam[j].E_unfolded; resparam[i].E_unfolded = E_rand; resparam[j].E_unfolded = E_rand; } } } } ++i; } free_memory(line); } /* reads forcefield information into atomparam, resparam, and rotamerlib arrays (allocated by calling function). sets up a lot of cross-pointers between these data structures,initializes some matricies, etc. read SASA tables (requires MAX_ATOMS to be defined) and ligand parameter files as well. Called by input_stuff.cpp: input_stuff. */ void read_forcefield(char *atomdatafile, ATOMPARAM *atomparam, RESPARAM *resparam, ROTAMERLIB *rotamerlib) { FILE *input; char *line, *torsionfilename, *resparamfilename, *dummystring, *SASA_points_file, *energy_dir; int i, j, k; int atom_classes, torsion_classes; /* Atom class counter */ extern char *RESPARAMFILE, *ROTAMERFILE, **LIGANDFILENAME; double beta; extern double VDW_CUTOFF, OVERALL_ENERGY_SCALE, VDW_ATTRACTIVE_FACTOR, VDW_CLASH_ENERGY; extern double WEIGHT_VDW, WEIGHT_ELECTROSTATICS, WEIGHT_1_4, GENERAL_ASP, HYDROPHOBIC_ASP, CORESURF_MEA_BORN,CORESURF_MEA_SASA,VDW_REPULSIVE_FACTOR, DEBYE_SCALE_FACTOR, TEMPERATURE, RT; extern double VDW_LINEAR_START_POINT, UNSAT_HBOND_REFSTATE; extern int INTRA_ROTAMER_FLAG, INTRA_ROTAMER_FLAG_14, LOCAL_MINIMIZATION_FLAG; static int warning_flag=0; extern BACKBONE DUMMY_BKBN, PROLINE_BKBN; extern int HBOND_FUNCTION_TYPE, HBOND_FUNCTION_FLAG, SPECIFICITY_FLAG; MAX_RES_SIZE=0; PSEUDO_SIDECHAIN_ATOMRESPARAM.atomnumber=ENDFLAG; if(atomdatafile==NULL) { fprintf(stderr,"ERROR Must define FORCEFIELD in the input file\n"); exit(1); } // extract directory of energy function energy_dir = (char *)calloc(MAXLINE,sizeof(char)); i=0; while(atomdatafile[i]!='\0') // advance to end-of-line ++i; --i; // last character while(atomdatafile[i]!='/') // go back to the last '/' --i; ++i; j=0; while(j!=i) { energy_dir[j] = atomdatafile[j]; ++j; } energy_dir[j]='\0'; /* BACKBONE for a single residue */ DUMMY_BKBN.N.x = 27.200; DUMMY_BKBN.N.y = 23.980; DUMMY_BKBN.N.z = 2.804; DUMMY_BKBN.H.x = 27.966; DUMMY_BKBN.H.y = 24.066; DUMMY_BKBN.H.z = 2.106; DUMMY_BKBN.HA.x = 26.062; DUMMY_BKBN.HA.y = 25.375; DUMMY_BKBN.HA.z = 1.687; DUMMY_BKBN.CA.x = 26.309; DUMMY_BKBN.CA.y = 25.182; DUMMY_BKBN.CA.z = 2.721; DUMMY_BKBN.C.x = 27.035; DUMMY_BKBN.C.y = 26.409; DUMMY_BKBN.C.z = 3.303; DUMMY_BKBN.CB.x = 25.017; DUMMY_BKBN.CB.y = 24.918; DUMMY_BKBN.CB.z = 3.515; DUMMY_BKBN.O.x = 28.189; DUMMY_BKBN.O.y = 26.329; DUMMY_BKBN.O.z = 3.672; /* coordinates from proline from insight2, amber minimized */ PROLINE_BKBN.N.x = -0.556; PROLINE_BKBN.N.y = -0.031; PROLINE_BKBN.N.z = 7.183; PROLINE_BKBN.H.x = -0.948; PROLINE_BKBN.H.y = -0.962; PROLINE_BKBN.H.z = 7.201; PROLINE_BKBN.HA.x = 0.949; PROLINE_BKBN.HA.y = -0.197; PROLINE_BKBN.HA.z = 5.740; PROLINE_BKBN.CA.x = 0.789; PROLINE_BKBN.CA.y = 0.245; PROLINE_BKBN.CA.z = 6.725; PROLINE_BKBN.C.x = 1.796; PROLINE_BKBN.C.y = -0.344; PROLINE_BKBN.C.z = 7.716; PROLINE_BKBN.CB.x = 0.820; PROLINE_BKBN.CB.y = 1.780; PROLINE_BKBN.CB.z = 6.637; PROLINE_BKBN.O.x = 2.145; PROLINE_BKBN.O.y = 0.281; PROLINE_BKBN.O.z = 8.716; line = (char *)calloc(MAXLINE, sizeof(char)); dummystring = (char *)calloc(MAXLINE, sizeof(char)); torsionfilename = NULL; resparamfilename = NULL; SASA_points_file = NULL; MAX_RESIMERS=3; MAX_ROTAMERS=0; atom_classes=0; torsion_classes=1; input = fopen_file(atomdatafile,"r"); fgets(line,MAXLINE,input); while(strncmp(line,"START",5)!=0) fgets(line,MAXLINE,input); fgets(line,MAXLINE,input); /* read data for atomtypes */ while(strncmp(line,"STOP",4)!=0) { if(word_count(line)!=0) // if 0, is a comment line { ++atom_classes; sscanf(line,"%s %s %lf %lf %lf %lf", atomparam[atom_classes].atomtype, atomparam[atom_classes].torsiontype, &atomparam[atom_classes].sigma, &atomparam[atom_classes].welldepth, &atomparam[atom_classes].united_radius, &atomparam[atom_classes].asp); /* for surface area calcs */ if(atomparam[atom_classes].united_radius > 0) { if(strcmp(atomparam[atom_classes].atomtype, "BHY")==0) /* for identifying buried hbond'ing atoms */ { if(warning_flag==0) fprintf(stderr,"BHY atoms in %s no longer supported .....ignoring\n",atomdatafile); warning_flag=1; } else atomparam[atom_classes].sasa_radius = atomparam[atom_classes].united_radius + WATER_RADIUS; } else atomparam[atom_classes].sasa_radius = 0; atomparam[atom_classes].sasa_diameter = 2*atomparam[atom_classes].sasa_radius; atomparam[atom_classes].sasa_radius2 = atomparam[atom_classes].sasa_radius*atomparam[atom_classes].sasa_radius; atomparam[atom_classes].sasa_sphere = PI_4*atomparam[atom_classes].sasa_radius2; atomparam[atom_classes].sa_sphere = PI_4*pow(atomparam[atom_classes].united_radius, 2); atomparam[atom_classes].volume = (4.0/3.0)*PI*pow(atomparam[atom_classes].united_radius, 3); atomparam[atom_classes].atom_class = atom_classes; atomparam[atom_classes].torsion=NULL; atomparam[atom_classes].first_torsion=NULL; if(atom_classes>1) { if(strcmp(atomparam[atom_classes].torsiontype,atomparam[atom_classes-1].torsiontype)!=0) ++torsion_classes; } /* set up some pointers for pseudoatoms */ if(strcmp(atomparam[atom_classes].torsiontype,"U")==0) { if(strcmp(atomparam[atom_classes].atomtype, "NBK")==0) { BKBN_N_ATOMRESPARAM.atom_ptr = &atomparam[atom_classes]; BKBN_N_ATOMRESPARAM.donorflag[0]=' '; BKBN_N_ATOMRESPARAM.E_hbond_specificity=0; } else if(strcmp(atomparam[atom_classes].atomtype, "CBK")==0) { BKBN_C_ATOMRESPARAM.atom_ptr = &atomparam[atom_classes]; BKBN_C_ATOMRESPARAM.donorflag[0]=' '; BKBN_C_ATOMRESPARAM.E_hbond_specificity=0; } else if(strcmp(atomparam[atom_classes].atomtype, "OBK")==0) { BKBN_O_ATOMRESPARAM.atom_ptr = &atomparam[atom_classes]; BKBN_O_ATOMRESPARAM.donorflag[0]=' '; BKBN_O_ATOMRESPARAM.E_hbond_specificity=0; } else if(strcmp(atomparam[atom_classes].atomtype, "CABK")==0) { BKBN_CA_ATOMRESPARAM.atom_ptr = &atomparam[atom_classes]; BKBN_CA_ATOMRESPARAM.donorflag[0]=' '; BKBN_CA_ATOMRESPARAM.E_hbond_specificity=0; } else if(strcmp(atomparam[atom_classes].atomtype, "HABK")==0) { BKBN_HA_ATOMRESPARAM.atom_ptr = &atomparam[atom_classes]; BKBN_HA_ATOMRESPARAM.donorflag[0]=' '; BKBN_HA_ATOMRESPARAM.E_hbond_specificity=0; } else if(strcmp(atomparam[atom_classes].atomtype, "PSE")==0) { PSEUDO_SIDECHAIN_ATOMRESPARAM.atom_ptr = &atomparam[atom_classes]; PSEUDO_SIDECHAIN_ATOMRESPARAM.donorflag[0]=' '; PSEUDO_SIDECHAIN_ATOMRESPARAM.E_hbond_specificity=0; } else if(strcmp(atomparam[atom_classes].atomtype, "PSEB")==0) { PSEB = PSEUDO_SIDECHAIN_ATOMRESPARAM; PSEB.atom_ptr = &atomparam[atom_classes]; PSEB.donorflag[0]=' '; PSEB.E_hbond_specificity=0; } } atomparam[atom_classes].torsion_class = torsion_classes; } fgets(line,MAXLINE,input); } /* define BHY buried hbond group probe atom */ ++atom_classes; BHY_CLASS = atom_classes; strcpy(atomparam[atom_classes].atomtype,"BHY"); strcpy(atomparam[atom_classes].torsiontype,"U"); atomparam[atom_classes].sigma = 0.0; atomparam[atom_classes].welldepth = 0.0; atomparam[atom_classes].united_radius = WATER_RADIUS; atomparam[atom_classes].asp = 0.0; atomparam[atom_classes].sasa_radius = atomparam[atom_classes].united_radius; BHY_ATOMRESPARAM = (ATOMRESPARAM *)malloc(sizeof(ATOMRESPARAM)); BHY_ATOMRESPARAM->atom_ptr = &atomparam[atom_classes]; sprintf(BHY_ATOMRESPARAM->atomname,"BHY"); sprintf(BHY_ATOMRESPARAM->atomtype,"BHY"); atomparam[atom_classes].sasa_diameter = 2*atomparam[atom_classes].sasa_radius; atomparam[atom_classes].sasa_radius2 = atomparam[atom_classes].sasa_radius*atomparam[atom_classes].sasa_radius; atomparam[atom_classes].sasa_sphere = PI_4*atomparam[atom_classes].sasa_radius2; atomparam[atom_classes].sa_sphere = PI_4*pow(atomparam[atom_classes].united_radius, 2); atomparam[atom_classes].volume = (4.0/3.0)*PI*pow(atomparam[atom_classes].united_radius, 3); atomparam[atom_classes].atom_class = atom_classes; atomparam[atom_classes].torsion_class = torsion_classes; atomparam[atom_classes].torsion=NULL; atomparam[atom_classes].first_torsion=NULL; /* read in atom-independent forcefield parameters */ while(fgets(line,MAXLINE,input) != NULL) { sscanf(line, "%s", dummystring); convert_string_to_all_caps(dummystring); if(strcmp(dummystring, "VDW_RADII_SCALE")==0) { if(VDW_RADII_SCALE == 0) /* if !=0, then it was defined in the user input file */ sscanf(line, "%s %lf", dummystring, &VDW_RADII_SCALE); } else if(strcmp(dummystring, "EP")==0) { if(INTERNAL_DIELECTRIC == 0) /* if !=0, then it was defined in the user input file */ sscanf(line, "%s %lf", dummystring, &INTERNAL_DIELECTRIC); } else if(strcmp(dummystring, "EHB")==0) { if(HBOND_DIELECTRIC == 0) /* if !=0, then it was defined in the user input file */ sscanf(line, "%s %lf", dummystring, &HBOND_DIELECTRIC); } /* these are deprecated, but still may be in some older files */ else if(strcmp(dummystring, "EXPOSED_HBOND_REWARD")==0) sscanf(line, "%s", dummystring); else if(strcmp(dummystring, "FRACTION_HYDROPHOBIC_SASA_CUTOFF")==0) { sscanf(line, "%s %lf", dummystring, &FRACTION_HYDROPHOBIC_SASA_CUTOFF); } else if(strcmp(dummystring, "TRANSFER_FREE_ENERGY_DENSITY_CUTOFF")==0) { sscanf(line, "%s %lf", dummystring, &TRANSFER_FREE_ENERGY_DENSITY_CUTOFF); } else if(strcmp(dummystring, "CORE_INTERFACIAL_BORN_RADIUS_CUTOFF")==0) { sscanf(line, "%s %lf", dummystring, &CORE_INTERFACIAL_BORN_RADIUS_CUTOFF); } else if(strcmp(dummystring, "INTERFACIAL_SURFACE_SASA_CUTOFF")==0) { sscanf(line, "%s %lf", dummystring, &INTERFACIAL_SURFACE_SASA_CUTOFF); } else if(strcmp(dummystring, "CORESURF_MEA_SASA")==0) { if(CORESURF_MEA_SASA == ENDFLAG) /* not defined in user input file */ sscanf(line, "%s %lf", dummystring, &CORESURF_MEA_SASA); } else if(strcmp(dummystring, "CORESURF_MEA_BORN")==0) { if(CORESURF_MEA_BORN == ENDFLAG) /* not defined in user input file */ sscanf(line, "%s %lf", dummystring, &CORESURF_MEA_SASA); } else if(strcmp(dummystring, "HYDROPHOBIC_ASP")==0) { if(HYDROPHOBIC_ASP == ENDFLAG) /* if !=ENDFLAG, then it was defined in the user input file */ sscanf(line, "%s %lf", dummystring, &HYDROPHOBIC_ASP); } else if(strcmp(dummystring, "NONBOND_FACTOR_1_4")==0) sscanf(line, "%s %lf", dummystring, &NONBOND_FACTOR_1_4); else if(strstr(dummystring, "UNSAT")!=0 && strstr(dummystring, "HBOND")!=0 && (strstr(dummystring, "REF")!=0 || strstr(dummystring, "SPECIF")!=0 ) ) { if(strstr(dummystring, "H_")!=0 || strstr(dummystring, "HYDROGEN")!=0) sscanf(line, "%s %lf", dummystring, &H_UNSAT_HBOND_REFSTATE); else if(strstr(dummystring, "N_")!=0 || strstr(dummystring, "DONOR")!=0) sscanf(line, "%s %lf", dummystring, &DONOR_UNSAT_HBOND_REFSTATE); else if(strstr(dummystring, "O_")!=0 || strstr(dummystring, "ACCEPTOR")!=0) sscanf(line, "%s %lf", dummystring, &ACCEPTOR_UNSAT_HBOND_REFSTATE); else if(UNSAT_HBOND_REFSTATE == ENDFLAG) sscanf(line, "%s %lf",dummystring, &UNSAT_HBOND_REFSTATE); /* else { sprintf(line,"ERROR Do not recognize %s in forcefield file",dummystring); failure_report(line,"exit"); } */ } else if(strcmp(dummystring, "GENERAL_ASP")==0) { if(GENERAL_ASP == ENDFLAG) sscanf(line, "%s %lf", dummystring, &GENERAL_ASP); } else if(strcmp(dummystring, "BORN_P1")==0) sscanf(line, "%s %lf", dummystring, &BORN_P1); else if(strcmp(dummystring, "BORN_P2")==0) sscanf(line, "%s %lf", dummystring, &BORN_P2); else if(strcmp(dummystring, "BORN_P3")==0) sscanf(line, "%s %lf", dummystring, &BORN_P3); else if(strcmp(dummystring, "BORN_P4")==0) sscanf(line, "%s %lf", dummystring, &BORN_P4); else if(strcmp(dummystring, "BORN_P5")==0) { sscanf(line, "%s %lf", dummystring, &BORN_P5); ONE_OVER_BORN_P5 = 1/BORN_P5; PI_TIMES_BORN_P5 = PI*BORN_P5; } else if(strcmp(dummystring, "DEBYE_SCALE_FACTOR")==0) { if(DEBYE_SCALE_FACTOR==0.000140) /* if it is still 0.00014, it was not defined in the inputfile */ sscanf(line, "%s %lf", dummystring, &DEBYE_SCALE_FACTOR); } else if(strcmp(dummystring, "BORN_LAMBDA")==0) { sscanf(line, "%s %lf", dummystring, &BORN_LAMBDA); ONE_OVER_BORN_LAMBDA = 1/BORN_LAMBDA; } else if(strcmp(dummystring, "VAR_FIX_BORN_SCALE")==0) sscanf(line, "%s %lf", dummystring, &VAR_FIX_BORN_SCALE); else if(strcmp(dummystring, "INTRA_ROTAMER_FLAG")==0) { if(INTRA_ROTAMER_FLAG==ENDFLAG) sscanf(line, "%s %d", dummystring, &INTRA_ROTAMER_FLAG); } else if(strcmp(dummystring, "MINIMIZE_ROTAMERS_FLAG")==0) { if(LOCAL_MINIMIZATION_FLAG==ENDFLAG) /* if not ENDFLAG, was user-specified */ { sscanf(line, "%s %d", dummystring, &LOCAL_MINIMIZATION_FLAG); if(LOCAL_MINIMIZATION_FLAG==1) LOCAL_MINIMIZATION_FLAG=2; } } else if(strcmp(dummystring, "VDW_CUTOFF")==0) { if(VDW_CUTOFF == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &VDW_CUTOFF); } else if(strstr(dummystring, "VDW_CLASH")!=0) { if(VDW_CUTOFF == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &VDW_CLASH_ENERGY); } else if(strcmp(dummystring, "WEIGHT_VDW")==0) { if(WEIGHT_VDW == ENDFLAG) sscanf(line, "%s %lf",dummystring, &WEIGHT_VDW); } else if(strcmp(dummystring, "WEIGHT_ELECTROSTATICS")==0) { if(WEIGHT_ELECTROSTATICS == ENDFLAG) sscanf(line, "%s %lf",dummystring, &WEIGHT_ELECTROSTATICS); } else if(strcmp(dummystring, "WEIGHT_1_4")==0 || strcmp(dummystring, "WEIGHT_TORSION")==0) { if(WEIGHT_1_4 == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &WEIGHT_1_4); } else if(strcmp(dummystring, "HBOND_FUNCTION_TYPE")==0) { sscanf(line, "%s %s",dummystring, dummystring); convert_string_to_all_caps(dummystring); if(strstr(dummystring,"10")!=0 && strstr(dummystring,"12")!=0) HBOND_FUNCTION_TYPE = 2; HBOND_FUNCTION_FLAG=1; } else if(strcmp(dummystring, "HBOND_FUNCTION_FLAG")==0) { if(HBOND_FUNCTION_FLAG == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %d",dummystring, &HBOND_FUNCTION_FLAG); } else if(strcmp(dummystring, "VDW_ATTRACTIVE_FACTOR")==0) { if(VDW_ATTRACTIVE_FACTOR == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &VDW_ATTRACTIVE_FACTOR); } else if(strcmp(dummystring, "VDW_REPULSIVE_FACTOR")==0) { if(VDW_REPULSIVE_FACTOR == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &VDW_REPULSIVE_FACTOR); } else if(strcmp(dummystring, "VDW_LINEAR_START_POINT")==0) { if(VDW_LINEAR_START_POINT == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &VDW_LINEAR_START_POINT); } else if(strcmp(dummystring, "OVERALL_ENERGY_SCALE")==0) { if(OVERALL_ENERGY_SCALE == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &OVERALL_ENERGY_SCALE); } else if(strstr(dummystring, "ENERGY")!=0 && (strstr(dummystring, "ERROR")!=0 || strstr(dummystring, "UNCERT")!=0 ) ) { if(ENERGY_ERROR == ENDFLAG) /* if it's not ENDFLAG, was defined in input file manually */ sscanf(line, "%s %lf",dummystring, &ENERGY_ERROR); } else if(strcmp(dummystring, "TORSION_FILE")==0) { torsionfilename = (char *)calloc(MAXLINE, sizeof(char)); sscanf(line, "%s %s",dummystring, torsionfilename); if(does_this_file_exist(torsionfilename)==0) { sprintf(dummystring,"%s%s",energy_dir,torsionfilename); strcpy(torsionfilename,dummystring); } } else if(strstr(dummystring, "RESPARAM")!=0) { resparamfilename = (char *)calloc(MAXLINE, sizeof(char)); sscanf(line, "%s %s",dummystring, resparamfilename); if(does_this_file_exist(resparamfilename)==0) { sprintf(dummystring,"%s%s",energy_dir,resparamfilename); strcpy(resparamfilename,dummystring); } } else if(strstr(dummystring, "ROTAMER")!=0 && strstr(dummystring, "LIB")!=0) { if(ROTAMERFILE==NULL) /* if not NULL, was defined in the inputfile, so ignore this */ { ROTAMERFILE = (char *)calloc(MAXLINE, sizeof(char)); sscanf(line, "%s %s",dummystring, ROTAMERFILE); } } else if(strcmp(dummystring, "SASA_POINTS_FILE")==0) { SASA_points_file= (char *)calloc(MAXLINE, sizeof(char)); sscanf(line, "%s %s",dummystring, SASA_points_file); if(does_this_file_exist(SASA_points_file)==0) { sprintf(dummystring,"%s%s",energy_dir,SASA_points_file); strcpy(SASA_points_file,dummystring); } } else if(word_count(line)!=0) // if 0, is a comment line { sprintf(line,"ERROR Do not recognize %s",dummystring); fprintf(stderr,"%s\n",line); exit(1); } } fclose(input); /* set some variables to default values if they are not specified here or in the inputfile */ if(INTRA_ROTAMER_FLAG==1 || INTRA_ROTAMER_FLAG==ENDFLAG) { INTRA_ROTAMER_FLAG=2; INTRA_ROTAMER_FLAG_14=1; } else { INTRA_ROTAMER_FLAG=3; INTRA_ROTAMER_FLAG_14=3; } if(LOCAL_MINIMIZATION_FLAG==ENDFLAG) LOCAL_MINIMIZATION_FLAG=0; if(CORESURF_MEA_SASA == ENDFLAG) CORESURF_MEA_SASA = 26.0; if(CORESURF_MEA_BORN == ENDFLAG) CORESURF_MEA_BORN = 3.0; if(VDW_CUTOFF == ENDFLAG) VDW_CUTOFF = 1000.0; if(VDW_CLASH_ENERGY == ENDFLAG) VDW_CLASH_ENERGY = 2.75; if(UNSAT_HBOND_REFSTATE == ENDFLAG) UNSAT_HBOND_REFSTATE = -1.0*VDW_CUTOFF; if(DONOR_UNSAT_HBOND_REFSTATE == ENDFLAG) DONOR_UNSAT_HBOND_REFSTATE = UNSAT_HBOND_REFSTATE; if(H_UNSAT_HBOND_REFSTATE == ENDFLAG) H_UNSAT_HBOND_REFSTATE = UNSAT_HBOND_REFSTATE; if(ACCEPTOR_UNSAT_HBOND_REFSTATE == ENDFLAG) ACCEPTOR_UNSAT_HBOND_REFSTATE = UNSAT_HBOND_REFSTATE; if(WEIGHT_VDW == ENDFLAG) WEIGHT_VDW = 1.0; if(WEIGHT_ELECTROSTATICS == ENDFLAG) WEIGHT_ELECTROSTATICS = 1.0; if(WEIGHT_1_4 == ENDFLAG) WEIGHT_1_4 = 1.0; if(VDW_ATTRACTIVE_FACTOR == ENDFLAG) VDW_ATTRACTIVE_FACTOR = 1.0; if(VDW_REPULSIVE_FACTOR == ENDFLAG) VDW_REPULSIVE_FACTOR = 1.0; if(OVERALL_ENERGY_SCALE == ENDFLAG) OVERALL_ENERGY_SCALE = 1.0; if(HYDROPHOBIC_ASP == ENDFLAG) HYDROPHOBIC_ASP = 0; if(GENERAL_ASP == ENDFLAG) GENERAL_ASP = 0; if(HYDROPHOBIC_ASP == 0 && GENERAL_ASP == 0) SASA_FLAG = 0; if(VDW_RADII_SCALE == 0) VDW_RADII_SCALE = 1.0; if(INTERNAL_DIELECTRIC == 0) { fprintf(stderr,"ERROR Need to define Ep in forcefield file or input file\n"); exit(1); } if(HBOND_DIELECTRIC == 0) HBOND_DIELECTRIC = 1.0; if(HBOND_FUNCTION_FLAG == ENDFLAG) HBOND_FUNCTION_FLAG = 0; /* constants for electrostatics calculations */ ONE_OVER_WATER_DIELECTRIC = 1.0/WATER_DIELECTRIC; KAPPA = -sqrt(DEBYE_SCALE_FACTOR*DEBYE_CONST*IONIC_STRENGTH/TEMPERATURE); RT = R_univ*TEMPERATURE; ONE_OVER_INTERNAL_DIELECTRIC = 1.0/INTERNAL_DIELECTRIC; BORN_TAU = ONE_OVER_WATER_DIELECTRIC - ONE_OVER_INTERNAL_DIELECTRIC; COULOMB_CONST_OVER_INTERNAL_DIELECTRIC = WEIGHT_ELECTROSTATICS*COULOMB_CONST/INTERNAL_DIELECTRIC; COUL_CONST_BORN_TAU = WEIGHT_ELECTROSTATICS*COULOMB_CONST*BORN_TAU; HALF_COUL_CONST_BORN_TAU = COUL_CONST_BORN_TAU/2.0; COULOMB_CONST_HBOND_INTERNAL_DIELECTRIC = COULOMB_CONST*(1.0/HBOND_DIELECTRIC-ONE_OVER_INTERNAL_DIELECTRIC); /* allocate memory for and precalculate vdw welldepths and sigmas for all possible atompairs */ beta=0; if(VDW_LINEAR_START_POINT>=10) VDW_LINEAR_START_POINT=ENDFLAG; if(ENERGY_ERROR == ENDFLAG) ENERGY_ERROR = RT; if(VDW_WELLDEPTH==NULL) { VDW_WELLDEPTH = (double **)calloc(MAXATOMTYPES, sizeof(double *)); VDW_WELLDEPTH_14 = (double **)calloc(MAXATOMTYPES, sizeof(double *)); if(VDW_LINEAR_START_POINT!=ENDFLAG) { beta = 1- VDW_LINEAR_START_POINT*( pow(2.0,(1.0/6.0)) - 1 ); VDW_LINEAR_SLOPE = (double **)calloc(MAXATOMTYPES, sizeof(double *)); VDW_LINEAR_INTERCEPT = (double **)calloc(MAXATOMTYPES, sizeof(double *)); VDW_LINEAR_SWITCHPOINT = (double **)calloc(MAXATOMTYPES, sizeof(double *)); } VDW_CONTACT_DISTANCE = (double **)calloc(MAXATOMTYPES, sizeof(double *)); VDW_SIGMA = (double **)calloc(MAXATOMTYPES, sizeof(double *)); VDW_SIGMA6_14 = (double **)calloc(MAXATOMTYPES, sizeof(double *)); VDW_SIGMA6_SCALE = (double **)calloc(MAXATOMTYPES, sizeof(double *)); for(i=0;iatom_ptr = &atomparam[i]; U_ATOMRESPARAM->donorflag[0]=' '; U_ATOMRESPARAM->E_hbond_specificity=0; /* read SASA_points_file; lookuptable for SASA calcs */ SASA_readfiles(SASA_points_file, atomparam); if(SASA_points_file!=NULL) free_memory(SASA_points_file); if(torsionfilename!=NULL) { read_torsiondata(torsionfilename, atomparam); free_memory(torsionfilename); torsionfilename=NULL; } else { fprintf(stderr,"ERROR TORSION_FILE not defined!\n"); exit(1); } if(RESPARAMFILE == NULL) /* resparam file not user-specified; use the one listed in the forcefield file */ { if(resparamfilename!=NULL) { /* determine MAX_RES_SIZE */ read_residuedata(resparamfilename, resparam, atomparam,1); if(LIGANDFILENAME != NULL) /* ligands! */ { i=1; while(LIGANDFILENAME[i]!=NULL) { read_residuedata(LIGANDFILENAME[i], resparam, atomparam,1); /* determine MAX_RES_SIZE */ ++i; } } /* now, actually read residue data */ read_residuedata(resparamfilename, resparam, atomparam, 0); free_memory(resparamfilename); } else { fprintf(stderr,"ERROR RESPARAM not defined!\n"); exit(1); } } else /* user-defined RESPARAM file */ { /* determine MAX_RES_SIZE */ read_residuedata(RESPARAMFILE, resparam, atomparam,1); /* determine MAX_RES_SIZE */ if(LIGANDFILENAME != NULL) /* ligands! */ { i=1; while(LIGANDFILENAME[i]!=NULL) { read_residuedata(LIGANDFILENAME[i], resparam, atomparam,1); /* determine MAX_RES_SIZE */ ++i; } } /* now, actually read residue data */ read_residuedata(RESPARAMFILE, resparam, atomparam, 0); free_memory(RESPARAMFILE); RESPARAMFILE=NULL; if(resparamfilename!=NULL) free_memory(resparamfilename); } if(ROTAMERFILE == NULL) { failure_report("ERROR ROTAMER_LIBRARY not defined","exit"); } else { if(does_this_file_exist(ROTAMERFILE)==0) { sprintf(dummystring,"%s%s",energy_dir,ROTAMERFILE); strcpy(ROTAMERFILE,dummystring); } read_rotamers(ROTAMERFILE, rotamerlib); free_memory(ROTAMERFILE); ROTAMERFILE=NULL; } /* read ligand data */ if(LIGANDFILENAME != NULL) { i=1; while(LIGANDFILENAME[i]!=NULL) { read_ligand_data(LIGANDFILENAME[i], resparam, rotamerlib, atomparam); ++i; } /* free memory */ i=1; while(LIGANDFILENAME[i]!=NULL) { free_memory(LIGANDFILENAME[i]); ++i; } } /* sort resparam and rotamerlib alphabetically */ qsort(resparam,MAXRESTYPES,sizeof(RESPARAM),restypecompare_resparam); qsort(rotamerlib,MAXRESTYPES,sizeof(ROTAMERLIB),restypecompare_rotamerlib); free_memory(line); free_memory(dummystring); free_memory(energy_dir); /* set up pointers between rotamerlib and resparam and check for inconsistencies */ i=1; while(strcmp(resparam[i].residuetype,"zzz")!=0) { if(strcmp(resparam[i].residuetype,"MEA")==0) MEA = &resparam[i]; j=1; while(strcmp(resparam[i].residuetype, rotamerlib[j].residuetype)!=0 && strcmp(rotamerlib[j].residuetype, "zzz")!=0) ++j; if(strcmp(rotamerlib[j].residuetype, "zzz")!=0) /* found rotamerlib */ { if(abs(resparam[i].residuecode)<100) /* sidechains have abs(resparam[i].residuecode)<100 */ { if(rotamerlib[j].numChi < resparam[i].num_rotatable_bonds) { fprintf(stderr,"ERROR %d chi-angles defined for %s, which has %d rotatable bonds\n", rotamerlib[j].numChi, rotamerlib[j].residuetype, resparam[i].num_rotatable_bonds); exit(1); } resparam[i].rotamerlib_ptr = &rotamerlib[j]; rotamerlib[j].residuecode = resparam[i].residuecode; } else /* backbone or something else fixed; no rotamerlib, so need a virtual one for numChi */ { resparam[i].rotamerlib_ptr = (ROTAMERLIB *)malloc(sizeof(ROTAMERLIB)); resparam[i].rotamerlib_ptr->numChi = 0; } } else /* can't find rotamerlib */ { resparam[i].rotamerlib_ptr = NULL; if(abs(resparam[i].residuecode)<100) /* this is sidechain; can't find rotamerlib is fatal */ { fprintf(stderr, "ERROR Cannot find rotamers for %s\n", resparam[i].residuetype); exit(1); } else /* backbone or something else fixed; no rotamerlib, so need a virtual one for numChi */ { resparam[i].rotamerlib_ptr = (ROTAMERLIB *)malloc(sizeof(ROTAMERLIB)); resparam[i].rotamerlib_ptr->numChi = 0; } } ++i; } volume_bonded(resparam); /* calculate atomic volumes for born radius calculations */ /* adjust radii of zero-radius hydrogens per Schaefer & Karplus J Phys Chem 1996 */ i=1; while(strcmp(resparam[i].residuetype,"zzz")!=0) { for(j=4;j<=(resparam[i].numberofAtoms+3);++j) { for(k=0;kunited_radius == 0 && strstr(resparam[i].atom[j].atomtype, "U")==0) { resparam[i].atom[j].atom_ptr->united_radius = resparam[i].atom[resparam[i].atom[j].contactAtomNumber].atom_ptr->united_radius - resparam[i].atom[j].bondlength; } } } ++i; } /* set up UNITED_RADII_SUM_SQUARED matrix for born radius calcs */ if(UNITED_RADII_SUM_SQUARED==NULL) { UNITED_RADII_SUM_SQUARED = (double **)calloc(MAXATOMTYPES, sizeof(double *)); for(i=0;i