Mercurial > repos > guerler > springsuite
view spring_package/pulchra/pulchra.c @ 17:c790d25086dc draft
"planemo upload commit b0ede77caf410ab69043d33a44e190054024d340-dirty"
| author | guerler |
|---|---|
| date | Wed, 28 Oct 2020 05:11:56 +0000 |
| parents | |
| children |
line wrap: on
line source
// // PULCHRA // Protein Chain Restoration Algorithm // // Version 3.04 // July 2007 // Contact: Piotr Rotkiewicz, piotr -at- pirx -dot- com // // to compile: // cc -O3 pulchra.c pulchra_data.c -lm -o pulchra // // to use: // ./pulchra file.pdb // // to display available options: // ./pulchra // #define COMPILE_BB #define COMPILE_ROT #include <math.h> #include <stdio.h> #include <stdlib.h> #include <string.h> #include <sys/timeb.h> #include <time.h> #define uchar unsigned char #define uint unsigned int #define real double #include "pulchra_common.h" #define PULCHRA_VERSION 3.04 #define MAX_BUF_SIZE 1000 #define FILE_SUCCESS 0 #define FILE_NOT_FOUND -1 #define FILE_WARNING -2 #define FATAL_MAE -1 #define FLAG_BACKBONE 1 #define FLAG_CALPHA 2 #define FLAG_SIDECHAIN 4 #define FLAG_SCM 8 #define FLAG_PROTEIN 1 #define FLAG_DNA 2 #define FLAG_RNA 4 #define FLAG_CHYDRO 8 #define RADDEG 180./M_PI #define DEGRAD M_PI/180. int _VERBOSE = 0; int _BB_REARRANGE = 1; int _BB_OPTIMIZE = 0; int _CA_OPTIMIZE = 1; int _CA_RANDOM = 0; int _CA_ITER = 100; int _CA_TRAJECTORY = 0; int _CISPRO = 0; int _CHIRAL = 1; int _CENTER_CHAIN = 0; int _REBUILD_BB = 1; int _REBUILD_SC = 1; int _REBUILD_H = 0; int _PDB_SG = 0; int _TIME_SEED = 0; int _XVOLUME = 1; int _XVOL_ITER = 3; real _CA_START_DIST = 3.0; real _CA_XVOL_DIST = 3.5; real _SG_XVOL_DIST = 1.6; #define CALC_C_ALPHA #define CALC_C_ALPHA_ANGLES #define CALC_C_ALPHA_START #define CALC_C_ALPHA_XVOL real CA_K=10.0; real CA_ANGLE_K=20.0; real CA_START_K=0.01; real CA_XVOL_K=10.00; #define CA_DIST 3.8 #define CA_DIST_TOL 0.1 #define CA_DIST_CISPRO 2.9 #define CA_DIST_CISPRO_TOL 0.1 #define E_EPS 1e-10 // grid resolution (used for fast clash detection) #define GRID_RES 6.0 int chain_length = 0; char AA_NAMES[21][4] = { "GLY", "ALA", "SER", "CYS", "VAL", "THR", "ILE", "PRO", "MET", "ASP", "ASN", "LEU", "LYS", "GLU", "GLN", "ARG", "HIS", "PHE", "TYR", "TRP", "UNK" }; char SHORT_AA_NAMES[22] = { "GASCVTIPMDNLKEQRHFYWX" }; int AA_NUMS[256]; int nheavy[20] = { 0, 1, 2, 2, 3, 3, 4, 3, 4, 4, 4, 4, 5, 5, 5, 7, 6, 7, 8, 10}; char *backbone_atoms[4] = { "N ", "CA ", "C ", "O " }; char *heavy_atoms[200]= { /* GLY */ NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* ALA */ "CB ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* SER */ "CB ", "OG ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* CYS */ "CB ", "SG ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* VAL */ "CB ", "CG1", "CG2", NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* THR */ "CB ", "OG1", "CG2", NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* ILE */ "CB ", "CG1", "CG2", "CD1", NULL, NULL, NULL, NULL, NULL, NULL, /* PRO */ "CB ", "CG ", "CD ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, /* MET */ "CB ", "CG ", "SD ", "CE ", NULL, NULL, NULL, NULL, NULL, NULL, /* ASP */ "CB ", "CG ", "OD1", "OD2", NULL, NULL, NULL, NULL, NULL, NULL, /* ASN */ "CB ", "CG ", "OD1", "ND2", NULL, NULL, NULL, NULL, NULL, NULL, /* LEU */ "CB ", "CG ", "CD1", "CD2", NULL, NULL, NULL, NULL, NULL, NULL, /* LYS */ "CB ", "CG ", "CD ", "CE ", "NZ ", NULL, NULL, NULL, NULL, NULL, /* GLU */ "CB ", "CG ", "CD ", "OE1", "OE2", NULL, NULL, NULL, NULL, NULL, /* GLN */ "CB ", "CG ", "CD ", "OE1", "NE2", NULL, NULL, NULL, NULL, NULL, /* ARG */ "CB ", "CG ", "CD ", "NE ", "CZ ", "NH1", "NH2", NULL, NULL, NULL, /* HIS */ "CB ", "CG ", "ND1", "CD2", "CE1", "NE2", NULL, NULL, NULL, NULL, /* PHE */ "CB ", "CG ", "CD1", "CD2", "CE1", "CE2", "CZ ", NULL, NULL, NULL, /* TYR */ "CB ", "CG ", "CD1", "CD2", "CE1", "CE2", "CZ ", "OH ", NULL, NULL, /* TRP */ "CB ", "CG ", "CD1", "CD2", "NE1", "CE2", "CE3", "CZ2", "CZ3", "CH2"}; /* reads full-atom pdb file */ struct _res_type; typedef struct _atom_type { struct _atom_type *next; real x, y, z; char *name; int num, locnum; int flag; char cispro; int gx, gy, gz; struct _res_type *res; struct _atom_type *prev; } atom_type; typedef struct _res_type { struct _res_type *next; atom_type *atoms; int num, locnum, natoms; int type; char pdbsg; char protein; char *name; char chain; real sgx, sgy, sgz; real cmx, cmy, cmz; struct _res_type *prev; } res_type; typedef struct _mol_type { struct _mol_type *next; res_type *residua; int nres; unsigned char *r14; char *name; uchar *seq; char **contacts; real **cutoffs; struct _mol_type *prev; } mol_type; #define MIN(a,b) (a<b?a:b) #define MAX(a,b) (a>b?a:b) mol_type *chain = NULL; real rnd(void) { return 0.001*(real)(rand()%1000); } atom_type *new_atom(void) { atom_type *tmpatom; tmpatom = (atom_type*) calloc(sizeof(atom_type),1); if (tmpatom) { tmpatom->x=tmpatom->y=tmpatom->z=0.; tmpatom->name=NULL; tmpatom->num=tmpatom->locnum=tmpatom->flag=0; tmpatom->next=tmpatom->prev=NULL; } return tmpatom; } res_type* new_res(void) { res_type *tmpres; tmpres = (res_type*) calloc(sizeof(res_type),1); if (tmpres) { tmpres->num=0; tmpres->name=NULL; tmpres->atoms=NULL; tmpres->chain=' '; tmpres->next=tmpres->prev=NULL; } return tmpres; } mol_type *new_mol(void) { mol_type *tmpmol; tmpmol = (mol_type*) calloc(sizeof(mol_type),1); if (tmpmol) { tmpmol->name=NULL; tmpmol->residua=NULL; tmpmol->next=tmpmol->prev=NULL; } return tmpmol; } void add_atom(atom_type* atomlist, atom_type* newatom) { atom_type *tmpatom; if (!atomlist) atomlist=newatom; else { tmpatom=atomlist->next; atomlist->next=newatom; newatom->prev=atomlist; newatom->next=tmpatom; if (tmpatom) tmpatom->prev=newatom; } } void add_res(res_type* reslist, res_type* newres) { res_type *tmpres; if (!reslist) reslist=newres; else { tmpres=reslist->next; reslist->next=newres; newres->prev=reslist; newres->next=tmpres; if (tmpres) tmpres->prev=newres; } } void add_mol(mol_type* mollist, mol_type* newmol) { mol_type *tmpmol; if (!mollist) mollist=newmol; else { tmpmol=mollist->next; mollist->next=newmol; newmol->prev=mollist; newmol->next=tmpmol; if (tmpmol) tmpmol->prev=newmol; } } void delete_atom(atom_type* atom) { atom_type *tmpatom; if (atom->prev) atom->prev->next=atom->next; if (atom->next) atom->next->prev=atom->prev; if (atom->name) free(atom->name); free(atom); atom=NULL; } void delete_res(res_type* res) { res_type *tmpres; atom_type *tmpatom; if (res->prev) res->prev->next=res->next; if (res->next) res->next->prev=res->prev; if (res->name) free(res->name); if (res->atoms) { while (res->atoms) { tmpatom = res->atoms->next; delete_atom(res->atoms); res->atoms=tmpatom; } } free(res); res=NULL; } void delete_mol(mol_type* mol) { mol_type *tmpmol; res_type *tmpres; int i; if (mol->prev) mol->prev->next=mol->next; if (mol->next) mol->next->prev=mol->prev; if (mol->name) free(mol->name); if (mol->residua) { while (mol->residua) { tmpres = mol->residua->next; delete_res(mol->residua); mol->residua=tmpres; } } if (mol->contacts) { for (i=0; i<mol->nres; i++) free(mol->contacts[i]); free(mol->contacts); } if (mol->cutoffs) { for (i=0; i<mol->nres; i++) free(mol->cutoffs[i]); free(mol->cutoffs); } free(mol); mol=NULL; } atom_type* get_last_atom(atom_type* atom) { while (atom->next) atom=atom->next; return atom; } res_type* get_last_res(res_type* res) { while (res->next) res=res->next; return res; } mol_type *get_last_mol(mol_type* mol) { while (mol->next) mol=mol->next; return mol; } char setseq(char* aaname) { int i; for (i=0; i<21; i++) { if ((aaname[0]==AA_NAMES[i][0]) && (aaname[1]==AA_NAMES[i][1]) && (aaname[2]==AA_NAMES[i][2])) break; } if (i==21) { if (!strcmp(aaname, "GLX")) return 'E'; if (!strcmp(aaname, "ASX")) return 'D'; if (!strcmp(aaname, "HID")) return 'H'; if (!strcmp(aaname, "MSE")) return 'M'; if (!strcmp(aaname, "SEP")) return 'S'; if (!strcmp(aaname, "TPO")) return 'T'; if (!strcmp(aaname, "PTR")) return 'Y'; i--; } return SHORT_AA_NAMES[i]; } int orient(res_type *res1, res_type *res2) { real x1, y1, z1; real x2, y2, z2; real cax, cay, caz; real len, vect, angle; atom_type *atom; if (!res1 || !res2) return 0; atom=res1->atoms; cax=cay=caz=0.; while (atom) { if (!strncmp(atom->name,"CA",2)) { cax=atom->x; cay=atom->y; caz=atom->z; } atom=atom->next; } x1=res1->sgx-cax; y1=res1->sgy-cay; z1=res1->sgz-caz; if (x1==0. && y1==0. && z1==0.) x1+=1.0; atom=res2->atoms; cax=cay=caz=0.; while (atom) { if (!strncmp(atom->name,"CA",2)) { cax=atom->x; cay=atom->y; caz=atom->z; } atom=atom->next; } x2=res2->sgx-cax; y2=res2->sgy-cay; z2=res2->sgz-caz; if (x2==0. && y2==0. && z2==0.) x2+=1.0; vect = x1*x2+y1*y2+z1*z2; len = sqrt(x1*x1+y1*y1+z1*z1)*sqrt(x2*x2+y2*y2+z2*z2); if (len) vect /= len; angle=RADDEG*acos(vect); if (angle>120.) return 1; /*anti*/ if (angle>60.) return 2; /*mid*/ return 3; /*par*/ } int res_contact(res_type *res1, res_type *res2) { atom_type *atoms1, *atoms2; real dx, dy, dz; atoms1 = res1->atoms; while (atoms1) { atoms2 = res2->atoms; while (atoms2) { dx=atoms1->x-atoms2->x; dy=atoms1->y-atoms2->y; dz=atoms1->z-atoms2->z; if ((atoms1->flag & FLAG_SIDECHAIN) && (atoms2->flag & FLAG_SIDECHAIN) && (dx*dx+dy*dy+dz*dz<20.25)) { return 1; } atoms2=atoms2->next; } atoms1=atoms1->next; } return 0; } int read_pdb_file(char* filename, mol_type* molecules, char *realname) { FILE *inp; char buffer[1000]; char atmname[10]; char resname[10]; char version; int prevresnum, resnum, atmnum, locatmnum, num, locnum=0, i, j; atom_type *prevatom1, *prevatom2, *prevatom3, *prevatom4; int sgnum, cc, nres, ok, natom; real sgx, sgy, sgz; res_type *res, *test1, *test2; atom_type *atom; real x, y, z; real dist; unsigned char bin; int warn=0; real cutoff; if (_VERBOSE) printf("Reading input file %s...\n", filename); inp = fopen(filename, "r"); if (!inp) { if (_VERBOSE) printf("ERROR: can't open %s !!!\n", filename); return FILE_NOT_FOUND; } molecules->nres=0; molecules->name=(char*)calloc(strlen(realname)+1,1); strcpy(molecules->name, realname); atmname[3]=0; resname[3]=0; prevresnum=-666; locatmnum=0; sgnum=0; sgx=sgy=sgz=0.; res=NULL; while (!feof(inp)) { if (fgets(buffer, 1000, inp)!=buffer) break; if (!strncmp(buffer, "END", 3) || !strncmp(buffer, "TER", 3)) break; // end of file; only single molecule read if (!strncmp(buffer, "ATOM", 4) || !strncmp(buffer, "HETATM", 6)) { if (buffer[16]!=' ' && buffer[16]!='A') continue; sscanf(&buffer[22], "%d", &resnum); strncpy(resname, &buffer[17], 3); strncpy(atmname, &buffer[13], 3); if (resnum==prevresnum && !strncmp(atmname, "N ", 2)) { if (_VERBOSE) printf("WARNING: fault in numeration at residuum %s[%d]\n", resname, resnum); warn=1; } if (atmname[0]=='H') continue; if (resnum!=prevresnum || !strncmp(atmname, "N ", 2)) { prevresnum=resnum; if (res) { if (sgnum) { res->sgx=sgx/(real)sgnum; res->sgy=sgy/(real)sgnum; res->sgz=sgz/(real)sgnum; } else { res->sgx=res->sgy=res->sgz=0.; } } locatmnum=0; version=' '; res = new_res(); sgnum=0; sgx=sgy=sgz=0.; molecules->nres++; res->name = calloc(strlen(resname)+1, 1); res->type = AA_NUMS[setseq(resname)]; res->locnum=locnum++; res->num = resnum; res->natoms=0; res->chain = buffer[21]; strcpy(res->name, resname); if (molecules->residua) { add_res(get_last_res(molecules->residua), res); } else { molecules->residua = res; } } atom = new_atom(); atom->res = res; res->natoms++; locatmnum++; sscanf(&buffer[7], "%d", &atmnum); sscanf(&buffer[30], "%lf", &x); sscanf(&buffer[38], "%lf", &y); sscanf(&buffer[46], "%lf", &z); version = buffer[16]; atom->name = calloc(strlen(atmname)+1,1); strcpy(atom->name, atmname); atom->x=x; atom->y=y; atom->z=z; atom->num = atmnum; atom->locnum = locatmnum; if ((atmname[0]=='S' && atmname[1]=='C')||(atmname[0]=='C' && atmname[1]=='M')) { res->cmx = x; res->cmy = y; res->cmz = z; res->pdbsg=1; if (res->type<20) { res->protein=1; } } else if (!( ((atmname[0]=='C' || atmname[0]=='N' || atmname[0]=='O') && atmname[1]==' ') || (atmname[0]=='H') || (atmname[0]=='C' && atmname[1]=='A') || (atmname[0]=='O' && atmname[1]=='X' && atmname[2]=='T') ) ) { sgx+=x; sgy+=y; sgz+=z; sgnum++; atom->flag |= FLAG_SIDECHAIN; } else atom->flag |= FLAG_BACKBONE; if (atmname[0]=='C' && atmname[1]=='A') { atom->flag |= FLAG_BACKBONE; if (res->type<20) { res->protein=1; } if (!res->pdbsg) { res->cmx = x; res->cmy = y; res->cmz = z; } } if (atmname[0]=='C' && atmname[1]=='M') { atom->flag |= FLAG_SCM; } if (atmname[0]=='S' && atmname[1]=='C') { atom->flag |= FLAG_SCM; } if (res->atoms) { add_atom(get_last_atom(res->atoms), atom); } else { res->atoms = atom; } } } if (res) { if (sgnum) { res->sgx=sgx/(real)sgnum; res->sgy=sgy/(real)sgnum; res->sgz=sgz/(real)sgnum; } else { res->sgx=res->sgy=res->sgz=0.; } } fclose(inp); molecules->seq = (uchar*)calloc(sizeof(uchar)*molecules->nres+1,1); res=molecules->residua; i=0; while (res) { molecules->seq[i++]=(uchar)AA_NUMS[(int)setseq(res->name)]; res=res->next; } if (!warn) return FILE_SUCCESS; else return FILE_WARNING; } #define bool int #define true 1 #define false 0 real calc_ca_energy(atom_type **c_alpha, real **new_c_alpha, real **init_c_alpha, real **gradient, real alpha, real *ene, bool calc_gradient) { int i, j; real dx, dy, dz; real dist, ddist, ddist2; real new_e_pot; real theta0, tdif, th, aa, bb, ab; real ff0, ff2, dth, m0, m2, grad, f0[3], f2[3]; real adiff[3], bdiff[3]; real deriv, theta, dtheta, len1, len2, cos_theta, sin_theta; real dx1, dy1, dz1; real dx2, dy2, dz2; real dx3, dy3, dz3; real vx1, vy1, vz1; real vx2, vy2, vz2; real vx3, vy3, vz3; real r12x, r12y, r12z; real r32x, r32y, r32z; real d12, d32, d12inv, d32inv, c1, c2, diff; real f1x, f1y, f1z; real f2x, f2y, f2z; real f3x, f3y, f3z; for (i=0; i<chain_length; i++) { new_c_alpha[i][0]=c_alpha[i]->x+alpha*gradient[i][0]; new_c_alpha[i][1]=c_alpha[i]->y+alpha*gradient[i][1]; new_c_alpha[i][2]=c_alpha[i]->z+alpha*gradient[i][2]; } new_e_pot = 0.0; ene[0]=ene[1]=ene[2]=ene[3]=0.0; for (i=0; i<chain_length; i++) { #ifdef CALC_C_ALPHA_START dx=new_c_alpha[i][0]-init_c_alpha[i][0]; dy=new_c_alpha[i][1]-init_c_alpha[i][1]; dz=new_c_alpha[i][2]-init_c_alpha[i][2]; dist=sqrt(dx*dx+dy*dy+dz*dz); ddist = -dist; if (dist>_CA_START_DIST) { ddist2=dist*dist; new_e_pot+=CA_START_K*ddist2; ene[1] += CA_START_K*ddist2; if (calc_gradient) { grad = ddist * (-2.0*CA_START_K)/dist; gradient[i][0]-=grad*dx; gradient[i][1]-=grad*dy; gradient[i][2]-=grad*dz; } } #endif #ifdef CALC_C_ALPHA if (i>0) { dx=new_c_alpha[i][0]-new_c_alpha[i-1][0]; dy=new_c_alpha[i][1]-new_c_alpha[i-1][1]; dz=new_c_alpha[i][2]-new_c_alpha[i-1][2]; dist=sqrt(dx*dx+dy*dy+dz*dz); if (c_alpha[i]->cispro) { ddist=CA_DIST_CISPRO-dist; // if (fabs(ddist)<CA_DIST_CISPRO_TOL) ddist=0.0; } else { ddist=CA_DIST-dist; // if (fabs(ddist)<CA_DIST_TOL) ddist=0.0; } ddist2=ddist*ddist; new_e_pot+=CA_K*ddist2; ene[0] += CA_K*ddist2; if (calc_gradient) { grad = ddist * (-2.0*CA_K)/dist; gradient[i][0]-=grad*dx; gradient[i][1]-=grad*dy; gradient[i][2]-=grad*dz; gradient[i-1][0]+=grad*dx; gradient[i-1][1]+=grad*dy; gradient[i-1][2]+=grad*dz; } } #endif #ifdef CALC_C_ALPHA_XVOL for (j=0;j<i;j++) { if (abs(i-j)>2) { dx=new_c_alpha[i][0]-new_c_alpha[j][0]; dy=new_c_alpha[i][1]-new_c_alpha[j][1]; dz=new_c_alpha[i][2]-new_c_alpha[j][2]; dist=sqrt(dx*dx+dy*dy+dz*dz); ddist = dist-_CA_XVOL_DIST; if (dist<_CA_XVOL_DIST) { ddist2 = dist*dist; new_e_pot+=CA_XVOL_K*ddist2; ene[3] += CA_XVOL_K*ddist2; if (calc_gradient) { grad = ddist*(8.0*CA_XVOL_K)/dist; gradient[i][0]-=grad*dx; gradient[i][1]-=grad*dy; gradient[i][2]-=grad*dz; gradient[j][0]+=grad*dx; gradient[j][1]+=grad*dy; gradient[j][2]+=grad*dz; } } } } #endif #ifdef CALC_C_ALPHA_ANGLES if (i>0 && i<chain_length-1) { r12x=new_c_alpha[i-1][0]-new_c_alpha[i][0]; r12y=new_c_alpha[i-1][1]-new_c_alpha[i][1]; r12z=new_c_alpha[i-1][2]-new_c_alpha[i][2]; r32x=new_c_alpha[i+1][0]-new_c_alpha[i][0]; r32y=new_c_alpha[i+1][1]-new_c_alpha[i][1]; r32z=new_c_alpha[i+1][2]-new_c_alpha[i][2]; d12 = sqrt(r12x*r12x+r12y*r12y+r12z*r12z); d32 = sqrt(r32x*r32x+r32y*r32y+r32z*r32z); cos_theta = (r12x*r32x+r12y*r32y+r12z*r32z)/(d12*d32); if (cos_theta>1.0) cos_theta = 1.0; else if (cos_theta<-1.0) cos_theta = -1.0; sin_theta = sqrt(1.0-cos_theta*cos_theta); theta = acos(cos_theta); if (RADDEG*theta<80.) diff = theta-80.*DEGRAD; else if (RADDEG*theta>150.) diff = theta-150.*DEGRAD; else diff = 0.0; new_e_pot += CA_ANGLE_K*diff*diff; ene[2] += CA_ANGLE_K*diff*diff; if (calc_gradient) { d12inv = 1.0/d12; d32inv = 1.0/d32; diff *= (-2.0*CA_ANGLE_K)/sin_theta; c1 = diff*d12inv; c2 = diff*d32inv; f1x = c1*(r12x*(d12inv*cos_theta)-r32x*d32inv); f1y = c1*(r12y*(d12inv*cos_theta)-r32y*d32inv); f1z = c1*(r12z*(d12inv*cos_theta)-r32z*d32inv); f3x = c2*(r32x*(d32inv*cos_theta)-r12x*d12inv); f3y = c2*(r32y*(d32inv*cos_theta)-r12y*d12inv); f3z = c2*(r32z*(d32inv*cos_theta)-r12z*d12inv); f2x = -f1x-f3x; f2y = -f1y-f3y; f2z = -f1z-f3z; gradient[i-1][0]+=f1x; gradient[i-1][1]+=f1y; gradient[i-1][2]+=f1z; gradient[i][0]+=f2x; gradient[i][1]+=f2y; gradient[i][2]+=f2z; gradient[i+1][0]+=f3x; gradient[i+1][1]+=f3y; gradient[i+1][2]+=f3z; } } #endif } //printf("ene[3] = %f\n", ene[3]); return new_e_pot; } /* * Steepest gradient optimization using v=k*(r0-r)^2 * k = CA_K, r0 = CA_DIST */ void ca_optimize(char *tname, char *iname) { char buf[1000]; int i, j, hx, my_iter; real dx, dy, dz, dd, dist, dist2, dist3, ddist, ddist2; real e_pot, new_e_pot, grad, alpha, e_pot1, e_pot2, e_pot3; real adiff[3], bdiff[3]; real ff0, ff2, aa, ab, bb, th, tdif, dth, m0, m2; real theta0, deg_th, maxgrad, sum; real f0[3], f2[3]; real x, y, z; int numsteps, numsteps2, msteps; int *sec; real **new_c_alpha, **gradient, **init_c_alpha, last_alpha, tmp, last_good_alpha, d_alpha, last_e_pot; atom_type *atom, **c_alpha; res_type *res; FILE *inp, *out; int mnum, init, ok; real alpha1, alpha2, alpha3, a0; real ene1, ene2, ene3, e0; real energies[4]; real w1, w2, w3, eps; real gnorm, last_gnorm; int mode, fcnt; if (_CA_TRAJECTORY) { out = fopen(tname,"w"); if (out) fclose(out); } if (_VERBOSE) printf("Alpha carbons optimization...\n"); new_c_alpha = (real**)calloc(sizeof(real*)*(chain_length+1),1); init_c_alpha = (real**)calloc(sizeof(real*)*(chain_length+1),1); for (i=0;i<=chain_length;i++) { new_c_alpha[i] = (real*)calloc(sizeof(real)*3,1); init_c_alpha[i] = (real*)calloc(sizeof(real)*3,1); } gradient = (real**)calloc(sizeof(real*)*(chain_length+1),1); for (i=0;i<=chain_length;i++) { gradient[i] = (real*)calloc(sizeof(real)*3,1); } c_alpha = (atom_type**)calloc(sizeof(atom_type*)*(chain_length+1),1); i = 0; res = chain->residua; while (res) { atom = res->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]=='A') { if (i<chain_length) { c_alpha[i] = atom; i++; break; } else { if (_VERBOSE) printf("WARNING: number of C-alpha atoms exceeds the chain length!\n"); break; } } atom = atom->next; } res = res->next; } if (i<chain_length) chain_length = i; for (i=0; i<chain_length; i++) { init_c_alpha[i][0] = c_alpha[i]->x; init_c_alpha[i][1] = c_alpha[i]->y; init_c_alpha[i][2] = c_alpha[i]->z; } if (_CISPRO) { for (i=1; i<chain_length; i++) { dx = c_alpha[i]->x-c_alpha[i-1]->x; dy = c_alpha[i]->y-c_alpha[i-1]->y; dz = c_alpha[i]->z-c_alpha[i-1]->z; dd = sqrt(dx*dx+dy*dy+dz*dz); if ((setseq(c_alpha[i]->res->name)=='P') && (dd>CA_DIST_CISPRO-5*CA_DIST_CISPRO_TOL) && (dd<CA_DIST_CISPRO+5*CA_DIST_CISPRO_TOL)) { if (_VERBOSE) printf("Probable cis-proline found at postion %d\n", c_alpha[i]->res->num); c_alpha[i]->cispro = 1; } } } if (_CA_RANDOM) { if (_VERBOSE) printf("Generating random C-alpha coordinates...\n"); c_alpha[0]->x = 0.0; c_alpha[0]->y = 0.0; c_alpha[0]->z = 0.0; for (i=1;i<chain_length;i++) { dx = 0.01*(100-rand()%200); dy = 0.01*(100-rand()%200); dz = 0.01*(100-rand()%200); dd = 3.8/sqrt(dx*dx+dy*dy+dz*dz); dx *= dd; dy *= dd; dz *= dd; c_alpha[i]->x = c_alpha[i-1]->x+dx; c_alpha[i]->y = c_alpha[i-1]->y+dy; c_alpha[i]->z = c_alpha[i-1]->z+dz; } } if (iname) { inp = fopen(iname,"r"); if (inp) { if (_VERBOSE) printf("Reading initial structure %s...\n", iname); i = 0; while (!feof(inp)) { if (fgets(buf,1000,inp)==buf && buf[13]=='C' && buf[14]=='A') { if (i<chain_length) { if (sscanf(&buf[30],"%lf%lf%lf",&x,&y,&z)==3) { c_alpha[i]->x = x; c_alpha[i]->y = y; c_alpha[i]->z = z; i++; } } else { if (_VERBOSE) printf("WARNING: number of ini-file C-alpha atoms exceeds the chain length!\n"); break; } } } fclose(inp); } else if (_VERBOSE) printf("WARNING: can't read initial corrdinates %s\n", iname); } mnum = 1; mode = 0; init = 0; numsteps=numsteps2=0; last_alpha = 0.0; if (_VERBOSE) printf("Optimizing alpha carbons...\n"); eps = 0.5; fcnt=0; last_gnorm = 1000.; do { last_e_pot = e_pot; if (_CA_TRAJECTORY) { out = fopen(tname,"a"); if (out) { fprintf(out,"MODEL %d\n",mnum++); for (i=0; i<chain_length; i++) { fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n", i+1, "CA ", c_alpha[i]->res->name, ' ', c_alpha[i]->res->num, c_alpha[i]->x, c_alpha[i]->y, c_alpha[i]->z); } fprintf(out,"ENDMDL\n"); fclose(out); } } // calculate gradients e_pot=e_pot1=e_pot2=e_pot3=0.; for (i=0; i<chain_length; i++) gradient[i][0]=gradient[i][1]=gradient[i][2]=0.; e_pot = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, 0.0, energies, true); if (_VERBOSE && !init) { printf("Initial energy: bond=%.5lf angle=%.5f restraints=%.5f xvol=%.5f total=%.5f\n", energies[0], energies[2], energies[1], energies[3], e_pot); } if (!init) init=1; // LINE SEARCH alpha1 = -1.0; alpha2 = 0.0; alpha3 = 1.0; ene1 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha1, energies, false); ene2 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha2, energies, false); ene3 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha3, energies, false); msteps = 0; while (ene2>MIN(ene1,ene3) && msteps<_CA_ITER) { msteps++; alpha1 *= 2.0; alpha3 *= 2.0; ene1 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha1, energies, false); ene3 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha3, energies, false); } msteps = 0; do { if (alpha3-alpha2>alpha2-alpha1) { a0 = 0.5*(alpha2+alpha3); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0-1e-5, energies, false); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0+1e-5, energies, false); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false); if (e0<ene2) { alpha1 = alpha2; alpha2 = a0; ene1 = ene2; ene2 = e0; } else { alpha3 = a0; ene3 = e0; } } else { a0 = 0.5*(alpha1+alpha2); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0-1e-5, energies, false); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0+1e-5, energies, false); e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false); if (e0<ene2) { alpha3 = alpha2; alpha2 = a0; ene3 = ene2; ene2 = e0; } else { alpha1 = a0; ene1 = e0; } } msteps++; } while (alpha3-alpha1>1e-6 && msteps<20); last_alpha = alpha2; e_pot = ene2; for (i=0; i<chain_length; i++) { c_alpha[i]->x=c_alpha[i]->x+(last_alpha+last_alpha*(rnd()-0.5)*eps)*gradient[i][0]; c_alpha[i]->y=c_alpha[i]->y+(last_alpha+last_alpha*(rnd()-0.5)*eps)*gradient[i][1]; c_alpha[i]->z=c_alpha[i]->z+(last_alpha+last_alpha*(rnd()-0.5)*eps)*gradient[i][2]; } e_pot = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, 0.0, energies, false); eps *= 0.75; if (eps<1e-3) eps=0.0; numsteps++; gnorm = 0.0; for (i=0; i<chain_length; i++) { gnorm += gradient[i][0]*gradient[i][0] + gradient[i][1]*gradient[i][1] + gradient[i][2]*gradient[i][2]; } gnorm = sqrt(gnorm/(double)chain_length); if (last_gnorm-gnorm<1e-3) fcnt++; last_gnorm = gnorm; } while ( (fcnt<3) && (gnorm>0.01) && (numsteps<_CA_ITER)); if (_VERBOSE) { for (i=0; i<chain_length; i++) { #ifdef CALC_C_ALPHA if (i>0) { dx=c_alpha[i]->x-c_alpha[i-1]->x; dy=c_alpha[i]->y-c_alpha[i-1]->y; dz=c_alpha[i]->z-c_alpha[i-1]->z; dist=sqrt(dx*dx+dy*dy+dz*dz); if (c_alpha[i]->cispro) { ddist=CA_DIST_CISPRO-dist; if (fabs(ddist)<CA_DIST_CISPRO_TOL) ddist=0.0; } else { ddist=CA_DIST-dist; if (fabs(ddist)<CA_DIST_TOL) ddist=0.0; } ddist2=ddist*ddist; if (fabs(ddist)>=CA_DIST_TOL) printf("WARNING: distance %d = %.3lf A\n", i, dist); } #endif } for (i=0; i<chain_length; i++) { #ifdef CALC_C_ALPHA_ANGLES if (i>0 && i<chain_length-1) { aa=ab=bb=0.0; adiff[0]=c_alpha[i-1]->x-c_alpha[i]->x; bdiff[0]=c_alpha[i+1]->x-c_alpha[i]->x; aa+=adiff[0]*adiff[0]; ab+=adiff[0]*bdiff[0]; bb+=bdiff[0]*bdiff[0]; adiff[1]=c_alpha[i-1]->y-c_alpha[i]->y; bdiff[1]=c_alpha[i+1]->y-c_alpha[i]->y; aa+=adiff[1]*adiff[1]; ab+=adiff[1]*bdiff[1]; bb+=bdiff[1]*bdiff[1]; adiff[2]=c_alpha[i-1]->z-c_alpha[i]->z; bdiff[2]=c_alpha[i+1]->z-c_alpha[i]->z; aa+=adiff[2]*adiff[2]; ab+=adiff[2]*bdiff[2]; bb+=bdiff[2]*bdiff[2]; th=ab/sqrt(aa*bb); if (th<-1.0) th=-1.0; if (th>1.0) th=1.0; th=acos(th); deg_th=RADDEG*th; if (deg_th>150.) theta0=DEGRAD*150.; else if (deg_th<75.) theta0=DEGRAD*75.; else theta0=th; if (fabs(deg_th-RADDEG*theta0)>1.0) printf("WARNING: angle %d = %.3lf degrees\n", i, deg_th); } #endif } } if (_VERBOSE) printf("Optimization done after %d step(s).\nFinal energy: bond=%.5lf angle=%.5f restraints=%.5f xvol=%.5f total=%.5f\n", numsteps, energies[0], energies[2], energies[1], energies[3], e_pot); if (_CA_TRAJECTORY) { out = fopen(tname,"a"); if (out) { fprintf(out,"END\n"); } } for (i=0;i<chain_length+1;i++) { free(init_c_alpha[i]); free(new_c_alpha[i]); free(gradient[i]); } free(new_c_alpha); free(gradient); free(c_alpha); free(init_c_alpha); } void center_chain(mol_type *mol) { real cx, cy, cz; int natom; res_type *res; atom_type *atom; cx = cy = cz = 0.0; natom = 0; res = mol->residua; while (res) { atom = res->atoms; while (atom) { cx += atom->x; cy += atom->y; cz += atom->z; natom++; atom=atom->next; } res = res->next; } cx /= (real)natom; cy /= (real)natom; cz /= (real)natom; if (_VERBOSE) printf("Molecule center: %8.3f %8.3f %8.3f -> 0.000 0.000 0.000\n", cx, cy, cz); res = mol->residua; while (res) { atom = res->atoms; while (atom) { atom->x -= cx; atom->y -= cy; atom->z -= cz; natom++; atom=atom->next; } res = res->next; } } void write_pdb(char *name, mol_type *mol) { FILE *out; res_type *res; atom_type *atom; int anum; out = fopen(name,"w"); if (!out) { if (_VERBOSE) printf("Can't open output file!\n"); return; } fprintf(out,"REMARK 999 REBUILT BY PULCHRA V.%.2f\n", PULCHRA_VERSION); anum=1; res = mol->residua; while (res) { if (res->protein) { if (!_BB_REARRANGE) { atom = res->atoms; while (atom) { if (!(atom->name[0]=='D' && atom->name[1]=='U') && !(atom->name[0]=='S' && atom->name[1]=='C') && !(atom->name[0]=='C' && atom->name[1]=='M') && !(atom->name[0]=='H' && !_REBUILD_H)) fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n", anum++, atom->name, res->name, ' ', res->num, atom->x, atom->y, atom->z); atom=atom->next; } } else { atom = res->atoms; while (atom) { if (!(atom->name[0]=='D' && atom->name[1]=='U') && !(atom->name[0]=='S' && atom->name[1]=='C') && !(atom->name[0]=='C' && atom->name[1]==' ') && !(atom->name[0]=='O' && atom->name[1]==' ') && !(atom->name[0]=='C' && atom->name[1]=='M') && !(atom->name[0]=='H' && !_REBUILD_H)) fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n", anum++, atom->name, res->name, ' ', res->num, atom->x, atom->y, atom->z); atom=atom->next; } atom = res->atoms; while (atom) { if (((atom->name[0]=='C' && atom->name[1]==' ') || (atom->name[0]=='O' && atom->name[1]==' ')) && !(atom->name[0]=='H' && !_REBUILD_H)) fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n", anum++, atom->name, res->name, ' ', res->num, atom->x, atom->y, atom->z); atom=atom->next; } } } res = res->next; } fprintf(out,"TER\nEND\n"); fclose(out); } void write_pdb_sg(char *name, mol_type *mol) { FILE *out; res_type *res; atom_type *atom; int anum; out = fopen(name,"w"); if (!out) { if (_VERBOSE) printf("Can't open output file!\n"); return; } fprintf(out,"REMARK 999 REBUILT BY PULCHRA V.%.2f\n", PULCHRA_VERSION); anum=1; res = mol->residua; while (res) { if (res->protein) { atom = res->atoms; while (atom) { if ((atom->name[0]=='C' && atom->name[1]=='A')) fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n", anum++, atom->name, res->name, ' ', res->num, atom->x, atom->y, atom->z); atom=atom->next; } fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n", anum++, "CM ", res->name, ' ', res->num, res->cmx, res->cmy, res->cmz); } res = res->next; } fprintf(out,"TER\nEND\n"); fclose(out); } real calc_distance(real x1, real y1, real z1, real x2, real y2, real z2) { real dx,dy,dz; real dist2; dx = (x1) - (x2); dy = (y1) - (y2); dz = (z1) - (z2); if (dx || dy || dz ) { dist2 = dx*dx+dy*dy+dz*dz; return (sqrt(dist2)); } else return 0.0; } real calc_r14(real x1, real y1, real z1, real x2, real y2, real z2, real x3, real y3, real z3, real x4, real y4, real z4) { real r, dx, dy, dz; real vx1, vy1, vz1, vx2, vy2, vz2, vx3, vy3, vz3; real hand; dx = x4-x1; dy = y4-y1; dz = z4-z1; r = sqrt(dx*dx+dy*dy+dz*dz); vx1=x2-x1; vy1=y2-y1; vz1=z2-z1; vx2=x3-x2; vy2=y3-y2; vz2=z3-z2; vx3=x4-x3; vy3=y4-y3; vz3=z4-z3; hand = (vy1*vz2-vy2*vz1)*vx3+ (vz1*vx2-vz2*vx1)*vy3+ (vx1*vy2-vx2*vy1)*vz3; if (hand<0) r=-r; return r; } real superimpose2(real **coords1, real **coords2, int npoints, real **tpoints, int ntpoints) { real mat_s[3][3], mat_a[3][3], mat_b[3][3], mat_g[3][3]; real mat_u[3][3], tmp_mat[3][3]; real val, d, alpha, beta, gamma, x, y, z; real cx1, cy1, cz1, cx2, cy2, cz2, tmpx, tmpy, tmpz; int i, j, k, n; cx1=cy1=cz1=cx2=cy2=cz2=0.; for (i=0; i<npoints; i++) { cx1+=coords1[i][0]; cy1+=coords1[i][1]; cz1+=coords1[i][2]; cx2+=coords2[i][0]; cy2+=coords2[i][1]; cz2+=coords2[i][2]; } cx1/=(real)npoints; cy1/=(real)npoints; cz1/=(real)npoints; cx2/=(real)npoints; cy2/=(real)npoints; cz2/=(real)npoints; for (i=0; i<npoints; i++) { coords1[i][0]-=cx1; coords1[i][1]-=cy1; coords1[i][2]-=cz1; coords2[i][0]-=cx2; coords2[i][1]-=cy2; coords2[i][2]-=cz2; } for (i=0; i<ntpoints; i++) { tpoints[i][0]-=cx2; tpoints[i][1]-=cy2; tpoints[i][2]-=cz2; } for (i=0; i<3; i++) for (j=0; j<3; j++) { if (i==j) mat_s[i][j]=mat_a[i][j]=mat_b[i][j]=mat_g[i][j]=1.0; else mat_s[i][j]=mat_a[i][j]=mat_b[i][j]=mat_g[i][j]=0.0; mat_u[i][j]=0.; } for (n=0; n<npoints; n++) { mat_u[0][0]+=coords1[n][0]*coords2[n][0]; mat_u[0][1]+=coords1[n][0]*coords2[n][1]; mat_u[0][2]+=coords1[n][0]*coords2[n][2]; mat_u[1][0]+=coords1[n][1]*coords2[n][0]; mat_u[1][1]+=coords1[n][1]*coords2[n][1]; mat_u[1][2]+=coords1[n][1]*coords2[n][2]; mat_u[2][0]+=coords1[n][2]*coords2[n][0]; mat_u[2][1]+=coords1[n][2]*coords2[n][1]; mat_u[2][2]+=coords1[n][2]*coords2[n][2]; } for (i=0; i<3; i++) for (j=0; j<3; j++) tmp_mat[i][j]=0.; do { d=mat_u[2][1]-mat_u[1][2]; if (d==0) alpha=0; else alpha=atan(d/(mat_u[1][1]+mat_u[2][2])); if (cos(alpha)*(mat_u[1][1]+mat_u[2][2])+sin(alpha)*(mat_u[2][1]-mat_u[1][2])<0.0) alpha+=M_PI; mat_a[1][1]=mat_a[2][2]=cos(alpha); mat_a[2][1]=sin(alpha); mat_a[1][2]=-mat_a[2][1]; for (i=0; i<3; i++) for (j=0; j<3; j++) for (k=0; k<3; k++) tmp_mat[i][j]+=mat_u[i][k]*mat_a[j][k]; for (i=0; i<3; i++) for (j=0; j<3; j++) { mat_u[i][j]=tmp_mat[i][j]; tmp_mat[i][j]=0.; } for (i=0; i<3; i++) for (j=0; j<3; j++) for (k=0; k<3; k++) tmp_mat[i][j]+=mat_a[i][k]*mat_s[k][j]; for (i=0; i<3; i++) for (j=0; j<3; j++) { mat_s[i][j]=tmp_mat[i][j]; tmp_mat[i][j]=0.; } d=mat_u[0][2]-mat_u[2][0]; if (d==0) beta=0; else beta=atan(d/(mat_u[0][0]+mat_u[2][2])); if (cos(beta)*(mat_u[0][0]+mat_u[2][2])+sin(beta)*(mat_u[0][2]-mat_u[2][0])<0.0) beta+=M_PI; mat_b[0][0]=mat_b[2][2]=cos(beta); mat_b[0][2]=sin(beta); mat_b[2][0]=-mat_b[0][2]; for (i=0; i<3; i++) for (j=0; j<3; j++) for (k=0; k<3; k++) tmp_mat[i][j]+=mat_u[i][k]*mat_b[j][k]; for (i=0; i<3; i++) for (j=0; j<3; j++) { mat_u[i][j]=tmp_mat[i][j]; tmp_mat[i][j]=0.; } for (i=0; i<3; i++) for (j=0; j<3; j++) for (k=0; k<3; k++) tmp_mat[i][j]+=mat_b[i][k]*mat_s[k][j]; for (i=0; i<3; i++) for (j=0; j<3; j++) { mat_s[i][j]=tmp_mat[i][j]; tmp_mat[i][j]=0.; } d=mat_u[1][0]-mat_u[0][1]; if (d==0) gamma=0; else gamma=atan(d/(mat_u[0][0]+mat_u[1][1])); if (cos(gamma)*(mat_u[0][0]+mat_u[1][1])+sin(gamma)*(mat_u[1][0]-mat_u[0][1])<0.0) gamma+=M_PI; mat_g[0][0]=mat_g[1][1]=cos(gamma); mat_g[1][0]=sin(gamma); mat_g[0][1]=-mat_g[1][0]; for (i=0; i<3; i++) for (j=0; j<3; j++) for (k=0; k<3; k++) tmp_mat[i][j]+=mat_u[i][k]*mat_g[j][k]; for (i=0; i<3; i++) for (j=0; j<3; j++) { mat_u[i][j]=tmp_mat[i][j]; tmp_mat[i][j]=0.; } for (i=0; i<3; i++) for (j=0; j<3; j++) for (k=0; k<3; k++) tmp_mat[i][j]+=mat_g[i][k]*mat_s[k][j]; for (i=0; i<3; i++) for (j=0; j<3; j++) { mat_s[i][j]=tmp_mat[i][j]; tmp_mat[i][j]=0.; } val=fabs(alpha)+fabs(beta)+fabs(gamma); } while (val>0.001); val=0.; for (i=0; i<npoints; i++) { x=coords2[i][0]; y=coords2[i][1]; z=coords2[i][2]; tmpx=x*mat_s[0][0]+y*mat_s[0][1]+z*mat_s[0][2]; tmpy=x*mat_s[1][0]+y*mat_s[1][1]+z*mat_s[1][2]; tmpz=x*mat_s[2][0]+y*mat_s[2][1]+z*mat_s[2][2]; x=coords1[i][0]-tmpx; y=coords1[i][1]-tmpy; z=coords1[i][2]-tmpz; val+=x*x+y*y+z*z; } for (i=0; i<ntpoints; i++) { x=tpoints[i][0]; y=tpoints[i][1]; z=tpoints[i][2]; tpoints[i][0]=x*mat_s[0][0]+y*mat_s[0][1]+z*mat_s[0][2]; tpoints[i][1]=x*mat_s[1][0]+y*mat_s[1][1]+z*mat_s[1][2]; tpoints[i][2]=x*mat_s[2][0]+y*mat_s[2][1]+z*mat_s[2][2]; } for (i=0; i<npoints; i++) { coords1[i][0]+=cx1; coords1[i][1]+=cy1; coords1[i][2]+=cz1; coords2[i][0]+=cx2; coords2[i][1]+=cy2; coords2[i][2]+=cz2; } for (i=0; i<ntpoints; i++) { tpoints[i][0]+=cx1; tpoints[i][1]+=cy1; tpoints[i][2]+=cz1; } return sqrt(val/(real)npoints); } atom_type *find_atom(res_type *res, char *aname) { atom_type *atom; atom = res->atoms; while (atom) { if (atom->name[0]==aname[0] && atom->name[1]==aname[1] && atom->name[2]==aname[2]) { return atom; break; } atom = atom->next; } return NULL; } void add_replace(res_type *res, char *aname, real x, real y, real z, int flags) { atom_type *atom, *newatom; atom = res->atoms; while (atom) { if (atom->name[0]==aname[0] && atom->name[1]==aname[1] && atom->name[2]==aname[2]) { atom->x = x; atom->y = y; atom->z = z; atom->flag |= flags; break; } atom = atom->next; } if (!atom) { newatom = (atom_type*)calloc(sizeof(atom_type),1); newatom->x = x; newatom->y = y; newatom->z = z; newatom->flag |= flags; newatom->res = res; newatom->name = (char*)calloc(4,1); strcpy(newatom->name,aname); atom = res->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]=='A') break; atom = atom->next; } if (aname[0]=='N' && aname[1]==' ') { newatom->next = res->atoms; res->atoms = newatom; } else { while (atom->next) atom=atom->next; atom->next = newatom; } } } int **RBINS; real **X_COORDS, **C_ALPHA; #ifdef COMPILE_BB void rebuild_backbone(void) { res_type *res, *prevres; atom_type *atom; real **cacoords, **tmpcoords, **tmpstat; real x1, y1, z1; real x2, y2, z2; real x3, y3, z3; real x4, y4, z4; real r13_1, r13_2, r14; real besthit, hit; int bestpos; int i, j, k, l, m, bin13_1, bin13_2, bin14, found, pro; int b13_1, b13_2, b14; real rmsd, total, maxrms; FILE *debug, *out; if (_VERBOSE) printf("Rebuilding backbone...\n"); RBINS = (int**)calloc(sizeof(int*)*(chain_length+1),1); for (i=0;i<chain_length+1;i++) RBINS[i] = (int*)calloc(sizeof(int)*3,1); X_COORDS = (real**)calloc(sizeof(real*)*(chain_length+10),1); for (i=0;i<chain_length+10;i++) X_COORDS[i] = (real*)calloc(sizeof(real)*3,1); i = 5; res = chain->residua; while (res) { atom = res->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]=='A') { X_COORDS[i][0] = atom->x; X_COORDS[i][1] = atom->y; X_COORDS[i][2] = atom->z; i++; } atom = atom->next; } res = res->next; } cacoords = (real**)calloc(sizeof(real*)*(8),1); tmpcoords = (real**)calloc(sizeof(real*)*(8),1); tmpstat = (real**)calloc(sizeof(real*)*(8),1); for (i=0;i<8;i++) { cacoords[i] = (real*)calloc(sizeof(real)*3,1);; tmpcoords[i] = (real*)calloc(sizeof(real)*3,1);; tmpstat[i] = (real*)calloc(sizeof(real)*3,1);; } C_ALPHA = &X_COORDS[5]; // rebuild ends... for (i=0,j=0;i<5;i++,j++) for (k=0;k<3;k++) tmpcoords[j][k] = C_ALPHA[i][k]; for (i=2,j=0;i<5;i++,j++) for (k=0;k<3;k++) cacoords[j][k] = C_ALPHA[i][k]; for (i=0,j=0;i<3;i++,j++) for (k=0;k<3;k++) tmpstat[j][k] = C_ALPHA[i][k]; superimpose2(tmpstat,cacoords,3,tmpcoords,5); for (i=-2,j=0;i<0;i++,j++) for (k=0;k<3;k++) C_ALPHA[i][k] = tmpcoords[j][k]; for (i=chain_length-5,j=0;i<chain_length;i++,j++) for (k=0;k<3;k++) tmpcoords[j][k] = C_ALPHA[i][k]; for (i=chain_length-5,j=0;i<chain_length-2;i++,j++) for (k=0;k<3;k++) cacoords[j][k] = C_ALPHA[i][k]; for (i=chain_length-3,j=0;i<chain_length;i++,j++) for (k=0;k<3;k++) tmpstat[j][k] = C_ALPHA[i][k]; superimpose2(tmpstat,cacoords,3,tmpcoords,5); for (i=chain_length-3,j=0;i<chain_length;i++,j++) for (k=0;k<3;k++) C_ALPHA[i+3][k] = tmpcoords[j+3][k]; prevres = NULL; res = chain->residua; total = maxrms = 0.0; for (i=0;i<chain_length+1;i++) { x1 = C_ALPHA[i-2][0]; y1 = C_ALPHA[i-2][1]; z1 = C_ALPHA[i-2][2]; x2 = C_ALPHA[i-1][0]; y2 = C_ALPHA[i-1][1]; z2 = C_ALPHA[i-1][2]; x3 = C_ALPHA[i][0]; y3 = C_ALPHA[i][1]; z3 = C_ALPHA[i][2]; x4 = C_ALPHA[i+1][0]; y4 = C_ALPHA[i+1][1]; z4 = C_ALPHA[i+1][2]; r13_1 = calc_distance(x1, y1, z1, x3, y3, z3); r13_2 = calc_distance(x2, y2, z2, x4, y4, z4); r14 = calc_r14(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4); bin13_1 = (int)((r13_1-4.6)/0.3); bin13_2 = (int)((r13_2-4.6)/0.3); bin14 = (int)((r14+11.)/0.3); if (bin13_1<0) bin13_1=0; if (bin13_2<0) bin13_2=0; if (bin14<0) bin14=0; if (bin13_1>9) bin13_1=9; if (bin13_2>9) bin13_2=9; if (bin14>73) bin14=73; RBINS[i][0] = bin13_1; RBINS[i][1] = bin13_2; RBINS[i][2] = bin14; cacoords[0][0] = x1; cacoords[0][1] = y1; cacoords[0][2] = z1; cacoords[1][0] = x2; cacoords[1][1] = y2; cacoords[1][2] = z2; cacoords[2][0] = x3; cacoords[2][1] = y3; cacoords[2][2] = z3; cacoords[3][0] = x4; cacoords[3][1] = y4; cacoords[3][2] = z4; pro = 0; if (prevres && !strncmp(prevres->name,"PRO",3)) { j=0; besthit=1000.; bestpos=0; do { hit = abs(nco_stat_pro[j].bins[0]-bin13_1)+abs(nco_stat_pro[j].bins[1]-bin13_2)+0.2*abs(nco_stat_pro[j].bins[2]-bin14); if (hit<besthit) { besthit=hit; bestpos=j; } j++; } while (nco_stat_pro[j].bins[0]>=0 && hit>1e-3); for (j=0;j<4;j++) { for (k=0;k<3;k++) { tmpstat[j][k] = nco_stat_pro[bestpos].data[j][k]; } } for (j=0;j<8;j++) { for (k=0;k<3;k++) { tmpcoords[j][k] = nco_stat_pro[bestpos].data[j][k]; } } } else { j=0; besthit=1000.; bestpos=0; do { hit = abs(nco_stat[j].bins[0]-bin13_1)+abs(nco_stat[j].bins[1]-bin13_2)+0.2*abs(nco_stat[j].bins[2]-bin14); if (hit<besthit) { besthit=hit; bestpos=j; } j++; } while (nco_stat[j].bins[0]>=0 && hit>1e-3); for (j=0;j<4;j++) { for (k=0;k<3;k++) { tmpstat[j][k] = nco_stat[bestpos].data[j][k]; } } for (j=0;j<8;j++) { for (k=0;k<3;k++) { tmpcoords[j][k] = nco_stat[bestpos].data[j][k]; } } } rmsd=superimpose2(cacoords, tmpstat, 4, tmpcoords, 8); total += rmsd; if (rmsd>maxrms) maxrms=rmsd; // add-or-replace if (prevres) { add_replace(prevres, "C ", tmpcoords[4][0], tmpcoords[4][1], tmpcoords[4][2], FLAG_BACKBONE); add_replace(prevres, "O ", tmpcoords[5][0], tmpcoords[5][1], tmpcoords[5][2], FLAG_BACKBONE); } if (res) { add_replace(res, "N ", tmpcoords[6][0], tmpcoords[6][1], tmpcoords[6][2], FLAG_BACKBONE); } prevres = res; if (res) res = res->next; } if (_VERBOSE) printf("Backbone rebuilding deviation: average = %.3f, max = %.3f\n", total/(real)chain_length, maxrms); } #endif #ifdef COMPILE_ROT typedef struct _rot_struct { int r13_1, r13_2, r14; int nc; real ***coords; struct _rot_struct *next; } rot_struct; rot_struct *rotamers[20]; /* this is obsolete in a standalone version of PULCHRA */ void read_rotamers(void) { FILE *inp; char buf[1000]; char dum[100]; int aa, i, j, k, l, n; rot_struct *new_rot, *last_rot; real x, y, z; if (_VERBOSE) printf("Reading rotamer library...\n"); inp = fopen("NEWROT","r"); last_rot=NULL; while (!feof(inp)) { if (fgets(buf,1000,inp)==buf) { if (buf[0]=='A') { sscanf(buf,"%s %d", dum, &aa); if (last_rot) last_rot->next = NULL; last_rot = NULL; if (fgets(buf,1000,inp)!=buf) break; } // printf("aa: %d\n", aa); if (aa==20) break; sscanf(buf,"%d %d %d %s %d", &i, &j, &k, dum, &l); new_rot = (rot_struct*)calloc(sizeof(rot_struct),1); // printf("%d %d %d nc: %d\n", i, j, k, l); new_rot->r13_1 = i; new_rot->r13_2 = j; new_rot->r14 = k; new_rot->nc = l; new_rot->next = NULL; new_rot->coords = (real***)calloc(sizeof(real**)*l,1); for (i=0;i<l;i++) { new_rot->coords[i]=(real**)calloc(sizeof(real*)*(nheavy[aa]+1),1); for (j=0;j<(nheavy[aa]+1);j++) { new_rot->coords[i][j]=(real*)calloc(sizeof(real)*3,1); } } for (i=0;i<l;i++) { fgets(buf,1000,inp); for (j=0;j<(nheavy[aa]+1);j++) { fgets(buf,1000,inp); sscanf(buf,"%lf%lf%lf",&x, &y, &z); new_rot->coords[i][j][0]=x; new_rot->coords[i][j][1]=y; new_rot->coords[i][j][2]=z; } if (last_rot) { last_rot->next = new_rot; } else { rotamers[aa] = new_rot; } last_rot = new_rot; } } } fclose(inp); } void cross(real *v1, real *v2, real *v3) { v3[0] = v1[1]*v2[2]-v1[2]*v2[1]; v3[1] = v1[2]*v2[0]-v1[0]*v2[2]; v3[2] = v1[0]*v2[1]-v1[1]*v2[0]; } void norm(real *v) { real d; d = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]); v[0] /= d; v[1] /= d; v[2] /= d; } int check_xvol(res_type *res) { res_type *res2; atom_type *atom1, *atom2; real dx, dy, dz, dd; res2 = chain->residua; while (res2) { atom2 = res2->atoms; if (res!=res2) { while (atom2) { atom1 = res->atoms; while (atom1) { if (atom1->flag & FLAG_SIDECHAIN) { dx = atom1->x-atom2->x; dy = atom1->y-atom2->y; dz = atom1->z-atom2->z; dd = dx*dx+dy*dy+dz*dz; if (dd<(1.7*1.7)) { return 1; } } atom1=atom1->next; } atom2=atom2->next; } } res2=res2->next; } return 0; } real ***SORTED_ROTAMERS; void rebuild_sidechains(void) { FILE *out; res_type *res, *prevres, *testres; atom_type *atom, *atom1, *atom2; real **cacoords, **tmpcoords, **tmpstat; real x1, y1, z1; real x2, y2, z2; real x3, y3, z3; real x4, y4, z4; real x5, y5, z5; real r14, r13_1, r13_2; real dx, dy, dz, dd; real hit, besthit; int exvol, bestpos; int i, j, k, l, m, bin13_1, bin13_2, bin14; real rmsd, total; real v1[3], v2a[3], v2b[3], v2[3], v3[3]; int nsc, nca; real cax, cay, caz; real **lsys, **vv, **sc; char scn[12][4]; rot_struct *rot; int ok, last_a, last_b, last_c, last_d, jpos; int jx, jy, jz, jxi, jyi, jzi, b13_1, b13_2, b14, jm; int crot, bestrot, minexvol, totexvol, rtried, pos, cpos; real cmx, cmy, cmz, ddx, ddy, ddz, ddd, bestdd; real sort_rot[100][2]; if (_VERBOSE) printf("Rebuilding side chains...\n"); lsys = (real**)calloc(sizeof(real*)*3,1); vv = (real**)calloc(sizeof(real*)*3,1); sc = (real**)calloc(sizeof(real*)*12,1); for (i=0;i<12;i++) sc[i] = (real*)calloc(sizeof(real)*3,1); for (i=0;i<3;i++) { lsys[i] = (real*)calloc(sizeof(real)*3,1); vv[i] = (real*)calloc(sizeof(real)*3,1); } SORTED_ROTAMERS = (real***)calloc(sizeof(real**)*(chain_length+1),1); for (i=0;i<chain_length+1;i++) { SORTED_ROTAMERS[i] = (real**)calloc(sizeof(real*)*10,1); for (j=0;j<10;j++) { SORTED_ROTAMERS[i][j] = (real*)calloc(sizeof(real)*2,1); } } prevres = NULL; res = chain->residua; totexvol = 0; for (i=0;i<chain_length;i++) { if (!strncmp(res->name,"GLY",3) || !res->protein) { if (res->next) res = res->next; continue; } x1 = C_ALPHA[i-2][0]; y1 = C_ALPHA[i-2][1]; z1 = C_ALPHA[i-2][2]; x2 = C_ALPHA[i-1][0]; y2 = C_ALPHA[i-1][1]; z2 = C_ALPHA[i-1][2]; x3 = C_ALPHA[i][0]; y3 = C_ALPHA[i][1]; z3 = C_ALPHA[i][2]; x4 = C_ALPHA[i+1][0]; y4 = C_ALPHA[i+1][1]; z4 = C_ALPHA[i+1][2]; bin13_1 = RBINS[i][0]; bin13_2 = RBINS[i][1]; bin14 = RBINS[i][2]; v1[0] = x4-x2; v1[1] = y4-y2; v1[2] = z4-z2; v2a[0] = x4-x3; v2a[1] = y4-y3; v2a[2] = z4-z3; v2b[0] = x3-x2; v2b[1] = y3-y2; v2b[2] = z3-z2; cross(v2a, v2b, v2); cross(v1, v2, v3); norm(v1); norm(v2); norm(v3); // gather 10 closest rotamer conformations... for (j=0;j<10;j++) SORTED_ROTAMERS[i][j][0] = 500.; j = 0; besthit = 1000.; bestpos = 0; do { if (rot_stat_idx[j][0]==res->type) { hit = abs(rot_stat_idx[j][1]-bin13_1)+abs(rot_stat_idx[j][2]-bin13_2)+0.2*abs(rot_stat_idx[j][3]-bin14); if (hit<SORTED_ROTAMERS[i][9][0]) { k = 9; while (k>=0 && hit<SORTED_ROTAMERS[i][k][0]) { k--; } k++; // k = hit for (l=9;l>k;l--) { SORTED_ROTAMERS[i][l][0]=SORTED_ROTAMERS[i][l-1][0]; SORTED_ROTAMERS[i][l][1]=SORTED_ROTAMERS[i][l-1][1]; } SORTED_ROTAMERS[i][k][0]=hit; SORTED_ROTAMERS[i][k][1]=j; } } j++; } while (rot_stat_idx[j][0]>=0); besthit = SORTED_ROTAMERS[i][0][0]; bestpos = SORTED_ROTAMERS[i][0][1]; // new rebuild... pos = rot_stat_idx[bestpos][5]; nsc = nheavy[res->type]+1; if (_PDB_SG) { // more than one rotamer - check SC bestdd = 100.; crot = 0; for (l=0;l<2;l++) { // check two closest conformations cpos = SORTED_ROTAMERS[i][l][1]; for (m=0;m<rot_stat_idx[cpos][4];m++) { for (j=0;j<3;j++) { vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j]; for (k=0;k<3;k++) { if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.; } } pos = rot_stat_idx[cpos][5]+nsc*m; for (j=0;j<nsc;j++) { for (k=0;k<3;k++) { sc[j][k] = rot_stat_coords[pos+j][k]; } } superimpose2(vv,lsys,3,sc,nsc); for (j=0;j<nsc;j++) { sc[j][0] += x3; sc[j][1] += y3; sc[j][2] += z3; } cmx = 0.; cmy = 0.; cmz = 0.; for (j=0;j<nsc;j++) { cmx += sc[j][0]; cmy += sc[j][1]; cmz += sc[j][2]; } cmx /= (real) nsc; cmy /= (real) nsc; cmz /= (real) nsc; ddx = res->cmx-cmx; ddy = res->cmy-cmy; ddz = res->cmz-cmz; ddx *= ddx; ddy *= ddy; ddz *= ddz; ddd = ddx+ddy+ddz; if (ddd<bestdd) { bestdd = ddd; crot = pos; // closest rotamer position } } } pos = crot; } // PDB_SG for (j=0;j<3;j++) { vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j]; for (k=0;k<3;k++) { if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.; } } for (j=0;j<nsc;j++) { for (k=0;k<3;k++) { sc[j][k] = rot_stat_coords[pos+j][k]; } } superimpose2(vv,lsys,3,sc,nsc); for (j=0;j<nsc;j++) { sc[j][0] += x3; sc[j][1] += y3; sc[j][2] += z3; } for (j=1;j<nsc;j++) { add_replace(res, heavy_atoms[10*res->type+j-1], sc[j][0], sc[j][1], sc[j][2], FLAG_SIDECHAIN); } if (res->next) res = res->next; } // i++, next res for (i=0;i<12;i++) free(sc[i]); for (i=0;i<3;i++) { free(lsys[i]); free(vv[i]); } free(sc); free(lsys); free(vv); } typedef struct _atom_list { atom_type *atom; struct _atom_list *next; } atom_list; int get_conflicts(res_type *res, atom_list ****grid, int xgrid, int ygrid, int zgrid) { atom_list *llist; atom_type *atom, *atom2; int i, j, k, x, y, z; int ii, jj, kk, con, iter, maxcon, merged; real dx, dy, dz, dd; con = 0; atom = res->atoms; while (atom) { i = atom->gx; j = atom->gy; k = atom->gz; for (ii=i-2;ii<=i+2;ii++) for (jj=j-2;jj<=j+2;jj++) for (kk=k-2;kk<=k+2;kk++) { if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<zgrid) { llist = grid[ii][jj][kk]; while (llist) { atom2 = llist->atom; if (atom && atom2 && res && atom2->res) { merged=0; if (res==atom2->res) { // self-xvol if (atom->flag & FLAG_SIDECHAIN && atom2->flag & FLAG_SIDECHAIN) merged=1; if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1; if (atom->name[0]=='C' && atom->name[1]=='A' && atom2->name[0]=='C' && atom2->name[1]=='B') merged=1; if (atom->name[0]=='C' && atom->name[1]=='B' && atom2->name[0]=='C' && atom2->name[1]=='A') merged=1; if (res->name[0]=='P') { if (atom->name[0]=='C' && atom->name[1]=='D' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1; if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]=='D') merged=1; } if (!merged) { // printf("merged: %s[%d] %s-%s %d %d\n", res->name,res->num,atom->name,atom2->name,atom->flag,atom2->flag); } } else if (res->next==atom2->res || res==atom2->res->next) { if (atom->name[0]=='C' && atom->name[1]==' ' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1; if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]==' ') merged=1; } if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1; // for now if (atom->flag & FLAG_SCM || atom2->flag & FLAG_SCM) merged=1; // for now if (!merged) { dx = atom->x-atom2->x; dx*=dx; dy = atom->y-atom2->y; dy*=dy; dz = atom->z-atom2->z; dz*=dz; dd = dx+dy+dz; if (dd<_SG_XVOL_DIST*_SG_XVOL_DIST) { con++; } } } llist = llist->next; } } } atom = atom->next; } return con; } int display_conflicts(res_type *res, atom_list ****grid, int xgrid, int ygrid, int zgrid) { atom_list *llist; atom_type *atom, *atom2; int i, j, k, x, y, z; int ii, jj, kk, con, iter, maxcon, merged; real dx, dy, dz, dd; con = 0; atom = res->atoms; while (atom) { i = atom->gx; j = atom->gy; k = atom->gz; for (ii=i-2;ii<=i+2;ii++) for (jj=j-2;jj<=j+2;jj++) for (kk=k-2;kk<=k+2;kk++) { if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<zgrid) { llist = grid[ii][jj][kk]; while (llist) { atom2 = llist->atom; if (atom && atom2 && res && atom2->res) { merged=0; if (res==atom2->res) { // self-xvol if (atom->flag & FLAG_SIDECHAIN && atom2->flag & FLAG_SIDECHAIN) merged=1; if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1; if (atom->name[0]=='C' && atom->name[1]=='A' && atom2->name[0]=='C' && atom2->name[1]=='B') merged=1; if (atom->name[0]=='C' && atom->name[1]=='B' && atom2->name[0]=='C' && atom2->name[1]=='A') merged=1; if (res->name[0]=='P') { if (atom->name[0]=='C' && atom->name[1]=='D' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1; if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]=='D') merged=1; } if (!merged) { // printf("merged: %s[%d] %s-%s %d %d\n", res->name,res->num,atom->name,atom2->name,atom->flag,atom2->flag); } } else if (res->next==atom2->res || res==atom2->res->next) { if (atom->name[0]=='C' && atom->name[1]==' ' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1; if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]==' ') merged=1; } if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1; // for now if (atom->flag & FLAG_SCM || atom2->flag & FLAG_SCM) merged=1; // for now if (!merged) { dx = atom->x-atom2->x; dx*=dx; dy = atom->y-atom2->y; dy*=dy; dz = atom->z-atom2->z; dz*=dz; dd = dx+dy+dz; if (dd<1.6*1.6) { printf("STERIC CONFLICT: %s[%d]%s-%s[%d]%s\n", atom->res->name,atom->res->num,atom->name,atom2->res->name,atom2->res->num,atom2->name); con++; } } } llist = llist->next; } } } atom = atom->next; } return con; } void allocate_grid(atom_list *****grid_, int *xgrid_, int *ygrid_, int *zgrid_) { static int xgrid, ygrid, zgrid; static atom_list ****grid = NULL; atom_list *llist, *alist; real min[3], max[3]; res_type *res, *worst; atom_type *atom, *atom2; int i, j, x, y, z; if (!grid && chain->residua && chain->residua->atoms) { res = chain->residua; min[0]=max[0]=res->atoms->x; min[1]=max[1]=res->atoms->y; min[2]=max[2]=res->atoms->z; while (res) { atom = res->atoms; while (atom) { if (atom->x<min[0]) min[0]=atom->x; if (atom->y<min[1]) min[1]=atom->y; if (atom->z<min[2]) min[2]=atom->z; if (atom->x>max[0]) max[0]=atom->x; if (atom->y>max[1]) max[1]=atom->y; if (atom->z>max[2]) max[2]=atom->z; atom = atom->next; } res = res->next; } xgrid = (max[0]-min[0])/GRID_RES; ygrid = (max[1]-min[1])/GRID_RES; zgrid = (max[2]-min[2])/GRID_RES; if (_VERBOSE) printf("Allocating grid (%d %d %d)...\n", xgrid, ygrid, zgrid); grid = (atom_list****)calloc(sizeof(atom_list***)*(xgrid+1),1); for (i=0;i<xgrid+1;i++) { grid[i] = (atom_list***)calloc(sizeof(atom_list**)*(ygrid+1),1); for (j=0;j<ygrid+1;j++) { grid[i][j] = (atom_list**)calloc(sizeof(atom_list*)*(zgrid+1),1); } } res = chain->residua; while (res) { atom = res->atoms; while (atom) { x = xgrid*(atom->x-min[0])/(max[0]-min[0]); y = ygrid*(atom->y-min[1])/(max[1]-min[1]); z = zgrid*(atom->z-min[2])/(max[2]-min[2]); alist = (atom_list*)calloc(sizeof(atom_list),1); alist->atom = atom; atom->gx = x; atom->gy = y; atom->gz = z; if (grid[x][y][z]!=NULL) { llist = grid[x][y][z]; while (llist->next) llist=llist->next; llist->next = alist; } else { grid[x][y][z]=alist; } atom = atom->next; } res = res->next; } } else { if (_VERBOSE) printf("Grid already allocated (%d %d %d)\n", xgrid, ygrid, zgrid); } *grid_ = grid; *xgrid_ = xgrid; *ygrid_ = ygrid; *zgrid_ = zgrid; } void optimize_exvol(void) { real min[3], max[3]; res_type *res, *worst; atom_type *atom, *atom2; int xgrid, ygrid, zgrid; atom_list ****grid, *llist, *alist; int i, j, k, l, m, x, y, z; int ii, jj, kk, con, iter, maxcon, totcon; int cpos, bestpos, pos, con0; real v1[3], v2a[3], v2b[3], v2[3], v3[3]; int nsc, nca; real cax, cay, caz; real **lsys, **vv, **sc; real x1, y1, z1; real x2, y2, z2; real x3, y3, z3; real x4, y4, z4; min[0]=1e5; min[1]=1e5; min[2]=1e5; max[0]=-1e5; max[1]=-1e5; max[2]=-1e5; lsys = (real**)calloc(sizeof(real*)*3,1); vv = (real**)calloc(sizeof(real*)*3,1); sc = (real**)calloc(sizeof(real*)*12,1); for (i=0;i<12;i++) sc[i] = (real*)calloc(sizeof(real)*3,1); for (i=0;i<3;i++) { lsys[i] = (real*)calloc(sizeof(real)*3,1); vv[i] = (real*)calloc(sizeof(real)*3,1); } allocate_grid(&grid, &xgrid, &ygrid, &zgrid); if (_VERBOSE) printf("Finding excluded volume conflicts...\n"); iter = 0; do { //printf("ITER: %d\n", iter); maxcon = 0; totcon=0; res = chain->residua; while (res) { if (res->protein) { con = get_conflicts(res, grid, xgrid, ygrid, zgrid); if (con>0) { totcon+=con; if (con>maxcon) { maxcon = con; worst = res; } } } res = res->next; } if (_VERBOSE && iter==0) { printf("Total number of conflicts: %d\n", totcon); } if (totcon==0) break; if (_VERBOSE && iter==0) { printf("Maximum number of conflicts: %s[%d] : %d\n", worst->name, worst->num, maxcon); } totcon=0; if (maxcon>0) { // try to fix... res = chain->residua; for (i=0;i<chain_length;i++) { if (!strncmp(res->name,"GLY",3) || !res->protein) { if (res->next) res = res->next; continue; } nsc = nheavy[res->type]+1; x1 = C_ALPHA[i-2][0]; y1 = C_ALPHA[i-2][1]; z1 = C_ALPHA[i-2][2]; x2 = C_ALPHA[i-1][0]; y2 = C_ALPHA[i-1][1]; z2 = C_ALPHA[i-1][2]; x3 = C_ALPHA[i][0]; y3 = C_ALPHA[i][1]; z3 = C_ALPHA[i][2]; x4 = C_ALPHA[i+1][0]; y4 = C_ALPHA[i+1][1]; z4 = C_ALPHA[i+1][2]; v1[0] = x4-x2; v1[1] = y4-y2; v1[2] = z4-z2; v2a[0] = x4-x3; v2a[1] = y4-y3; v2a[2] = z4-z3; v2b[0] = x3-x2; v2b[1] = y3-y2; v2b[2] = z3-z2; cross(v2a, v2b, v2); cross(v1, v2, v3); norm(v1); norm(v2); norm(v3); con = get_conflicts(res, grid, xgrid, ygrid, zgrid); if (con>0) { bestpos=0; con0 = 100; for (l=0;l<10;l++) { // check two closest conformations cpos = SORTED_ROTAMERS[i][l][1]; for (m=0;m<rot_stat_idx[cpos][4];m++) { for (j=0;j<3;j++) { vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j]; for (k=0;k<3;k++) { if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.; } } pos = rot_stat_idx[cpos][5]+nsc*m; for (j=0;j<nsc;j++) { for (k=0;k<3;k++) { sc[j][k] = rot_stat_coords[pos+j][k]; } } superimpose2(vv,lsys,3,sc,nsc); for (j=0;j<nsc;j++) { sc[j][0] += x3; sc[j][1] += y3; sc[j][2] += z3; } for (j=1;j<nsc;j++) { add_replace(res, heavy_atoms[10*res->type+j-1], sc[j][0], sc[j][1], sc[j][2], FLAG_SIDECHAIN); } con = get_conflicts(res, grid, xgrid, ygrid, zgrid); //printf("test: %d\n", con); if (con<con0) { con0 = con; bestpos = pos; } if (con==0) break; } if (con==0) break; } totcon += con0; for (j=0;j<3;j++) { vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j]; for (k=0;k<3;k++) { if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.; } } pos = bestpos; for (j=0;j<nsc;j++) { for (k=0;k<3;k++) { sc[j][k] = rot_stat_coords[pos+j][k]; } } superimpose2(vv,lsys,3,sc,nsc); for (j=0;j<nsc;j++) { sc[j][0] += x3; sc[j][1] += y3; sc[j][2] += z3; } for (j=1;j<nsc;j++) { add_replace(res, heavy_atoms[10*res->type+j-1], sc[j][0], sc[j][1], sc[j][2], FLAG_SIDECHAIN); } } res=res->next; } // i } iter++; } while (iter<_XVOL_ITER); if (_VERBOSE) { if (totcon>0) printf("WARNING: %d steric conflict(s) are still there.\n", totcon); else printf("All steric conflicts removed.\n"); } for (i=0;i<12;i++) free(sc[i]); for (i=0;i<3;i++) { free(lsys[i]); free(vv[i]); } free(sc); free(lsys); free(vv); } void vcross(real ax,real ay,real az,real bx,real by,real bz,real *cx,real *cy,real *cz) { *cx = ay * bz - by * az; *cy = az * bx - bz * ax; *cz = ax * by - bx * ay; } real vdot(real ax,real ay,real az,real bx,real by,real bz) { return ax*bx+ay*by+az*bz; } real calc_torsion(atom_type *a1, atom_type *a2, atom_type *a3, atom_type *a4) { real v12x, v12y, v12z; real v43x, v43y, v43z; real zx, zy, zz; real px, py, pz; real xx, xy, xz; real yx, yy, yz; real u, v, angle; v12x = a1->x-a2->x; v12y = a1->y-a2->y; v12z = a1->z-a2->z; v43x = a4->x-a3->x; v43y = a4->y-a3->y; v43z = a4->z-a3->z; zx = a2->x-a3->x; zy = a2->y-a3->y; zz = a2->z-a3->z; vcross(zx,zy,zz,v12x,v12y,v12z,&px,&py,&pz); vcross(zx,zy,zz,v43x,v43y,v43z,&xx,&xy,&xz); vcross(zx,zy,zz,xx,xy,xz,&yx,&yy,&yz); u = vdot(xx,xy,xz,xx,xy,xz); v = vdot(yx,yy,yz,yx,yy,yz); angle = 360.; if (u<0. || v<0.) return angle; u = vdot(px,py,pz,xx,xy,xz) / sqrt(u); v = vdot(px,py,pz,yx,yy,yz) / sqrt(v); if (u != 0.0 || v != 0.0) angle = atan2(v, u) * RADDEG; return angle; } // Ca-N-C-Cb angle should be close to 34 deg void chirality_check(void) { int i; atom_type *a_ca, *a_n, *a_c, *a_cb; atom_type *atom; res_type *res; real angle; real nx, ny, nz; real px, py, pz; real qx, qy, qz; real rx, ry, rz; real xx, xy, xz; real yx, yy, yz; real dd, costheta, sintheta; if (_VERBOSE) printf("Checking chirality...\n"); res = chain->residua; while (res) { a_ca = a_n = a_c = a_cb = NULL; a_ca = find_atom(res,"CA "); a_n = find_atom(res,"N "); a_c = find_atom(res,"C "); a_cb = find_atom(res,"CB "); if (a_ca && a_n && a_c && a_cb) { angle = calc_torsion(a_ca, a_n, a_c, a_cb); if (angle<0.) { if (_VERBOSE) printf("WARNING: D-aa detected at %s %3d : %5.2f", res->name, res->num, angle); xx = a_ca->x-a_n->x; xy = a_ca->y-a_n->y; xz = a_ca->z-a_n->z; yx = a_c->x-a_ca->x; yy = a_c->y-a_ca->y; yz = a_c->z-a_ca->z; vcross(xx,xy,xz,yx,yy,yz,&nx,&ny,&nz); dd = sqrt(nx*nx+ny*ny+nz*nz); nx /= dd; ny /= dd; nz /= dd; // nx, ny, nz = reflection plane normal rx = xx-yx; ry = xy-yy; rz = xz-yz; dd = sqrt(rx*rx+ry*ry+rz*rz); rx /= dd; ry /= dd; rz /= dd; costheta = -1.; sintheta = 0.; atom = res->atoms; while (atom) { if (atom->flag & FLAG_SIDECHAIN) { px = atom->x-a_ca->x; py = atom->y-a_ca->y; pz = atom->z-a_ca->z; qx = qy = qz = 0.; qx += (costheta + (1 - costheta) * rx * rx) * px; qx += ((1 - costheta) * rx * ry - rz * sintheta) * py; qx += ((1 - costheta) * rx * rz + ry * sintheta) * pz; qy += ((1 - costheta) * rx * ry + rz * sintheta) * px; qy += (costheta + (1 - costheta) * ry * ry) * py; qy += ((1 - costheta) * ry * rz - rx * sintheta) * pz; qz += ((1 - costheta) * rx * rz - ry * sintheta) * px; qz += ((1 - costheta) * ry * rz + rx * sintheta) * py; qz += (costheta + (1 - costheta) * rz * rz) * pz; qx += a_ca->x; qy += a_ca->y; qz += a_ca->z; atom->x = qx; atom->y = qy; atom->z = qz; } atom = atom->next; } angle = calc_torsion(a_ca, a_n, a_c, a_cb); if (_VERBOSE) printf(", fixed : %5.2f\n", angle); } } res = res->next; } } #endif real hb_energy(res_type *res, atom_list ****grid, int xgrid, int ygrid, int zgrid) { atom_type *atom, *c_atom1, *o_atom1, *n_atom1, *c_atom2, *o_atom2, *n_atom2, *tmp_atom; atom_type h_atom; int i, j, k, ii, jj, kk; atom_list *llist, *alist; real dx, dy, dz, dist, min_dist1, min_dist2; real hx1, hy1, hz1, dd; real dno, dnc, dho, dhc; real ene, Q; ene = 1e3; if (!res || !res->prev) return ene; Q = -27888.0; // DSSP h-bond energy constant c_atom1 = o_atom1 = n_atom1 = NULL; atom = res->prev->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]==' ') c_atom1 = atom; if (atom->name[0]=='O' && atom->name[1]==' ') o_atom1 = atom; atom = atom->next; } atom = res->atoms; while (atom) { if (atom->name[0]=='N' && atom->name[1]==' ') { n_atom1 = atom; break; } atom = atom->next; } // first bond min_dist2 = 1e10; o_atom2 = c_atom2 = NULL; if (n_atom1) { i = n_atom1->gx; j = n_atom1->gy; k = n_atom1->gz; for (ii=i-1;ii<=i+1;ii++) { for (jj=j-1;jj<=j+1;jj++) { for (kk=k-1;kk<=k+1;kk++) { if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<=zgrid) { llist = grid[ii][jj][kk]; while (llist) { if (llist->atom->name[0]=='O' && llist->atom->name[1]==' ' && abs(llist->atom->res->locnum-n_atom1->res->locnum)>2) { tmp_atom = llist->atom; dx = n_atom1->x-tmp_atom->x; dy = n_atom1->y-tmp_atom->y; dz = n_atom1->z-tmp_atom->z; dist = dx*dx+dy*dy+dz*dz; if (dist<min_dist2 && dist<25.0) { o_atom2=tmp_atom; min_dist2 = dist; } } llist = llist->next; } } } } } } if (o_atom2) { atom = o_atom2->res->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]==' ') { c_atom2 = atom; break; } atom = atom->next; } if (c_atom2) { hx1 = o_atom1->x-c_atom1->x; hy1 = o_atom1->y-c_atom1->y; hz1 = o_atom1->z-c_atom1->z; dd = -1.081f/sqrt(hx1*hx1+hy1*hy1+hz1*hz1); hx1 *= dd; hy1 *= dd; hz1 *= dd; hx1 += n_atom1->x; hy1 += n_atom1->y; hz1 += n_atom1->z; add_replace(n_atom1->res, "H ", hx1, hy1, hz1, FLAG_BACKBONE); // dno dx = n_atom1->x-o_atom2->x; dy = n_atom1->y-o_atom2->y; dz = n_atom1->z-o_atom2->z; dno = sqrt(dx*dx+dy*dy+dz*dz); // dnc dx = n_atom1->x-c_atom2->x; dy = n_atom1->y-c_atom2->y; dz = n_atom1->z-c_atom2->z; dnc = sqrt(dx*dx+dy*dy+dz*dz); // dho dx = hx1-o_atom2->x; dy = hy1-o_atom2->y; dz = hz1-o_atom2->z; dho = sqrt(dx*dx+dy*dy+dz*dz); // dhc dx = hx1-c_atom2->x; dy = hy1-c_atom2->y; dz = hz1-c_atom2->z; dhc = sqrt(dx*dx+dy*dy+dz*dz); if (dho<0.01F || dhc<0.01F || dnc<0.01F || dno<0.01F) { ene = -10.0; } else { ene = 0.001*(Q/dho - Q/dhc + Q/dnc - Q/dno); } } } /****** // second bond min_dist2 = 1e10; n_atom2 = NULL; if (n_atom1) { i = o_atom1->gx; j = o_atom1->gy; k = o_atom1->gz; for (ii=i-1;ii<=i+1;ii++) { for (jj=j-1;jj<=j+1;jj++) { for (kk=k-1;kk<=k+1;kk++) { if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<=zgrid) { llist = grid[ii][jj][kk]; while (llist) { if (llist->atom->name[0]=='N' && llist->atom->name[1]==' ' && (abs(llist->atom->res->locnum-n_atom1->res->locnum)>2)) { tmp_atom = llist->atom; if (tmp_atom->res!=c_atom2->res) { dx = o_atom1->x-tmp_atom->x; dy = o_atom1->y-tmp_atom->y; dz = o_atom1->z-tmp_atom->z; dist = dx*dx+dy*dy+dz*dz; if (dist<min_dist2 && dist<25.0) { n_atom2=tmp_atom; min_dist2 = dist; } } } llist = llist->next; } } } } } } if (n_atom2) { c_atom2 = o_atom2 = NULL; atom = n_atom2->res->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]==' ') { c_atom2 = atom; } if (atom->name[0]=='O' && atom->name[1]==' ') { c_atom2 = atom; } atom = atom->next; } if (c_atom2) { hx1 = o_atom1->x-c_atom1->x; hy1 = o_atom1->y-c_atom1->y; hz1 = o_atom1->z-c_atom1->z; dd = -1.081f/sqrt(hx1*hx1+hy1*hy1+hz1*hz1); hx1 *= dd; hy1 *= dd; hz1 *= dd; hx1 += n_atom1->x; hy1 += n_atom1->y; hz1 += n_atom1->z; } *******/ return ene; } // rotates a point around a vector void rot_point_vector(real *x, real *y, real *z, real u, real v, real w, real angle) { real ux, uy, uz, vx, vy, vz, wx, wy, wz, sa, ca; sa = sinf(10.0*M_PI*angle/180.0); ca = cosf(10.0*M_PI*angle/180.0); ux = u**x; uy = u**y; uz = u**z; vx = v**x; vy = v**y; vz = v**z; wx = w**x; wy = w**y; wz = w**z; *x = u*(ux+vy+wz)+(*x*(v*v+w*w)-u*(vy+wz))*ca+(-wy+vz)*sa; *y = v*(ux+vy+wz)+(*y*(u*u+w*w)-v*(ux+wz))*ca+( wx-uz)*sa; *z = w*(ux+vy+wz)+(*z*(u*u+v*v)-w*(ux+vy))*ca+(-vx+uy)*sa; } // rotates a peptide plate void rot_peptide(res_type *res, real angle) { atom_type *atom, *c_atom, *o_atom, *n_atom, *ca_atom1, *ca_atom2; real u, v, w, x, y, z, dd; if (!res || !res->prev) return; c_atom = o_atom = n_atom = ca_atom1 = ca_atom2 = NULL; atom = res->prev->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]=='A') ca_atom1 = atom; if (atom->name[0]=='C' && atom->name[1]==' ') c_atom = atom; if (atom->name[0]=='O' && atom->name[1]==' ') o_atom = atom; atom = atom->next; } atom = res->atoms; while (atom) { if (atom->name[0]=='C' && atom->name[1]=='A') ca_atom2 = atom; if (atom->name[0]=='N' && atom->name[1]==' ') n_atom = atom; atom = atom->next; } if (c_atom && o_atom && n_atom && ca_atom1 && ca_atom2) { u = ca_atom2->x-ca_atom1->x; v = ca_atom2->y-ca_atom1->y; w = ca_atom2->z-ca_atom1->z; dd = 1.0f/sqrt(u*u+v*v+w*w); u*=dd; v*=dd; w*=dd; // normalize ca-ca vector x = n_atom->x-ca_atom1->x; y = n_atom->y-ca_atom1->y; z = n_atom->z-ca_atom1->z; rot_point_vector(&x, &y, &z, u, v, w, angle); n_atom->x = x+ca_atom1->x; n_atom->y = y+ca_atom1->y; n_atom->z = z+ca_atom1->z; x = c_atom->x-ca_atom1->x; y = c_atom->y-ca_atom1->y; z = c_atom->z-ca_atom1->z; rot_point_vector(&x, &y, &z, u, v, w, angle); c_atom->x = x+ca_atom1->x; c_atom->y = y+ca_atom1->y; c_atom->z = z+ca_atom1->z; x = o_atom->x-ca_atom1->x; y = o_atom->y-ca_atom1->y; z = o_atom->z-ca_atom1->z; rot_point_vector(&x, &y, &z, u, v, w, angle); o_atom->x = x+ca_atom1->x; o_atom->y = y+ca_atom1->y; o_atom->z = z+ca_atom1->z; } } void optimize_backbone(mol_type *chain) { int xgrid, ygrid, zgrid; atom_list ****grid; atom_type *atom; res_type *res; real ene, min_ene, tot1, tot2; int i, k, best; FILE *out; if (_VERBOSE) printf("Optimizing backbone...\n"); allocate_grid(&grid, &xgrid, &ygrid, &zgrid); tot1 = tot2 = 0.0; res = chain->residua; while (res) { ene = hb_energy(res, grid, xgrid, ygrid, zgrid); if (ene<-0.5) tot1 += ene; res = res->next; } res = chain->residua; while (res) { if (res->type!=7) { ene = hb_energy(res, grid, xgrid, ygrid, zgrid); if (ene<1.0) { // try to optimize min_ene = ene; rot_peptide(res, -1.1); best = 0; for (i=-10;i<10;i++) { rot_peptide(res, 0.1); ene = hb_energy(res, grid, xgrid, ygrid, zgrid); if (ene<min_ene) { best = i; min_ene = ene; } } rot_peptide(res,-0.9); ene = hb_energy(res, grid, xgrid, ygrid, zgrid); if (min_ene<ene) { rot_peptide(res,0.1*best); ene = hb_energy(res, grid, xgrid, ygrid, zgrid); } } } res = res->next; } res = chain->residua; while (res) { ene = hb_energy(res, grid, xgrid, ygrid, zgrid); if (ene<-0.5) tot2 += ene; res = res->next; } if (_VERBOSE) printf("Backbone HB energy: before %g, after: %g, difference: %g\n", tot1, tot2, tot2-tot1); } int main(int argc, char **argv) { int i, j, next; char buf[100]; char *name=NULL, *ini_name=NULL; char *ptr, out_name[1000]; real f; mol_type *mol; struct timeb time0, time1; for (i=1; i<argc; i++) { if (argv[i][0]=='-') { next=0; for (j=1; j<(int)strlen(argv[i]); j++) { switch(argv[i][j]) { case 'v': _VERBOSE=1; break; case 'c': _CA_OPTIMIZE=0; break; case 'e': _BB_REARRANGE=1; break; case 'r': _CA_RANDOM=1; break; case 'z': _CHIRAL=0; break; case 't': _CA_TRAJECTORY=1; break; case 'n': _CENTER_CHAIN=1; break; case 'b': _REBUILD_BB=0; break; case 's': _REBUILD_SC=0; break; case 'i': ini_name = argv[++i]; next=1; break; case 'g': _PDB_SG=1; break; case 'x': _TIME_SEED=1; break; case 'o': _XVOLUME=0; break; case 'h': _REBUILD_H=0; break; case 'q': _BB_OPTIMIZE=1; break; case 'p': _CISPRO=1; break; case 'u': if (sscanf(argv[++i],"%lf",&f)==1) { _CA_START_DIST = f; } next=1; break; default: { printf("Unknown option: %c\n", argv[i][j]); return -1; } } if (next) break; } } else { if (!name) name=argv[i]; } } if (!name) { printf("PULCHRA Protein Chain Restoration Algorithm version %4.2f\n", PULCHRA_VERSION); printf("Usage: %s [options] <pdb_file>\n", argv[0]); printf("The program default input is a PDB file.\n"); printf("Output file <pdb_file.rebuild.pdb> will be created as a result.\n"); printf("Valid options are:\n\n"); printf(" -v : verbose output (default: off)\n"); printf(" -n : center chain (default: off)\n"); printf(" -x : time-seed random number generator (default: off)\n"); printf(" -g : use PDBSG as an input format (CA=C-alpha, SC or CM=side chain c.m.)\n\n"); printf(" -c : skip C-alpha positions optimization (default: on)\n"); printf(" -p : detect cis-prolins (default: off)\n"); printf(" -r : start from a random chain (default: off)\n"); printf(" -i pdbfile : read the initial C-alpha coordinates from a PDB file\n"); printf(" -t : save chain optimization trajectory to file <pdb_file.pdb.trajectory>\n"); printf(" -u value : maximum shift from the restraint coordinates (default: 0.5A)\n\n"); printf(" -e : rearrange backbone atoms (C, O are output after side chain) (default: off)\n"); #ifdef COMPILE_BB printf(" -b : skip backbone reconstruction (default: on)\n"); printf(" -q : optimize backbone hydrogen bonds pattern (default: off)\n"); printf(" -h : outputs hydrogen atoms (default: off)\n"); #endif #ifdef COMPILE_ROT printf(" -s : skip side chains reconstruction (default: on)\n"); printf(" -o : don't attempt to fix excluded volume conflicts (default: on)\n"); printf(" -z : don't check amino acid chirality (default: on)\n"); #endif printf("\n"); return -1; } for (i=0; i<255; i++) /* prepare hash table*/ AA_NUMS[i] = 20; /* dummy aa code */ for (i=0; i<20; i++) AA_NUMS[(int)SHORT_AA_NAMES[i]] = i; setbuf(stdout,0); if (_TIME_SEED) srand(time(NULL)); else srand(1234); if (_VERBOSE) printf("PULCHRA Protein Chain Restoration Algorithm version %4.2f\n", PULCHRA_VERSION); ftime(&time0); chain = new_mol(); if (read_pdb_file(name,chain,"chain")==FILE_NOT_FOUND) { if (_VERBOSE) printf("Can't read the input file!\n"); return -1; } if (_VERBOSE) printf("%d residua read.\n", chain->nres); chain_length = chain->nres; if (_CA_OPTIMIZE) { snprintf(out_name,1000,"%s.tra",name); ca_optimize(out_name, ini_name); } #ifdef COMPILE_BB if (_REBUILD_BB) { rebuild_backbone(); if (_BB_OPTIMIZE) { optimize_backbone(chain); } } #endif #ifdef COMPILE_ROT if (_REBUILD_BB && _REBUILD_SC) { rebuild_sidechains(); if (_XVOLUME) optimize_exvol(); if (_CHIRAL) chirality_check(); } #endif if (_CENTER_CHAIN) { center_chain(chain); } if (_BB_REARRANGE) { if (_VERBOSE) printf("Rearranging backbone atoms...\n"); } ptr = strstr(name,".pdb"); if (ptr) ptr[0]=0; snprintf(out_name,1000,"%s.rebuilt.pdb",name); if (_VERBOSE) printf("Writing output file %s...\n", out_name); write_pdb(out_name, chain); ftime(&time1); if (_VERBOSE) printf("Done. Reconstruction finished in %.3f s.\n", (real)0.001*(1000.*(time1.time-time0.time)+(time1.millitm-time0.millitm))); return 0; }