Mercurial > repos > guerler > springsuite

comparison spring_package/pulchra/pulchra.c @ 17:c790d25086dc draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit b0ede77caf410ab69043d33a44e190054024d340-dirty"
author guerler
date Wed, 28 Oct 2020 05:11:56 +0000
parents
children
comparison
equal deleted inserted replaced
16:16eb2acaaa20 17:c790d25086dc
1
2 //
3 // PULCHRA
4 // Protein Chain Restoration Algorithm
5 //
6 // Version 3.04
7 // July 2007
8 // Contact: Piotr Rotkiewicz, piotr -at- pirx -dot- com
9 //
10 // to compile:
11 // cc -O3 pulchra.c pulchra_data.c -lm -o pulchra
12 //
13 // to use:
14 // ./pulchra file.pdb
15 //
16 // to display available options:
17 // ./pulchra
18 //
19
20 #define COMPILE_BB
21 #define COMPILE_ROT
22
23 #include <math.h>
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <string.h>
27 #include <sys/timeb.h>
28 #include <time.h>
29
30 #define uchar unsigned char
31 #define uint unsigned int
32 #define real double
33
34 #include "pulchra_common.h"
35
36 #define PULCHRA_VERSION 3.04
37 #define MAX_BUF_SIZE 1000
38
39 #define FILE_SUCCESS 0
40 #define FILE_NOT_FOUND -1
41 #define FILE_WARNING -2
42
43 #define FATAL_MAE -1
44
45 #define FLAG_BACKBONE 1
46 #define FLAG_CALPHA 2
47 #define FLAG_SIDECHAIN 4
48 #define FLAG_SCM 8
49
50 #define FLAG_PROTEIN 1
51 #define FLAG_DNA 2
52 #define FLAG_RNA 4
53 #define FLAG_CHYDRO 8
54
55 #define RADDEG 180./M_PI
56 #define DEGRAD M_PI/180.
57
58 int _VERBOSE = 0;
59 int _BB_REARRANGE = 1;
60 int _BB_OPTIMIZE = 0;
61 int _CA_OPTIMIZE = 1;
62 int _CA_RANDOM = 0;
63 int _CA_ITER = 100;
64 int _CA_TRAJECTORY = 0;
65 int _CISPRO = 0;
66 int _CHIRAL = 1;
67 int _CENTER_CHAIN = 0;
68 int _REBUILD_BB = 1;
69 int _REBUILD_SC = 1;
70 int _REBUILD_H = 0;
71 int _PDB_SG = 0;
72 int _TIME_SEED = 0;
73 int _XVOLUME = 1;
74 int _XVOL_ITER = 3;
75 real _CA_START_DIST = 3.0;
76 real _CA_XVOL_DIST = 3.5;
77 real _SG_XVOL_DIST = 1.6;
78
79 #define CALC_C_ALPHA
80 #define CALC_C_ALPHA_ANGLES
81 #define CALC_C_ALPHA_START
82 #define CALC_C_ALPHA_XVOL
83
84 real CA_K=10.0;
85 real CA_ANGLE_K=20.0;
86 real CA_START_K=0.01;
87 real CA_XVOL_K=10.00;
88
89 #define CA_DIST 3.8
90 #define CA_DIST_TOL 0.1
91 #define CA_DIST_CISPRO 2.9
92 #define CA_DIST_CISPRO_TOL 0.1
93 #define E_EPS 1e-10
94
95 // grid resolution (used for fast clash detection)
96 #define GRID_RES 6.0
97
98 int chain_length = 0;
99
100 char AA_NAMES[21][4] =
101 { "GLY", "ALA", "SER", "CYS", "VAL",
102 "THR", "ILE", "PRO", "MET", "ASP",
103 "ASN", "LEU", "LYS", "GLU", "GLN",
104 "ARG", "HIS", "PHE", "TYR", "TRP",
105 "UNK" };
106
107 char SHORT_AA_NAMES[22] = { "GASCVTIPMDNLKEQRHFYWX" };
108
109 int AA_NUMS[256];
110
111 int nheavy[20] = { 0, 1, 2, 2, 3, 3, 4, 3, 4, 4, 4, 4, 5, 5, 5, 7, 6, 7, 8, 10};
112
113 char *backbone_atoms[4] = { "N ", "CA ", "C ", "O " };
114
115 char *heavy_atoms[200]= {
116 /* GLY */ NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
117 /* ALA */ "CB ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
118 /* SER */ "CB ", "OG ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
119 /* CYS */ "CB ", "SG ", NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL,
120 /* VAL */ "CB ", "CG1", "CG2", NULL, NULL, NULL, NULL, NULL, NULL, NULL,
121 /* THR */ "CB ", "OG1", "CG2", NULL, NULL, NULL, NULL, NULL, NULL, NULL,
122 /* ILE */ "CB ", "CG1", "CG2", "CD1", NULL, NULL, NULL, NULL, NULL, NULL,
123 /* PRO */ "CB ", "CG ", "CD ", NULL, NULL, NULL, NULL, NULL, NULL, NULL,
124 /* MET */ "CB ", "CG ", "SD ", "CE ", NULL, NULL, NULL, NULL, NULL, NULL,
125 /* ASP */ "CB ", "CG ", "OD1", "OD2", NULL, NULL, NULL, NULL, NULL, NULL,
126 /* ASN */ "CB ", "CG ", "OD1", "ND2", NULL, NULL, NULL, NULL, NULL, NULL,
127 /* LEU */ "CB ", "CG ", "CD1", "CD2", NULL, NULL, NULL, NULL, NULL, NULL,
128 /* LYS */ "CB ", "CG ", "CD ", "CE ", "NZ ", NULL, NULL, NULL, NULL, NULL,
129 /* GLU */ "CB ", "CG ", "CD ", "OE1", "OE2", NULL, NULL, NULL, NULL, NULL,
130 /* GLN */ "CB ", "CG ", "CD ", "OE1", "NE2", NULL, NULL, NULL, NULL, NULL,
131 /* ARG */ "CB ", "CG ", "CD ", "NE ", "CZ ", "NH1", "NH2", NULL, NULL, NULL,
132 /* HIS */ "CB ", "CG ", "ND1", "CD2", "CE1", "NE2", NULL, NULL, NULL, NULL,
133 /* PHE */ "CB ", "CG ", "CD1", "CD2", "CE1", "CE2", "CZ ", NULL, NULL, NULL,
134 /* TYR */ "CB ", "CG ", "CD1", "CD2", "CE1", "CE2", "CZ ", "OH ", NULL, NULL,
135 /* TRP */ "CB ", "CG ", "CD1", "CD2", "NE1", "CE2", "CE3", "CZ2", "CZ3", "CH2"};
136
137 /* reads full-atom pdb file */
138
139 struct _res_type;
140
141 typedef struct _atom_type {
142 struct _atom_type *next;
143 real x, y, z;
144 char *name;
145 int num, locnum;
146 int flag;
147 char cispro;
148 int gx, gy, gz;
149 struct _res_type *res;
150 struct _atom_type *prev;
151 } atom_type;
152
153 typedef struct _res_type {
154 struct _res_type *next;
155 atom_type *atoms;
156 int num, locnum, natoms;
157 int type;
158 char pdbsg;
159 char protein;
160 char *name;
161 char chain;
162 real sgx, sgy, sgz;
163 real cmx, cmy, cmz;
164 struct _res_type *prev;
165 } res_type;
166
167 typedef struct _mol_type {
168 struct _mol_type *next;
169 res_type *residua;
170 int nres;
171 unsigned char *r14;
172 char *name;
173 uchar *seq;
174 char **contacts;
175 real **cutoffs;
176 struct _mol_type *prev;
177 } mol_type;
178
179 #define MIN(a,b) (a<b?a:b)
180 #define MAX(a,b) (a>b?a:b)
181
182 mol_type *chain = NULL;
183
184
185 real rnd(void)
186 {
187 return 0.001*(real)(rand()%1000);
188 }
189
190 atom_type *new_atom(void)
191 {
192 atom_type *tmpatom;
193
194 tmpatom = (atom_type*) calloc(sizeof(atom_type),1);
195 if (tmpatom) {
196 tmpatom->x=tmpatom->y=tmpatom->z=0.;
197 tmpatom->name=NULL;
198 tmpatom->num=tmpatom->locnum=tmpatom->flag=0;
199 tmpatom->next=tmpatom->prev=NULL;
200 }
201
202 return tmpatom;
203 }
204
205 res_type* new_res(void)
206 {
207 res_type *tmpres;
208
209 tmpres = (res_type*) calloc(sizeof(res_type),1);
210 if (tmpres) {
211 tmpres->num=0;
212 tmpres->name=NULL;
213 tmpres->atoms=NULL;
214 tmpres->chain=' ';
215 tmpres->next=tmpres->prev=NULL;
216 }
217
218 return tmpres;
219 }
220
221 mol_type *new_mol(void)
222 {
223 mol_type *tmpmol;
224
225 tmpmol = (mol_type*) calloc(sizeof(mol_type),1);
226 if (tmpmol) {
227 tmpmol->name=NULL;
228 tmpmol->residua=NULL;
229 tmpmol->next=tmpmol->prev=NULL;
230 }
231
232 return tmpmol;
233 }
234
235 void add_atom(atom_type* atomlist, atom_type* newatom)
236 {
237 atom_type *tmpatom;
238
239 if (!atomlist)
240 atomlist=newatom;
241 else {
242 tmpatom=atomlist->next;
243 atomlist->next=newatom;
244 newatom->prev=atomlist;
245 newatom->next=tmpatom;
246 if (tmpatom) tmpatom->prev=newatom;
247 }
248 }
249
250 void add_res(res_type* reslist, res_type* newres)
251 {
252 res_type *tmpres;
253
254 if (!reslist)
255 reslist=newres;
256 else {
257 tmpres=reslist->next;
258 reslist->next=newres;
259 newres->prev=reslist;
260 newres->next=tmpres;
261 if (tmpres) tmpres->prev=newres;
262 }
263 }
264
265 void add_mol(mol_type* mollist, mol_type* newmol)
266 {
267 mol_type *tmpmol;
268
269 if (!mollist)
270 mollist=newmol;
271 else {
272 tmpmol=mollist->next;
273 mollist->next=newmol;
274 newmol->prev=mollist;
275 newmol->next=tmpmol;
276 if (tmpmol) tmpmol->prev=newmol;
277 }
278 }
279
280 void delete_atom(atom_type* atom)
281 {
282 atom_type *tmpatom;
283
284 if (atom->prev) atom->prev->next=atom->next;
285 if (atom->next) atom->next->prev=atom->prev;
286 if (atom->name) free(atom->name);
287 free(atom);
288 atom=NULL;
289 }
290
291 void delete_res(res_type* res)
292 {
293 res_type *tmpres;
294 atom_type *tmpatom;
295
296 if (res->prev) res->prev->next=res->next;
297 if (res->next) res->next->prev=res->prev;
298 if (res->name) free(res->name);
299 if (res->atoms) {
300 while (res->atoms) {
301 tmpatom = res->atoms->next;
302 delete_atom(res->atoms);
303 res->atoms=tmpatom;
304 }
305 }
306 free(res);
307 res=NULL;
308 }
309
310 void delete_mol(mol_type* mol)
311 {
312 mol_type *tmpmol;
313 res_type *tmpres;
314 int i;
315
316 if (mol->prev) mol->prev->next=mol->next;
317 if (mol->next) mol->next->prev=mol->prev;
318 if (mol->name) free(mol->name);
319 if (mol->residua) {
320 while (mol->residua) {
321 tmpres = mol->residua->next;
322 delete_res(mol->residua);
323 mol->residua=tmpres;
324 }
325 }
326 if (mol->contacts) {
327 for (i=0; i<mol->nres; i++) free(mol->contacts[i]);
328 free(mol->contacts);
329 }
330 if (mol->cutoffs) {
331 for (i=0; i<mol->nres; i++) free(mol->cutoffs[i]);
332 free(mol->cutoffs);
333 }
334 free(mol);
335 mol=NULL;
336 }
337
338
339 atom_type* get_last_atom(atom_type* atom)
340 {
341 while (atom->next) atom=atom->next;
342
343 return atom;
344 }
345
346 res_type* get_last_res(res_type* res)
347 {
348 while (res->next) res=res->next;
349
350 return res;
351 }
352
353 mol_type *get_last_mol(mol_type* mol)
354 {
355 while (mol->next) mol=mol->next;
356
357 return mol;
358 }
359
360 char setseq(char* aaname)
361 {
362 int i;
363
364 for (i=0; i<21; i++) {
365 if ((aaname[0]==AA_NAMES[i][0]) &&
366 (aaname[1]==AA_NAMES[i][1]) &&
367 (aaname[2]==AA_NAMES[i][2]))
368 break;
369 }
370
371 if (i==21) {
372 if (!strcmp(aaname, "GLX"))
373 return 'E';
374 if (!strcmp(aaname, "ASX"))
375 return 'D';
376 if (!strcmp(aaname, "HID"))
377 return 'H';
378 if (!strcmp(aaname, "MSE"))
379 return 'M';
380 if (!strcmp(aaname, "SEP"))
381 return 'S';
382 if (!strcmp(aaname, "TPO"))
383 return 'T';
384 if (!strcmp(aaname, "PTR"))
385 return 'Y';
386 i--;
387 }
388
389 return SHORT_AA_NAMES[i];
390 }
391
392 int orient(res_type *res1, res_type *res2)
393 {
394 real x1, y1, z1;
395 real x2, y2, z2;
396 real cax, cay, caz;
397 real len, vect, angle;
398 atom_type *atom;
399
400 if (!res1 || !res2) return 0;
401
402 atom=res1->atoms;
403 cax=cay=caz=0.;
404 while (atom) {
405 if (!strncmp(atom->name,"CA",2)) {
406 cax=atom->x; cay=atom->y; caz=atom->z;
407 }
408 atom=atom->next;
409 }
410 x1=res1->sgx-cax; y1=res1->sgy-cay; z1=res1->sgz-caz;
411 if (x1==0. && y1==0. && z1==0.) x1+=1.0;
412
413 atom=res2->atoms;
414 cax=cay=caz=0.;
415 while (atom) {
416 if (!strncmp(atom->name,"CA",2)) {
417 cax=atom->x; cay=atom->y; caz=atom->z;
418 }
419 atom=atom->next;
420 }
421 x2=res2->sgx-cax; y2=res2->sgy-cay; z2=res2->sgz-caz;
422 if (x2==0. && y2==0. && z2==0.) x2+=1.0;
423
424 vect = x1*x2+y1*y2+z1*z2;
425 len = sqrt(x1*x1+y1*y1+z1*z1)*sqrt(x2*x2+y2*y2+z2*z2);
426 if (len) vect /= len;
427
428 angle=RADDEG*acos(vect);
429
430 if (angle>120.) return 1; /*anti*/
431 if (angle>60.) return 2; /*mid*/
432
433 return 3; /*par*/
434 }
435
436 int res_contact(res_type *res1, res_type *res2) {
437 atom_type *atoms1, *atoms2;
438 real dx, dy, dz;
439
440 atoms1 = res1->atoms;
441 while (atoms1) {
442 atoms2 = res2->atoms;
443 while (atoms2) {
444 dx=atoms1->x-atoms2->x;
445 dy=atoms1->y-atoms2->y;
446 dz=atoms1->z-atoms2->z;
447 if ((atoms1->flag & FLAG_SIDECHAIN) && (atoms2->flag & FLAG_SIDECHAIN) && (dx*dx+dy*dy+dz*dz<20.25)) {
448 return 1;
449 }
450 atoms2=atoms2->next;
451 }
452 atoms1=atoms1->next;
453 }
454
455 return 0;
456 }
457
458 int read_pdb_file(char* filename, mol_type* molecules, char *realname)
459 {
460 FILE *inp;
461 char buffer[1000];
462 char atmname[10];
463 char resname[10];
464 char version;
465 int prevresnum, resnum, atmnum, locatmnum, num, locnum=0, i, j;
466 atom_type *prevatom1, *prevatom2, *prevatom3, *prevatom4;
467 int sgnum, cc, nres, ok, natom;
468 real sgx, sgy, sgz;
469 res_type *res, *test1, *test2;
470 atom_type *atom;
471 real x, y, z;
472 real dist;
473 unsigned char bin;
474 int warn=0;
475 real cutoff;
476
477 if (_VERBOSE) printf("Reading input file %s...\n", filename);
478
479 inp = fopen(filename, "r");
480 if (!inp) {
481 if (_VERBOSE) printf("ERROR: can't open %s !!!\n", filename);
482 return FILE_NOT_FOUND;
483 }
484
485 molecules->nres=0;
486 molecules->name=(char*)calloc(strlen(realname)+1,1);
487 strcpy(molecules->name, realname);
488
489 atmname[3]=0;
490 resname[3]=0;
491 prevresnum=-666;
492 locatmnum=0;
493 sgnum=0;
494 sgx=sgy=sgz=0.;
495 res=NULL;
496 while (!feof(inp)) {
497 if (fgets(buffer, 1000, inp)!=buffer) break;
498 if (!strncmp(buffer, "END", 3) || !strncmp(buffer, "TER", 3)) break; // end of file; only single molecule read
499 if (!strncmp(buffer, "ATOM", 4) || !strncmp(buffer, "HETATM", 6)) {
500 if (buffer[16]!=' ' && buffer[16]!='A') continue;
501 sscanf(&buffer[22], "%d", &resnum);
502 strncpy(resname, &buffer[17], 3);
503 strncpy(atmname, &buffer[13], 3);
504 if (resnum==prevresnum && !strncmp(atmname, "N ", 2)) {
505 if (_VERBOSE) printf("WARNING: fault in numeration at residuum %s[%d]\n", resname, resnum);
506 warn=1;
507 }
508 if (atmname[0]=='H') continue;
509 if (resnum!=prevresnum || !strncmp(atmname, "N ", 2)) {
510 prevresnum=resnum;
511 if (res) {
512 if (sgnum) {
513 res->sgx=sgx/(real)sgnum;
514 res->sgy=sgy/(real)sgnum;
515 res->sgz=sgz/(real)sgnum;
516 } else {
517 res->sgx=res->sgy=res->sgz=0.;
518 }
519 }
520 locatmnum=0;
521 version=' ';
522 res = new_res();
523 sgnum=0;
524 sgx=sgy=sgz=0.;
525 molecules->nres++;
526 res->name = calloc(strlen(resname)+1, 1);
527 res->type = AA_NUMS[setseq(resname)];
528 res->locnum=locnum++;
529 res->num = resnum;
530 res->natoms=0;
531 res->chain = buffer[21];
532 strcpy(res->name, resname);
533 if (molecules->residua) {
534 add_res(get_last_res(molecules->residua), res);
535 } else {
536 molecules->residua = res;
537 }
538 }
539 atom = new_atom();
540 atom->res = res;
541 res->natoms++;
542 locatmnum++;
543 sscanf(&buffer[7], "%d", &atmnum);
544 sscanf(&buffer[30], "%lf", &x);
545 sscanf(&buffer[38], "%lf", &y);
546 sscanf(&buffer[46], "%lf", &z);
547 version = buffer[16];
548 atom->name = calloc(strlen(atmname)+1,1);
549 strcpy(atom->name, atmname);
550 atom->x=x; atom->y=y; atom->z=z;
551 atom->num = atmnum;
552 atom->locnum = locatmnum;
553 if ((atmname[0]=='S' && atmname[1]=='C')||(atmname[0]=='C' && atmname[1]=='M')) {
554 res->cmx = x;
555 res->cmy = y;
556 res->cmz = z;
557 res->pdbsg=1;
558 if (res->type<20) {
559 res->protein=1;
560 }
561 } else
562 if (!( ((atmname[0]=='C' || atmname[0]=='N' || atmname[0]=='O') && atmname[1]==' ') ||
563 (atmname[0]=='H') ||
564 (atmname[0]=='C' && atmname[1]=='A') ||
565 (atmname[0]=='O' && atmname[1]=='X' && atmname[2]=='T') ) ) {
566 sgx+=x;
567 sgy+=y;
568 sgz+=z;
569 sgnum++;
570 atom->flag |= FLAG_SIDECHAIN;
571 } else
572 atom->flag |= FLAG_BACKBONE;
573 if (atmname[0]=='C' && atmname[1]=='A') {
574 atom->flag |= FLAG_BACKBONE;
575 if (res->type<20) {
576 res->protein=1;
577 }
578 if (!res->pdbsg) {
579 res->cmx = x;
580 res->cmy = y;
581 res->cmz = z;
582 }
583 }
584 if (atmname[0]=='C' && atmname[1]=='M') {
585 atom->flag |= FLAG_SCM;
586 }
587 if (atmname[0]=='S' && atmname[1]=='C') {
588 atom->flag |= FLAG_SCM;
589 }
590 if (res->atoms) {
591 add_atom(get_last_atom(res->atoms), atom);
592 } else {
593 res->atoms = atom;
594 }
595 }
596 }
597
598 if (res) {
599 if (sgnum) {
600 res->sgx=sgx/(real)sgnum;
601 res->sgy=sgy/(real)sgnum;
602 res->sgz=sgz/(real)sgnum;
603 } else {
604 res->sgx=res->sgy=res->sgz=0.;
605 }
606 }
607 fclose(inp);
608
609 molecules->seq = (uchar*)calloc(sizeof(uchar)*molecules->nres+1,1);
610 res=molecules->residua; i=0;
611 while (res) {
612 molecules->seq[i++]=(uchar)AA_NUMS[(int)setseq(res->name)];
613 res=res->next;
614 }
615
616 if (!warn) return FILE_SUCCESS; else return FILE_WARNING;
617 }
618
619 #define bool int
620 #define true 1
621 #define false 0
622
623 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)
624 {
625 int i, j;
626 real dx, dy, dz;
627 real dist, ddist, ddist2;
628 real new_e_pot;
629 real theta0, tdif, th, aa, bb, ab;
630 real ff0, ff2, dth, m0, m2, grad, f0[3], f2[3];
631 real adiff[3], bdiff[3];
632 real deriv, theta, dtheta, len1, len2, cos_theta, sin_theta;
633 real dx1, dy1, dz1;
634 real dx2, dy2, dz2;
635 real dx3, dy3, dz3;
636 real vx1, vy1, vz1;
637 real vx2, vy2, vz2;
638 real vx3, vy3, vz3;
639
640 real r12x, r12y, r12z;
641 real r32x, r32y, r32z;
642 real d12, d32, d12inv, d32inv, c1, c2, diff;
643 real f1x, f1y, f1z;
644 real f2x, f2y, f2z;
645 real f3x, f3y, f3z;
646
647 for (i=0; i<chain_length; i++) {
648 new_c_alpha[i][0]=c_alpha[i]->x+alpha*gradient[i][0];
649 new_c_alpha[i][1]=c_alpha[i]->y+alpha*gradient[i][1];
650 new_c_alpha[i][2]=c_alpha[i]->z+alpha*gradient[i][2];
651 }
652
653 new_e_pot = 0.0;
654
655 ene[0]=ene[1]=ene[2]=ene[3]=0.0;
656
657 for (i=0; i<chain_length; i++) {
658 #ifdef CALC_C_ALPHA_START
659 dx=new_c_alpha[i][0]-init_c_alpha[i][0];
660 dy=new_c_alpha[i][1]-init_c_alpha[i][1];
661 dz=new_c_alpha[i][2]-init_c_alpha[i][2];
662 dist=sqrt(dx*dx+dy*dy+dz*dz);
663 ddist = -dist;
664 if (dist>_CA_START_DIST) {
665 ddist2=dist*dist;
666 new_e_pot+=CA_START_K*ddist2;
667 ene[1] += CA_START_K*ddist2;
668 if (calc_gradient) {
669 grad = ddist * (-2.0*CA_START_K)/dist;
670 gradient[i][0]-=grad*dx;
671 gradient[i][1]-=grad*dy;
672 gradient[i][2]-=grad*dz;
673 }
674 }
675
676 #endif
677
678
679 #ifdef CALC_C_ALPHA
680 if (i>0) {
681 dx=new_c_alpha[i][0]-new_c_alpha[i-1][0];
682 dy=new_c_alpha[i][1]-new_c_alpha[i-1][1];
683 dz=new_c_alpha[i][2]-new_c_alpha[i-1][2];
684 dist=sqrt(dx*dx+dy*dy+dz*dz);
685 if (c_alpha[i]->cispro) {
686 ddist=CA_DIST_CISPRO-dist;
687 // if (fabs(ddist)<CA_DIST_CISPRO_TOL) ddist=0.0;
688 } else {
689 ddist=CA_DIST-dist;
690 // if (fabs(ddist)<CA_DIST_TOL) ddist=0.0;
691 }
692 ddist2=ddist*ddist;
693 new_e_pot+=CA_K*ddist2;
694 ene[0] += CA_K*ddist2;
695 if (calc_gradient) {
696 grad = ddist * (-2.0*CA_K)/dist;
697 gradient[i][0]-=grad*dx;
698 gradient[i][1]-=grad*dy;
699 gradient[i][2]-=grad*dz;
700 gradient[i-1][0]+=grad*dx;
701 gradient[i-1][1]+=grad*dy;
702 gradient[i-1][2]+=grad*dz;
703 }
704 }
705 #endif
706
707 #ifdef CALC_C_ALPHA_XVOL
708 for (j=0;j<i;j++) {
709 if (abs(i-j)>2) {
710 dx=new_c_alpha[i][0]-new_c_alpha[j][0];
711 dy=new_c_alpha[i][1]-new_c_alpha[j][1];
712 dz=new_c_alpha[i][2]-new_c_alpha[j][2];
713 dist=sqrt(dx*dx+dy*dy+dz*dz);
714 ddist = dist-_CA_XVOL_DIST;
715 if (dist<_CA_XVOL_DIST) {
716 ddist2 = dist*dist;
717 new_e_pot+=CA_XVOL_K*ddist2;
718 ene[3] += CA_XVOL_K*ddist2;
719 if (calc_gradient) {
720 grad = ddist*(8.0*CA_XVOL_K)/dist;
721 gradient[i][0]-=grad*dx;
722 gradient[i][1]-=grad*dy;
723 gradient[i][2]-=grad*dz;
724 gradient[j][0]+=grad*dx;
725 gradient[j][1]+=grad*dy;
726 gradient[j][2]+=grad*dz;
727 }
728 }
729 }
730 }
731 #endif
732
733 #ifdef CALC_C_ALPHA_ANGLES
734
735 if (i>0 && i<chain_length-1) {
736 r12x=new_c_alpha[i-1][0]-new_c_alpha[i][0];
737 r12y=new_c_alpha[i-1][1]-new_c_alpha[i][1];
738 r12z=new_c_alpha[i-1][2]-new_c_alpha[i][2];
739 r32x=new_c_alpha[i+1][0]-new_c_alpha[i][0];
740 r32y=new_c_alpha[i+1][1]-new_c_alpha[i][1];
741 r32z=new_c_alpha[i+1][2]-new_c_alpha[i][2];
742 d12 = sqrt(r12x*r12x+r12y*r12y+r12z*r12z);
743 d32 = sqrt(r32x*r32x+r32y*r32y+r32z*r32z);
744 cos_theta = (r12x*r32x+r12y*r32y+r12z*r32z)/(d12*d32);
745 if (cos_theta>1.0)
746 cos_theta = 1.0;
747 else
748 if (cos_theta<-1.0)
749 cos_theta = -1.0;
750 sin_theta = sqrt(1.0-cos_theta*cos_theta);
751 theta = acos(cos_theta);
752
753 if (RADDEG*theta<80.)
754 diff = theta-80.*DEGRAD;
755 else
756 if (RADDEG*theta>150.)
757 diff = theta-150.*DEGRAD;
758 else
759 diff = 0.0;
760
761 new_e_pot += CA_ANGLE_K*diff*diff;
762 ene[2] += CA_ANGLE_K*diff*diff;
763 if (calc_gradient) {
764 d12inv = 1.0/d12;
765 d32inv = 1.0/d32;
766 diff *= (-2.0*CA_ANGLE_K)/sin_theta;
767 c1 = diff*d12inv;
768 c2 = diff*d32inv;
769 f1x = c1*(r12x*(d12inv*cos_theta)-r32x*d32inv);
770 f1y = c1*(r12y*(d12inv*cos_theta)-r32y*d32inv);
771 f1z = c1*(r12z*(d12inv*cos_theta)-r32z*d32inv);
772 f3x = c2*(r32x*(d32inv*cos_theta)-r12x*d12inv);
773 f3y = c2*(r32y*(d32inv*cos_theta)-r12y*d12inv);
774 f3z = c2*(r32z*(d32inv*cos_theta)-r12z*d12inv);
775 f2x = -f1x-f3x;
776 f2y = -f1y-f3y;
777 f2z = -f1z-f3z;
778 gradient[i-1][0]+=f1x;
779 gradient[i-1][1]+=f1y;
780 gradient[i-1][2]+=f1z;
781 gradient[i][0]+=f2x;
782 gradient[i][1]+=f2y;
783 gradient[i][2]+=f2z;
784 gradient[i+1][0]+=f3x;
785 gradient[i+1][1]+=f3y;
786 gradient[i+1][2]+=f3z;
787 }
788 }
789
790 #endif
791
792 }
793
794 //printf("ene[3] = %f\n", ene[3]);
795
796 return new_e_pot;
797 }
798
799 /*
800 * Steepest gradient optimization using v=k*(r0-r)^2
801 * k = CA_K, r0 = CA_DIST
802 */
803 void ca_optimize(char *tname, char *iname)
804 {
805 char buf[1000];
806 int i, j, hx, my_iter;
807 real dx, dy, dz, dd, dist, dist2, dist3, ddist, ddist2;
808 real e_pot, new_e_pot, grad, alpha, e_pot1, e_pot2, e_pot3;
809 real adiff[3], bdiff[3];
810 real ff0, ff2, aa, ab, bb, th, tdif, dth, m0, m2;
811 real theta0, deg_th, maxgrad, sum;
812 real f0[3], f2[3];
813 real x, y, z;
814 int numsteps, numsteps2, msteps;
815 int *sec;
816 real **new_c_alpha, **gradient, **init_c_alpha, last_alpha, tmp, last_good_alpha, d_alpha, last_e_pot;
817 atom_type *atom, **c_alpha;
818 res_type *res;
819 FILE *inp, *out;
820 int mnum, init, ok;
821 real alpha1, alpha2, alpha3, a0;
822 real ene1, ene2, ene3, e0;
823 real energies[4];
824 real w1, w2, w3, eps;
825 real gnorm, last_gnorm;
826 int mode, fcnt;
827
828
829 if (_CA_TRAJECTORY) {
830 out = fopen(tname,"w");
831 if (out) fclose(out);
832 }
833
834 if (_VERBOSE) printf("Alpha carbons optimization...\n");
835
836 new_c_alpha = (real**)calloc(sizeof(real*)*(chain_length+1),1);
837 init_c_alpha = (real**)calloc(sizeof(real*)*(chain_length+1),1);
838 for (i=0;i<=chain_length;i++) {
839 new_c_alpha[i] = (real*)calloc(sizeof(real)*3,1);
840 init_c_alpha[i] = (real*)calloc(sizeof(real)*3,1);
841 }
842 gradient = (real**)calloc(sizeof(real*)*(chain_length+1),1);
843 for (i=0;i<=chain_length;i++) {
844 gradient[i] = (real*)calloc(sizeof(real)*3,1);
845 }
846
847 c_alpha = (atom_type**)calloc(sizeof(atom_type*)*(chain_length+1),1);
848
849 i = 0;
850 res = chain->residua;
851 while (res) {
852 atom = res->atoms;
853 while (atom) {
854 if (atom->name[0]=='C' && atom->name[1]=='A') {
855 if (i<chain_length) {
856 c_alpha[i] = atom;
857 i++;
858 break;
859 } else {
860 if (_VERBOSE) printf("WARNING: number of C-alpha atoms exceeds the chain length!\n");
861 break;
862 }
863 }
864 atom = atom->next;
865 }
866 res = res->next;
867 }
868
869 if (i<chain_length) chain_length = i;
870
871 for (i=0; i<chain_length; i++) {
872 init_c_alpha[i][0] = c_alpha[i]->x;
873 init_c_alpha[i][1] = c_alpha[i]->y;
874 init_c_alpha[i][2] = c_alpha[i]->z;
875 }
876
877 if (_CISPRO) {
878 for (i=1; i<chain_length; i++) {
879 dx = c_alpha[i]->x-c_alpha[i-1]->x;
880 dy = c_alpha[i]->y-c_alpha[i-1]->y;
881 dz = c_alpha[i]->z-c_alpha[i-1]->z;
882 dd = sqrt(dx*dx+dy*dy+dz*dz);
883 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)) {
884 if (_VERBOSE) printf("Probable cis-proline found at postion %d\n", c_alpha[i]->res->num);
885 c_alpha[i]->cispro = 1;
886 }
887 }
888 }
889
890 if (_CA_RANDOM) {
891 if (_VERBOSE) printf("Generating random C-alpha coordinates...\n");
892 c_alpha[0]->x = 0.0;
893 c_alpha[0]->y = 0.0;
894 c_alpha[0]->z = 0.0;
895 for (i=1;i<chain_length;i++) {
896 dx = 0.01*(100-rand()%200);
897 dy = 0.01*(100-rand()%200);
898 dz = 0.01*(100-rand()%200);
899 dd = 3.8/sqrt(dx*dx+dy*dy+dz*dz);
900 dx *= dd;
901 dy *= dd;
902 dz *= dd;
903 c_alpha[i]->x = c_alpha[i-1]->x+dx;
904 c_alpha[i]->y = c_alpha[i-1]->y+dy;
905 c_alpha[i]->z = c_alpha[i-1]->z+dz;
906 }
907 }
908
909 if (iname) {
910 inp = fopen(iname,"r");
911 if (inp) {
912 if (_VERBOSE) printf("Reading initial structure %s...\n", iname);
913 i = 0;
914 while (!feof(inp)) {
915 if (fgets(buf,1000,inp)==buf && buf[13]=='C' && buf[14]=='A') {
916 if (i<chain_length) {
917 if (sscanf(&buf[30],"%lf%lf%lf",&x,&y,&z)==3) {
918 c_alpha[i]->x = x;
919 c_alpha[i]->y = y;
920 c_alpha[i]->z = z;
921 i++;
922 }
923 } else {
924 if (_VERBOSE) printf("WARNING: number of ini-file C-alpha atoms exceeds the chain length!\n");
925 break;
926 }
927 }
928 }
929 fclose(inp);
930 } else
931 if (_VERBOSE) printf("WARNING: can't read initial corrdinates %s\n", iname);
932 }
933
934 mnum = 1;
935 mode = 0;
936 init = 0;
937 numsteps=numsteps2=0;
938 last_alpha = 0.0;
939
940
941 if (_VERBOSE) printf("Optimizing alpha carbons...\n");
942
943 eps = 0.5;
944
945 fcnt=0;
946
947 last_gnorm = 1000.;
948
949 do {
950 last_e_pot = e_pot;
951
952 if (_CA_TRAJECTORY) {
953 out = fopen(tname,"a");
954 if (out) {
955 fprintf(out,"MODEL %d\n",mnum++);
956 for (i=0; i<chain_length; i++) {
957 fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n",
958 i+1, "CA ", c_alpha[i]->res->name, ' ', c_alpha[i]->res->num,
959 c_alpha[i]->x, c_alpha[i]->y, c_alpha[i]->z);
960
961 }
962 fprintf(out,"ENDMDL\n");
963 fclose(out);
964 }
965 }
966
967 // calculate gradients
968
969 e_pot=e_pot1=e_pot2=e_pot3=0.;
970
971 for (i=0; i<chain_length; i++)
972 gradient[i][0]=gradient[i][1]=gradient[i][2]=0.;
973
974 e_pot = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, 0.0, energies, true);
975
976 if (_VERBOSE && !init) {
977 printf("Initial energy: bond=%.5lf angle=%.5f restraints=%.5f xvol=%.5f total=%.5f\n", energies[0], energies[2], energies[1], energies[3], e_pot);
978 }
979
980 if (!init) init=1;
981
982 // LINE SEARCH
983
984 alpha1 = -1.0;
985 alpha2 = 0.0;
986 alpha3 = 1.0;
987
988 ene1 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha1, energies, false);
989 ene2 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha2, energies, false);
990 ene3 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha3, energies, false);
991
992 msteps = 0;
993 while (ene2>MIN(ene1,ene3) && msteps<_CA_ITER) {
994 msteps++;
995 alpha1 *= 2.0;
996 alpha3 *= 2.0;
997 ene1 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha1, energies, false);
998 ene3 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, alpha3, energies, false);
999 }
1000
1001 msteps = 0;
1002 do {
1003 if (alpha3-alpha2>alpha2-alpha1) {
1004 a0 = 0.5*(alpha2+alpha3);
1005 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false);
1006 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0-1e-5, energies, false);
1007 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0+1e-5, energies, false);
1008 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false);
1009 if (e0<ene2) {
1010 alpha1 = alpha2;
1011 alpha2 = a0;
1012 ene1 = ene2;
1013 ene2 = e0;
1014 } else {
1015 alpha3 = a0;
1016 ene3 = e0;
1017 }
1018 } else {
1019 a0 = 0.5*(alpha1+alpha2);
1020 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false);
1021 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0-1e-5, energies, false);
1022 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0+1e-5, energies, false);
1023 e0 = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, a0, energies, false);
1024 if (e0<ene2) {
1025 alpha3 = alpha2;
1026 alpha2 = a0;
1027 ene3 = ene2;
1028 ene2 = e0;
1029 } else {
1030 alpha1 = a0;
1031 ene1 = e0;
1032 }
1033 }
1034 msteps++;
1035 } while (alpha3-alpha1>1e-6 && msteps<20);
1036
1037 last_alpha = alpha2;
1038 e_pot = ene2;
1039
1040 for (i=0; i<chain_length; i++) {
1041 c_alpha[i]->x=c_alpha[i]->x+(last_alpha+last_alpha*(rnd()-0.5)*eps)*gradient[i][0];
1042 c_alpha[i]->y=c_alpha[i]->y+(last_alpha+last_alpha*(rnd()-0.5)*eps)*gradient[i][1];
1043 c_alpha[i]->z=c_alpha[i]->z+(last_alpha+last_alpha*(rnd()-0.5)*eps)*gradient[i][2];
1044 }
1045
1046 e_pot = calc_ca_energy(c_alpha, new_c_alpha, init_c_alpha, gradient, 0.0, energies, false);
1047
1048 eps *= 0.75;
1049 if (eps<1e-3) eps=0.0;
1050
1051 numsteps++;
1052
1053 gnorm = 0.0;
1054
1055 for (i=0; i<chain_length; i++) {
1056 gnorm += gradient[i][0]*gradient[i][0] + gradient[i][1]*gradient[i][1] + gradient[i][2]*gradient[i][2];
1057 }
1058
1059 gnorm = sqrt(gnorm/(double)chain_length);
1060
1061 if (last_gnorm-gnorm<1e-3) fcnt++;
1062
1063 last_gnorm = gnorm;
1064
1065 } while ( (fcnt<3) && (gnorm>0.01) && (numsteps<_CA_ITER));
1066
1067
1068 if (_VERBOSE) {
1069 for (i=0; i<chain_length; i++) {
1070
1071 #ifdef CALC_C_ALPHA
1072 if (i>0) {
1073 dx=c_alpha[i]->x-c_alpha[i-1]->x;
1074 dy=c_alpha[i]->y-c_alpha[i-1]->y;
1075 dz=c_alpha[i]->z-c_alpha[i-1]->z;
1076 dist=sqrt(dx*dx+dy*dy+dz*dz);
1077 if (c_alpha[i]->cispro) {
1078 ddist=CA_DIST_CISPRO-dist;
1079 if (fabs(ddist)<CA_DIST_CISPRO_TOL) ddist=0.0;
1080 } else {
1081 ddist=CA_DIST-dist;
1082 if (fabs(ddist)<CA_DIST_TOL) ddist=0.0;
1083 }
1084 ddist2=ddist*ddist;
1085 if (fabs(ddist)>=CA_DIST_TOL) printf("WARNING: distance %d = %.3lf A\n", i, dist);
1086 }
1087 #endif
1088 }
1089
1090 for (i=0; i<chain_length; i++) {
1091 #ifdef CALC_C_ALPHA_ANGLES
1092 if (i>0 && i<chain_length-1) {
1093 aa=ab=bb=0.0;
1094 adiff[0]=c_alpha[i-1]->x-c_alpha[i]->x;
1095 bdiff[0]=c_alpha[i+1]->x-c_alpha[i]->x;
1096 aa+=adiff[0]*adiff[0];
1097 ab+=adiff[0]*bdiff[0];
1098 bb+=bdiff[0]*bdiff[0];
1099 adiff[1]=c_alpha[i-1]->y-c_alpha[i]->y;
1100 bdiff[1]=c_alpha[i+1]->y-c_alpha[i]->y;
1101 aa+=adiff[1]*adiff[1];
1102 ab+=adiff[1]*bdiff[1];
1103 bb+=bdiff[1]*bdiff[1];
1104 adiff[2]=c_alpha[i-1]->z-c_alpha[i]->z;
1105 bdiff[2]=c_alpha[i+1]->z-c_alpha[i]->z;
1106 aa+=adiff[2]*adiff[2];
1107 ab+=adiff[2]*bdiff[2];
1108 bb+=bdiff[2]*bdiff[2];
1109
1110 th=ab/sqrt(aa*bb);
1111 if (th<-1.0) th=-1.0;
1112 if (th>1.0) th=1.0;
1113 th=acos(th);
1114 deg_th=RADDEG*th;
1115 if (deg_th>150.) theta0=DEGRAD*150.; else
1116 if (deg_th<75.) theta0=DEGRAD*75.; else
1117 theta0=th;
1118 if (fabs(deg_th-RADDEG*theta0)>1.0) printf("WARNING: angle %d = %.3lf degrees\n", i, deg_th);
1119 }
1120 #endif
1121 }
1122 }
1123
1124 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);
1125
1126 if (_CA_TRAJECTORY) {
1127 out = fopen(tname,"a");
1128 if (out) {
1129 fprintf(out,"END\n");
1130 }
1131 }
1132
1133 for (i=0;i<chain_length+1;i++) {
1134 free(init_c_alpha[i]);
1135 free(new_c_alpha[i]);
1136 free(gradient[i]);
1137 }
1138 free(new_c_alpha);
1139 free(gradient);
1140 free(c_alpha);
1141 free(init_c_alpha);
1142 }
1143
1144 void center_chain(mol_type *mol)
1145 {
1146 real cx, cy, cz;
1147 int natom;
1148 res_type *res;
1149 atom_type *atom;
1150
1151 cx = cy = cz = 0.0;
1152 natom = 0;
1153
1154 res = mol->residua;
1155 while (res) {
1156 atom = res->atoms;
1157 while (atom) {
1158 cx += atom->x;
1159 cy += atom->y;
1160 cz += atom->z;
1161 natom++;
1162 atom=atom->next;
1163 }
1164 res = res->next;
1165 }
1166
1167 cx /= (real)natom;
1168 cy /= (real)natom;
1169 cz /= (real)natom;
1170
1171 if (_VERBOSE) printf("Molecule center: %8.3f %8.3f %8.3f -> 0.000 0.000 0.000\n", cx, cy, cz);
1172
1173 res = mol->residua;
1174 while (res) {
1175 atom = res->atoms;
1176 while (atom) {
1177 atom->x -= cx;
1178 atom->y -= cy;
1179 atom->z -= cz;
1180 natom++;
1181 atom=atom->next;
1182 }
1183 res = res->next;
1184 }
1185
1186 }
1187
1188 void write_pdb(char *name, mol_type *mol)
1189 {
1190 FILE *out;
1191 res_type *res;
1192 atom_type *atom;
1193 int anum;
1194
1195 out = fopen(name,"w");
1196 if (!out) {
1197 if (_VERBOSE) printf("Can't open output file!\n");
1198 return;
1199 }
1200 fprintf(out,"REMARK 999 REBUILT BY PULCHRA V.%.2f\n", PULCHRA_VERSION);
1201 anum=1;
1202 res = mol->residua;
1203 while (res) {
1204 if (res->protein) {
1205 if (!_BB_REARRANGE) {
1206 atom = res->atoms;
1207 while (atom) {
1208 if (!(atom->name[0]=='D' && atom->name[1]=='U') &&
1209 !(atom->name[0]=='S' && atom->name[1]=='C') &&
1210 !(atom->name[0]=='C' && atom->name[1]=='M') &&
1211 !(atom->name[0]=='H' && !_REBUILD_H))
1212 fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n",
1213 anum++, atom->name, res->name, ' ', res->num,
1214 atom->x, atom->y, atom->z);
1215 atom=atom->next;
1216 }
1217 } else {
1218 atom = res->atoms;
1219 while (atom) {
1220 if (!(atom->name[0]=='D' && atom->name[1]=='U') &&
1221 !(atom->name[0]=='S' && atom->name[1]=='C') &&
1222 !(atom->name[0]=='C' && atom->name[1]==' ') &&
1223 !(atom->name[0]=='O' && atom->name[1]==' ') &&
1224 !(atom->name[0]=='C' && atom->name[1]=='M') &&
1225 !(atom->name[0]=='H' && !_REBUILD_H))
1226 fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n",
1227 anum++, atom->name, res->name, ' ', res->num,
1228 atom->x, atom->y, atom->z);
1229 atom=atom->next;
1230 }
1231 atom = res->atoms;
1232 while (atom) {
1233 if (((atom->name[0]=='C' && atom->name[1]==' ') ||
1234 (atom->name[0]=='O' && atom->name[1]==' ')) &&
1235 !(atom->name[0]=='H' && !_REBUILD_H))
1236 fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n",
1237 anum++, atom->name, res->name, ' ', res->num,
1238 atom->x, atom->y, atom->z);
1239 atom=atom->next;
1240 }
1241 }
1242 }
1243 res = res->next;
1244 }
1245 fprintf(out,"TER\nEND\n");
1246 fclose(out);
1247 }
1248
1249
1250 void write_pdb_sg(char *name, mol_type *mol)
1251 {
1252 FILE *out;
1253 res_type *res;
1254 atom_type *atom;
1255 int anum;
1256
1257 out = fopen(name,"w");
1258 if (!out) {
1259 if (_VERBOSE) printf("Can't open output file!\n");
1260 return;
1261 }
1262 fprintf(out,"REMARK 999 REBUILT BY PULCHRA V.%.2f\n", PULCHRA_VERSION);
1263 anum=1;
1264 res = mol->residua;
1265 while (res) {
1266 if (res->protein) {
1267 atom = res->atoms;
1268 while (atom) {
1269 if ((atom->name[0]=='C' && atom->name[1]=='A'))
1270 fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n",
1271 anum++, atom->name, res->name, ' ', res->num,
1272 atom->x, atom->y, atom->z);
1273 atom=atom->next;
1274 }
1275 fprintf(out, "ATOM %5d %-3s %3s %c%4d %8.3f%8.3f%8.3f\n",
1276 anum++, "CM ", res->name, ' ', res->num,
1277 res->cmx, res->cmy, res->cmz);
1278 }
1279 res = res->next;
1280 }
1281 fprintf(out,"TER\nEND\n");
1282 fclose(out);
1283 }
1284
1285 real calc_distance(real x1, real y1, real z1,
1286 real x2, real y2, real z2)
1287 {
1288 real dx,dy,dz;
1289 real dist2;
1290
1291 dx = (x1) - (x2);
1292 dy = (y1) - (y2);
1293 dz = (z1) - (z2);
1294 if (dx || dy || dz ) {
1295 dist2 = dx*dx+dy*dy+dz*dz;
1296 return (sqrt(dist2));
1297 } else
1298 return 0.0;
1299 }
1300
1301 real calc_r14(real x1, real y1, real z1,
1302 real x2, real y2, real z2,
1303 real x3, real y3, real z3,
1304 real x4, real y4, real z4)
1305 {
1306 real r, dx, dy, dz;
1307 real vx1, vy1, vz1, vx2, vy2, vz2, vx3, vy3, vz3;
1308 real hand;
1309
1310 dx = x4-x1;
1311 dy = y4-y1;
1312 dz = z4-z1;
1313
1314 r = sqrt(dx*dx+dy*dy+dz*dz);
1315
1316 vx1=x2-x1;
1317 vy1=y2-y1;
1318 vz1=z2-z1;
1319 vx2=x3-x2;
1320 vy2=y3-y2;
1321 vz2=z3-z2;
1322 vx3=x4-x3;
1323 vy3=y4-y3;
1324 vz3=z4-z3;
1325
1326 hand = (vy1*vz2-vy2*vz1)*vx3+
1327 (vz1*vx2-vz2*vx1)*vy3+
1328 (vx1*vy2-vx2*vy1)*vz3;
1329
1330 if (hand<0) r=-r;
1331
1332 return r;
1333 }
1334
1335 real superimpose2(real **coords1, real **coords2, int npoints, real **tpoints, int ntpoints)
1336 {
1337 real mat_s[3][3], mat_a[3][3], mat_b[3][3], mat_g[3][3];
1338 real mat_u[3][3], tmp_mat[3][3];
1339 real val, d, alpha, beta, gamma, x, y, z;
1340 real cx1, cy1, cz1, cx2, cy2, cz2, tmpx, tmpy, tmpz;
1341 int i, j, k, n;
1342
1343 cx1=cy1=cz1=cx2=cy2=cz2=0.;
1344
1345 for (i=0; i<npoints; i++) {
1346 cx1+=coords1[i][0];
1347 cy1+=coords1[i][1];
1348 cz1+=coords1[i][2];
1349 cx2+=coords2[i][0];
1350 cy2+=coords2[i][1];
1351 cz2+=coords2[i][2];
1352 }
1353
1354 cx1/=(real)npoints;
1355 cy1/=(real)npoints;
1356 cz1/=(real)npoints;
1357
1358 cx2/=(real)npoints;
1359 cy2/=(real)npoints;
1360 cz2/=(real)npoints;
1361
1362 for (i=0; i<npoints; i++) {
1363 coords1[i][0]-=cx1;
1364 coords1[i][1]-=cy1;
1365 coords1[i][2]-=cz1;
1366 coords2[i][0]-=cx2;
1367 coords2[i][1]-=cy2;
1368 coords2[i][2]-=cz2;
1369 }
1370
1371 for (i=0; i<ntpoints; i++) {
1372 tpoints[i][0]-=cx2;
1373 tpoints[i][1]-=cy2;
1374 tpoints[i][2]-=cz2;
1375 }
1376
1377 for (i=0; i<3; i++)
1378 for (j=0; j<3; j++) {
1379 if (i==j)
1380 mat_s[i][j]=mat_a[i][j]=mat_b[i][j]=mat_g[i][j]=1.0;
1381 else
1382 mat_s[i][j]=mat_a[i][j]=mat_b[i][j]=mat_g[i][j]=0.0;
1383 mat_u[i][j]=0.;
1384 }
1385
1386 for (n=0; n<npoints; n++) {
1387 mat_u[0][0]+=coords1[n][0]*coords2[n][0];
1388 mat_u[0][1]+=coords1[n][0]*coords2[n][1];
1389 mat_u[0][2]+=coords1[n][0]*coords2[n][2];
1390 mat_u[1][0]+=coords1[n][1]*coords2[n][0];
1391 mat_u[1][1]+=coords1[n][1]*coords2[n][1];
1392 mat_u[1][2]+=coords1[n][1]*coords2[n][2];
1393 mat_u[2][0]+=coords1[n][2]*coords2[n][0];
1394 mat_u[2][1]+=coords1[n][2]*coords2[n][1];
1395 mat_u[2][2]+=coords1[n][2]*coords2[n][2];
1396 }
1397
1398 for (i=0; i<3; i++)
1399 for (j=0; j<3; j++)
1400 tmp_mat[i][j]=0.;
1401
1402 do {
1403 d=mat_u[2][1]-mat_u[1][2];
1404 if (d==0) alpha=0; else alpha=atan(d/(mat_u[1][1]+mat_u[2][2]));
1405 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;
1406 mat_a[1][1]=mat_a[2][2]=cos(alpha);
1407 mat_a[2][1]=sin(alpha);
1408 mat_a[1][2]=-mat_a[2][1];
1409 for (i=0; i<3; i++)
1410 for (j=0; j<3; j++)
1411 for (k=0; k<3; k++)
1412 tmp_mat[i][j]+=mat_u[i][k]*mat_a[j][k];
1413 for (i=0; i<3; i++)
1414 for (j=0; j<3; j++) {
1415 mat_u[i][j]=tmp_mat[i][j];
1416 tmp_mat[i][j]=0.;
1417 }
1418 for (i=0; i<3; i++)
1419 for (j=0; j<3; j++)
1420 for (k=0; k<3; k++)
1421 tmp_mat[i][j]+=mat_a[i][k]*mat_s[k][j];
1422 for (i=0; i<3; i++)
1423 for (j=0; j<3; j++) {
1424 mat_s[i][j]=tmp_mat[i][j];
1425 tmp_mat[i][j]=0.;
1426 }
1427 d=mat_u[0][2]-mat_u[2][0];
1428 if (d==0) beta=0; else beta=atan(d/(mat_u[0][0]+mat_u[2][2]));
1429 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;
1430 mat_b[0][0]=mat_b[2][2]=cos(beta);
1431 mat_b[0][2]=sin(beta);
1432 mat_b[2][0]=-mat_b[0][2];
1433 for (i=0; i<3; i++)
1434 for (j=0; j<3; j++)
1435 for (k=0; k<3; k++)
1436 tmp_mat[i][j]+=mat_u[i][k]*mat_b[j][k];
1437 for (i=0; i<3; i++)
1438 for (j=0; j<3; j++) {
1439 mat_u[i][j]=tmp_mat[i][j];
1440 tmp_mat[i][j]=0.;
1441 }
1442 for (i=0; i<3; i++)
1443 for (j=0; j<3; j++)
1444 for (k=0; k<3; k++)
1445 tmp_mat[i][j]+=mat_b[i][k]*mat_s[k][j];
1446 for (i=0; i<3; i++)
1447 for (j=0; j<3; j++) {
1448 mat_s[i][j]=tmp_mat[i][j];
1449 tmp_mat[i][j]=0.;
1450 }
1451 d=mat_u[1][0]-mat_u[0][1];
1452 if (d==0) gamma=0; else gamma=atan(d/(mat_u[0][0]+mat_u[1][1]));
1453 if (cos(gamma)*(mat_u[0][0]+mat_u[1][1])+sin(gamma)*(mat_u[1][0]-mat_u[0][1])<0.0)
1454 gamma+=M_PI;
1455 mat_g[0][0]=mat_g[1][1]=cos(gamma);
1456 mat_g[1][0]=sin(gamma);
1457 mat_g[0][1]=-mat_g[1][0];
1458 for (i=0; i<3; i++)
1459 for (j=0; j<3; j++)
1460 for (k=0; k<3; k++)
1461 tmp_mat[i][j]+=mat_u[i][k]*mat_g[j][k];
1462 for (i=0; i<3; i++)
1463 for (j=0; j<3; j++) {
1464 mat_u[i][j]=tmp_mat[i][j];
1465 tmp_mat[i][j]=0.;
1466 }
1467 for (i=0; i<3; i++)
1468 for (j=0; j<3; j++)
1469 for (k=0; k<3; k++)
1470 tmp_mat[i][j]+=mat_g[i][k]*mat_s[k][j];
1471 for (i=0; i<3; i++)
1472 for (j=0; j<3; j++) {
1473 mat_s[i][j]=tmp_mat[i][j];
1474 tmp_mat[i][j]=0.;
1475 }
1476 val=fabs(alpha)+fabs(beta)+fabs(gamma);
1477 } while (val>0.001);
1478
1479 val=0.;
1480 for (i=0; i<npoints; i++) {
1481 x=coords2[i][0];
1482 y=coords2[i][1];
1483 z=coords2[i][2];
1484 tmpx=x*mat_s[0][0]+y*mat_s[0][1]+z*mat_s[0][2];
1485 tmpy=x*mat_s[1][0]+y*mat_s[1][1]+z*mat_s[1][2];
1486 tmpz=x*mat_s[2][0]+y*mat_s[2][1]+z*mat_s[2][2];
1487 x=coords1[i][0]-tmpx;
1488 y=coords1[i][1]-tmpy;
1489 z=coords1[i][2]-tmpz;
1490 val+=x*x+y*y+z*z;
1491 }
1492
1493 for (i=0; i<ntpoints; i++) {
1494 x=tpoints[i][0];
1495 y=tpoints[i][1];
1496 z=tpoints[i][2];
1497 tpoints[i][0]=x*mat_s[0][0]+y*mat_s[0][1]+z*mat_s[0][2];
1498 tpoints[i][1]=x*mat_s[1][0]+y*mat_s[1][1]+z*mat_s[1][2];
1499 tpoints[i][2]=x*mat_s[2][0]+y*mat_s[2][1]+z*mat_s[2][2];
1500 }
1501
1502 for (i=0; i<npoints; i++) {
1503 coords1[i][0]+=cx1;
1504 coords1[i][1]+=cy1;
1505 coords1[i][2]+=cz1;
1506 coords2[i][0]+=cx2;
1507 coords2[i][1]+=cy2;
1508 coords2[i][2]+=cz2;
1509 }
1510
1511 for (i=0; i<ntpoints; i++) {
1512 tpoints[i][0]+=cx1;
1513 tpoints[i][1]+=cy1;
1514 tpoints[i][2]+=cz1;
1515 }
1516
1517 return sqrt(val/(real)npoints);
1518 }
1519
1520
1521 atom_type *find_atom(res_type *res, char *aname)
1522 {
1523 atom_type *atom;
1524
1525 atom = res->atoms;
1526 while (atom) {
1527 if (atom->name[0]==aname[0] && atom->name[1]==aname[1] && atom->name[2]==aname[2]) {
1528 return atom;
1529 break;
1530 }
1531 atom = atom->next;
1532 }
1533
1534 return NULL;
1535 }
1536
1537
1538 void add_replace(res_type *res, char *aname, real x, real y, real z, int flags)
1539 {
1540 atom_type *atom, *newatom;
1541
1542 atom = res->atoms;
1543 while (atom) {
1544 if (atom->name[0]==aname[0] && atom->name[1]==aname[1] && atom->name[2]==aname[2]) {
1545 atom->x = x; atom->y = y; atom->z = z;
1546 atom->flag |= flags;
1547 break;
1548 }
1549 atom = atom->next;
1550 }
1551
1552 if (!atom) {
1553 newatom = (atom_type*)calloc(sizeof(atom_type),1);
1554 newatom->x = x;
1555 newatom->y = y;
1556 newatom->z = z;
1557 newatom->flag |= flags;
1558 newatom->res = res;
1559 newatom->name = (char*)calloc(4,1);
1560 strcpy(newatom->name,aname);
1561
1562 atom = res->atoms;
1563 while (atom) {
1564 if (atom->name[0]=='C' && atom->name[1]=='A')
1565 break;
1566 atom = atom->next;
1567 }
1568 if (aname[0]=='N' && aname[1]==' ') {
1569 newatom->next = res->atoms;
1570 res->atoms = newatom;
1571 } else {
1572 while (atom->next) atom=atom->next;
1573 atom->next = newatom;
1574 }
1575 }
1576 }
1577
1578
1579 int **RBINS;
1580 real **X_COORDS, **C_ALPHA;
1581
1582 #ifdef COMPILE_BB
1583
1584 void rebuild_backbone(void)
1585 {
1586
1587 res_type *res, *prevres;
1588 atom_type *atom;
1589 real **cacoords, **tmpcoords, **tmpstat;
1590 real x1, y1, z1;
1591 real x2, y2, z2;
1592 real x3, y3, z3;
1593 real x4, y4, z4;
1594 real r13_1, r13_2, r14;
1595 real besthit, hit;
1596 int bestpos;
1597 int i, j, k, l, m, bin13_1, bin13_2, bin14, found, pro;
1598 int b13_1, b13_2, b14;
1599 real rmsd, total, maxrms;
1600 FILE *debug, *out;
1601
1602 if (_VERBOSE) printf("Rebuilding backbone...\n");
1603
1604 RBINS = (int**)calloc(sizeof(int*)*(chain_length+1),1);
1605 for (i=0;i<chain_length+1;i++)
1606 RBINS[i] = (int*)calloc(sizeof(int)*3,1);
1607
1608 X_COORDS = (real**)calloc(sizeof(real*)*(chain_length+10),1);
1609 for (i=0;i<chain_length+10;i++)
1610 X_COORDS[i] = (real*)calloc(sizeof(real)*3,1);
1611
1612 i = 5;
1613
1614 res = chain->residua;
1615 while (res) {
1616 atom = res->atoms;
1617 while (atom) {
1618 if (atom->name[0]=='C' && atom->name[1]=='A') {
1619 X_COORDS[i][0] = atom->x;
1620 X_COORDS[i][1] = atom->y;
1621 X_COORDS[i][2] = atom->z;
1622 i++;
1623 }
1624 atom = atom->next;
1625 }
1626 res = res->next;
1627 }
1628
1629 cacoords = (real**)calloc(sizeof(real*)*(8),1);
1630 tmpcoords = (real**)calloc(sizeof(real*)*(8),1);
1631 tmpstat = (real**)calloc(sizeof(real*)*(8),1);
1632 for (i=0;i<8;i++) {
1633 cacoords[i] = (real*)calloc(sizeof(real)*3,1);;
1634 tmpcoords[i] = (real*)calloc(sizeof(real)*3,1);;
1635 tmpstat[i] = (real*)calloc(sizeof(real)*3,1);;
1636 }
1637
1638 C_ALPHA = &X_COORDS[5];
1639
1640 // rebuild ends...
1641
1642 for (i=0,j=0;i<5;i++,j++)
1643 for (k=0;k<3;k++)
1644 tmpcoords[j][k] = C_ALPHA[i][k];
1645 for (i=2,j=0;i<5;i++,j++)
1646 for (k=0;k<3;k++)
1647 cacoords[j][k] = C_ALPHA[i][k];
1648 for (i=0,j=0;i<3;i++,j++)
1649 for (k=0;k<3;k++)
1650 tmpstat[j][k] = C_ALPHA[i][k];
1651
1652 superimpose2(tmpstat,cacoords,3,tmpcoords,5);
1653
1654 for (i=-2,j=0;i<0;i++,j++)
1655 for (k=0;k<3;k++)
1656 C_ALPHA[i][k] = tmpcoords[j][k];
1657
1658 for (i=chain_length-5,j=0;i<chain_length;i++,j++)
1659 for (k=0;k<3;k++)
1660 tmpcoords[j][k] = C_ALPHA[i][k];
1661 for (i=chain_length-5,j=0;i<chain_length-2;i++,j++)
1662 for (k=0;k<3;k++)
1663 cacoords[j][k] = C_ALPHA[i][k];
1664 for (i=chain_length-3,j=0;i<chain_length;i++,j++)
1665 for (k=0;k<3;k++)
1666 tmpstat[j][k] = C_ALPHA[i][k];
1667
1668 superimpose2(tmpstat,cacoords,3,tmpcoords,5);
1669
1670 for (i=chain_length-3,j=0;i<chain_length;i++,j++)
1671 for (k=0;k<3;k++)
1672 C_ALPHA[i+3][k] = tmpcoords[j+3][k];
1673
1674
1675 prevres = NULL;
1676 res = chain->residua;
1677
1678
1679 total = maxrms = 0.0;
1680
1681 for (i=0;i<chain_length+1;i++) {
1682 x1 = C_ALPHA[i-2][0];
1683 y1 = C_ALPHA[i-2][1];
1684 z1 = C_ALPHA[i-2][2];
1685
1686 x2 = C_ALPHA[i-1][0];
1687 y2 = C_ALPHA[i-1][1];
1688 z2 = C_ALPHA[i-1][2];
1689
1690 x3 = C_ALPHA[i][0];
1691 y3 = C_ALPHA[i][1];
1692 z3 = C_ALPHA[i][2];
1693
1694 x4 = C_ALPHA[i+1][0];
1695 y4 = C_ALPHA[i+1][1];
1696 z4 = C_ALPHA[i+1][2];
1697
1698 r13_1 = calc_distance(x1, y1, z1, x3, y3, z3);
1699 r13_2 = calc_distance(x2, y2, z2, x4, y4, z4);
1700 r14 = calc_r14(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4);
1701
1702 bin13_1 = (int)((r13_1-4.6)/0.3);
1703 bin13_2 = (int)((r13_2-4.6)/0.3);
1704 bin14 = (int)((r14+11.)/0.3);
1705
1706 if (bin13_1<0) bin13_1=0;
1707 if (bin13_2<0) bin13_2=0;
1708 if (bin14<0) bin14=0;
1709 if (bin13_1>9) bin13_1=9;
1710 if (bin13_2>9) bin13_2=9;
1711 if (bin14>73) bin14=73;
1712
1713 RBINS[i][0] = bin13_1;
1714 RBINS[i][1] = bin13_2;
1715 RBINS[i][2] = bin14;
1716
1717 cacoords[0][0] = x1;
1718 cacoords[0][1] = y1;
1719 cacoords[0][2] = z1;
1720
1721 cacoords[1][0] = x2;
1722 cacoords[1][1] = y2;
1723 cacoords[1][2] = z2;
1724
1725 cacoords[2][0] = x3;
1726 cacoords[2][1] = y3;
1727 cacoords[2][2] = z3;
1728
1729 cacoords[3][0] = x4;
1730 cacoords[3][1] = y4;
1731 cacoords[3][2] = z4;
1732
1733 pro = 0;
1734
1735 if (prevres && !strncmp(prevres->name,"PRO",3)) {
1736 j=0;
1737 besthit=1000.;
1738 bestpos=0;
1739 do {
1740 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);
1741 if (hit<besthit) {
1742 besthit=hit;
1743 bestpos=j;
1744 }
1745 j++;
1746 } while (nco_stat_pro[j].bins[0]>=0 && hit>1e-3);
1747 for (j=0;j<4;j++) {
1748 for (k=0;k<3;k++) {
1749 tmpstat[j][k] = nco_stat_pro[bestpos].data[j][k];
1750 }
1751 }
1752 for (j=0;j<8;j++) {
1753 for (k=0;k<3;k++) {
1754 tmpcoords[j][k] = nco_stat_pro[bestpos].data[j][k];
1755 }
1756 }
1757 } else {
1758 j=0;
1759 besthit=1000.;
1760 bestpos=0;
1761 do {
1762 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);
1763 if (hit<besthit) {
1764 besthit=hit;
1765 bestpos=j;
1766 }
1767 j++;
1768 } while (nco_stat[j].bins[0]>=0 && hit>1e-3);
1769 for (j=0;j<4;j++) {
1770 for (k=0;k<3;k++) {
1771 tmpstat[j][k] = nco_stat[bestpos].data[j][k];
1772 }
1773 }
1774 for (j=0;j<8;j++) {
1775 for (k=0;k<3;k++) {
1776 tmpcoords[j][k] = nco_stat[bestpos].data[j][k];
1777 }
1778 }
1779 }
1780
1781 rmsd=superimpose2(cacoords, tmpstat, 4, tmpcoords, 8);
1782
1783 total += rmsd;
1784 if (rmsd>maxrms) maxrms=rmsd;
1785
1786 // add-or-replace
1787
1788 if (prevres) {
1789 add_replace(prevres, "C ", tmpcoords[4][0], tmpcoords[4][1], tmpcoords[4][2], FLAG_BACKBONE);
1790 add_replace(prevres, "O ", tmpcoords[5][0], tmpcoords[5][1], tmpcoords[5][2], FLAG_BACKBONE);
1791 }
1792
1793 if (res) {
1794 add_replace(res, "N ", tmpcoords[6][0], tmpcoords[6][1], tmpcoords[6][2], FLAG_BACKBONE);
1795 }
1796
1797 prevres = res;
1798 if (res)
1799 res = res->next;
1800 }
1801
1802 if (_VERBOSE) printf("Backbone rebuilding deviation: average = %.3f, max = %.3f\n", total/(real)chain_length, maxrms);
1803 }
1804
1805 #endif
1806
1807
1808 #ifdef COMPILE_ROT
1809
1810 typedef struct _rot_struct {
1811 int r13_1, r13_2, r14;
1812 int nc;
1813 real ***coords;
1814 struct _rot_struct *next;
1815 } rot_struct;
1816
1817 rot_struct *rotamers[20];
1818
1819 /* this is obsolete in a standalone version of PULCHRA */
1820 void read_rotamers(void)
1821 {
1822 FILE *inp;
1823 char buf[1000];
1824 char dum[100];
1825 int aa, i, j, k, l, n;
1826 rot_struct *new_rot, *last_rot;
1827 real x, y, z;
1828
1829 if (_VERBOSE) printf("Reading rotamer library...\n");
1830
1831 inp = fopen("NEWROT","r");
1832 last_rot=NULL;
1833 while (!feof(inp)) {
1834 if (fgets(buf,1000,inp)==buf) {
1835 if (buf[0]=='A') {
1836 sscanf(buf,"%s %d", dum, &aa);
1837 if (last_rot) last_rot->next = NULL;
1838 last_rot = NULL;
1839 if (fgets(buf,1000,inp)!=buf) break;
1840 }
1841 // printf("aa: %d\n", aa);
1842 if (aa==20) break;
1843 sscanf(buf,"%d %d %d %s %d", &i, &j, &k, dum, &l);
1844 new_rot = (rot_struct*)calloc(sizeof(rot_struct),1);
1845 // printf("%d %d %d nc: %d\n", i, j, k, l);
1846 new_rot->r13_1 = i;
1847 new_rot->r13_2 = j;
1848 new_rot->r14 = k;
1849 new_rot->nc = l;
1850 new_rot->next = NULL;
1851 new_rot->coords = (real***)calloc(sizeof(real**)*l,1);
1852 for (i=0;i<l;i++) {
1853 new_rot->coords[i]=(real**)calloc(sizeof(real*)*(nheavy[aa]+1),1);
1854 for (j=0;j<(nheavy[aa]+1);j++) {
1855 new_rot->coords[i][j]=(real*)calloc(sizeof(real)*3,1);
1856 }
1857 }
1858 for (i=0;i<l;i++) {
1859 fgets(buf,1000,inp);
1860 for (j=0;j<(nheavy[aa]+1);j++) {
1861 fgets(buf,1000,inp);
1862 sscanf(buf,"%lf%lf%lf",&x, &y, &z);
1863 new_rot->coords[i][j][0]=x;
1864 new_rot->coords[i][j][1]=y;
1865 new_rot->coords[i][j][2]=z;
1866 }
1867 if (last_rot) {
1868 last_rot->next = new_rot;
1869 } else {
1870 rotamers[aa] = new_rot;
1871 }
1872 last_rot = new_rot;
1873 }
1874 }
1875 }
1876 fclose(inp);
1877 }
1878
1879
1880 void cross(real *v1, real *v2, real *v3)
1881 {
1882 v3[0] = v1[1]*v2[2]-v1[2]*v2[1];
1883 v3[1] = v1[2]*v2[0]-v1[0]*v2[2];
1884 v3[2] = v1[0]*v2[1]-v1[1]*v2[0];
1885 }
1886
1887 void norm(real *v)
1888 {
1889 real d;
1890
1891 d = sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]);
1892 v[0] /= d;
1893 v[1] /= d;
1894 v[2] /= d;
1895 }
1896
1897
1898 int check_xvol(res_type *res)
1899 {
1900 res_type *res2;
1901 atom_type *atom1, *atom2;
1902 real dx, dy, dz, dd;
1903
1904 res2 = chain->residua;
1905
1906 while (res2) {
1907 atom2 = res2->atoms;
1908 if (res!=res2) {
1909 while (atom2) {
1910 atom1 = res->atoms;
1911 while (atom1) {
1912 if (atom1->flag & FLAG_SIDECHAIN) {
1913 dx = atom1->x-atom2->x;
1914 dy = atom1->y-atom2->y;
1915 dz = atom1->z-atom2->z;
1916 dd = dx*dx+dy*dy+dz*dz;
1917 if (dd<(1.7*1.7)) {
1918 return 1;
1919 }
1920 }
1921 atom1=atom1->next;
1922 }
1923 atom2=atom2->next;
1924 }
1925 }
1926 res2=res2->next;
1927 }
1928
1929 return 0;
1930 }
1931
1932
1933 real ***SORTED_ROTAMERS;
1934
1935 void rebuild_sidechains(void)
1936 {
1937 FILE *out;
1938 res_type *res, *prevres, *testres;
1939 atom_type *atom, *atom1, *atom2;
1940 real **cacoords, **tmpcoords, **tmpstat;
1941 real x1, y1, z1;
1942 real x2, y2, z2;
1943 real x3, y3, z3;
1944 real x4, y4, z4;
1945 real x5, y5, z5;
1946 real r14, r13_1, r13_2;
1947 real dx, dy, dz, dd;
1948 real hit, besthit;
1949 int exvol, bestpos;
1950 int i, j, k, l, m, bin13_1, bin13_2, bin14;
1951 real rmsd, total;
1952 real v1[3], v2a[3], v2b[3], v2[3], v3[3];
1953 int nsc, nca;
1954 real cax, cay, caz;
1955 real **lsys, **vv, **sc;
1956 char scn[12][4];
1957 rot_struct *rot;
1958 int ok, last_a, last_b, last_c, last_d, jpos;
1959 int jx, jy, jz, jxi, jyi, jzi, b13_1, b13_2, b14, jm;
1960 int crot, bestrot, minexvol, totexvol, rtried, pos, cpos;
1961 real cmx, cmy, cmz, ddx, ddy, ddz, ddd, bestdd;
1962 real sort_rot[100][2];
1963
1964 if (_VERBOSE) printf("Rebuilding side chains...\n");
1965
1966 lsys = (real**)calloc(sizeof(real*)*3,1);
1967 vv = (real**)calloc(sizeof(real*)*3,1);
1968 sc = (real**)calloc(sizeof(real*)*12,1);
1969 for (i=0;i<12;i++)
1970 sc[i] = (real*)calloc(sizeof(real)*3,1);
1971 for (i=0;i<3;i++) {
1972 lsys[i] = (real*)calloc(sizeof(real)*3,1);
1973 vv[i] = (real*)calloc(sizeof(real)*3,1);
1974 }
1975
1976 SORTED_ROTAMERS = (real***)calloc(sizeof(real**)*(chain_length+1),1);
1977 for (i=0;i<chain_length+1;i++) {
1978 SORTED_ROTAMERS[i] = (real**)calloc(sizeof(real*)*10,1);
1979 for (j=0;j<10;j++) {
1980 SORTED_ROTAMERS[i][j] = (real*)calloc(sizeof(real)*2,1);
1981 }
1982 }
1983
1984 prevres = NULL;
1985 res = chain->residua;
1986 totexvol = 0;
1987
1988 for (i=0;i<chain_length;i++) {
1989 if (!strncmp(res->name,"GLY",3) || !res->protein) {
1990 if (res->next) res = res->next;
1991 continue;
1992 }
1993
1994 x1 = C_ALPHA[i-2][0];
1995 y1 = C_ALPHA[i-2][1];
1996 z1 = C_ALPHA[i-2][2];
1997 x2 = C_ALPHA[i-1][0];
1998 y2 = C_ALPHA[i-1][1];
1999 z2 = C_ALPHA[i-1][2];
2000 x3 = C_ALPHA[i][0];
2001 y3 = C_ALPHA[i][1];
2002 z3 = C_ALPHA[i][2];
2003 x4 = C_ALPHA[i+1][0];
2004 y4 = C_ALPHA[i+1][1];
2005 z4 = C_ALPHA[i+1][2];
2006
2007 bin13_1 = RBINS[i][0];
2008 bin13_2 = RBINS[i][1];
2009 bin14 = RBINS[i][2];
2010
2011 v1[0] = x4-x2;
2012 v1[1] = y4-y2;
2013 v1[2] = z4-z2;
2014
2015 v2a[0] = x4-x3;
2016 v2a[1] = y4-y3;
2017 v2a[2] = z4-z3;
2018
2019 v2b[0] = x3-x2;
2020 v2b[1] = y3-y2;
2021 v2b[2] = z3-z2;
2022
2023 cross(v2a, v2b, v2);
2024 cross(v1, v2, v3);
2025
2026 norm(v1);
2027 norm(v2);
2028 norm(v3);
2029
2030 // gather 10 closest rotamer conformations...
2031
2032 for (j=0;j<10;j++)
2033 SORTED_ROTAMERS[i][j][0] = 500.;
2034
2035 j = 0;
2036 besthit = 1000.;
2037 bestpos = 0;
2038 do {
2039 if (rot_stat_idx[j][0]==res->type) {
2040 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);
2041 if (hit<SORTED_ROTAMERS[i][9][0]) {
2042 k = 9;
2043 while (k>=0 && hit<SORTED_ROTAMERS[i][k][0]) {
2044 k--;
2045 }
2046 k++;
2047 // k = hit
2048 for (l=9;l>k;l--) {
2049 SORTED_ROTAMERS[i][l][0]=SORTED_ROTAMERS[i][l-1][0];
2050 SORTED_ROTAMERS[i][l][1]=SORTED_ROTAMERS[i][l-1][1];
2051 }
2052 SORTED_ROTAMERS[i][k][0]=hit;
2053 SORTED_ROTAMERS[i][k][1]=j;
2054 }
2055 }
2056 j++;
2057 } while (rot_stat_idx[j][0]>=0);
2058
2059
2060 besthit = SORTED_ROTAMERS[i][0][0];
2061 bestpos = SORTED_ROTAMERS[i][0][1];
2062
2063
2064 // new rebuild...
2065
2066 pos = rot_stat_idx[bestpos][5];
2067 nsc = nheavy[res->type]+1;
2068
2069 if (_PDB_SG) { // more than one rotamer - check SC
2070 bestdd = 100.; crot = 0;
2071 for (l=0;l<2;l++) { // check two closest conformations
2072 cpos = SORTED_ROTAMERS[i][l][1];
2073 for (m=0;m<rot_stat_idx[cpos][4];m++) {
2074 for (j=0;j<3;j++) {
2075 vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j];
2076 for (k=0;k<3;k++) {
2077 if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.;
2078 }
2079 }
2080 pos = rot_stat_idx[cpos][5]+nsc*m;
2081 for (j=0;j<nsc;j++) {
2082 for (k=0;k<3;k++) {
2083 sc[j][k] = rot_stat_coords[pos+j][k];
2084 }
2085 }
2086 superimpose2(vv,lsys,3,sc,nsc);
2087 for (j=0;j<nsc;j++) {
2088 sc[j][0] += x3;
2089 sc[j][1] += y3;
2090 sc[j][2] += z3;
2091 }
2092 cmx = 0.; cmy = 0.; cmz = 0.;
2093 for (j=0;j<nsc;j++) {
2094 cmx += sc[j][0];
2095 cmy += sc[j][1];
2096 cmz += sc[j][2];
2097 }
2098 cmx /= (real) nsc;
2099 cmy /= (real) nsc;
2100 cmz /= (real) nsc;
2101 ddx = res->cmx-cmx;
2102 ddy = res->cmy-cmy;
2103 ddz = res->cmz-cmz;
2104 ddx *= ddx;
2105 ddy *= ddy;
2106 ddz *= ddz;
2107 ddd = ddx+ddy+ddz;
2108 if (ddd<bestdd) {
2109 bestdd = ddd;
2110 crot = pos; // closest rotamer position
2111 }
2112 }
2113 }
2114 pos = crot;
2115 } // PDB_SG
2116
2117 for (j=0;j<3;j++) {
2118 vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j];
2119 for (k=0;k<3;k++) {
2120 if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.;
2121 }
2122 }
2123
2124 for (j=0;j<nsc;j++) {
2125 for (k=0;k<3;k++) {
2126 sc[j][k] = rot_stat_coords[pos+j][k];
2127 }
2128 }
2129
2130 superimpose2(vv,lsys,3,sc,nsc);
2131
2132 for (j=0;j<nsc;j++) {
2133 sc[j][0] += x3;
2134 sc[j][1] += y3;
2135 sc[j][2] += z3;
2136 }
2137
2138 for (j=1;j<nsc;j++) {
2139 add_replace(res, heavy_atoms[10*res->type+j-1], sc[j][0], sc[j][1], sc[j][2], FLAG_SIDECHAIN);
2140 }
2141
2142 if (res->next) res = res->next;
2143
2144 } // i++, next res
2145
2146 for (i=0;i<12;i++)
2147 free(sc[i]);
2148 for (i=0;i<3;i++) {
2149 free(lsys[i]);
2150 free(vv[i]);
2151 }
2152 free(sc);
2153 free(lsys);
2154 free(vv);
2155 }
2156
2157
2158 typedef struct _atom_list {
2159 atom_type *atom;
2160 struct _atom_list *next;
2161 } atom_list;
2162
2163 int get_conflicts(res_type *res, atom_list ****grid, int xgrid, int ygrid, int zgrid)
2164 {
2165 atom_list *llist;
2166 atom_type *atom, *atom2;
2167 int i, j, k, x, y, z;
2168 int ii, jj, kk, con, iter, maxcon, merged;
2169 real dx, dy, dz, dd;
2170
2171 con = 0;
2172 atom = res->atoms;
2173 while (atom) {
2174 i = atom->gx;
2175 j = atom->gy;
2176 k = atom->gz;
2177 for (ii=i-2;ii<=i+2;ii++)
2178 for (jj=j-2;jj<=j+2;jj++)
2179 for (kk=k-2;kk<=k+2;kk++) {
2180 if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<zgrid) {
2181 llist = grid[ii][jj][kk];
2182 while (llist) {
2183 atom2 = llist->atom;
2184 if (atom && atom2 && res && atom2->res) {
2185 merged=0;
2186 if (res==atom2->res) { // self-xvol
2187 if (atom->flag & FLAG_SIDECHAIN && atom2->flag & FLAG_SIDECHAIN) merged=1;
2188 if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1;
2189 if (atom->name[0]=='C' && atom->name[1]=='A' && atom2->name[0]=='C' && atom2->name[1]=='B') merged=1;
2190 if (atom->name[0]=='C' && atom->name[1]=='B' && atom2->name[0]=='C' && atom2->name[1]=='A') merged=1;
2191 if (res->name[0]=='P') {
2192 if (atom->name[0]=='C' && atom->name[1]=='D' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1;
2193 if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]=='D') merged=1;
2194 }
2195
2196 if (!merged) {
2197 // printf("merged: %s[%d] %s-%s %d %d\n", res->name,res->num,atom->name,atom2->name,atom->flag,atom2->flag);
2198 }
2199 } else
2200 if (res->next==atom2->res || res==atom2->res->next) {
2201 if (atom->name[0]=='C' && atom->name[1]==' ' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1;
2202 if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]==' ') merged=1;
2203 }
2204 if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1; // for now
2205 if (atom->flag & FLAG_SCM || atom2->flag & FLAG_SCM) merged=1; // for now
2206 if (!merged) {
2207 dx = atom->x-atom2->x;
2208 dx*=dx;
2209 dy = atom->y-atom2->y;
2210 dy*=dy;
2211 dz = atom->z-atom2->z;
2212 dz*=dz;
2213 dd = dx+dy+dz;
2214 if (dd<_SG_XVOL_DIST*_SG_XVOL_DIST) {
2215 con++;
2216 }
2217 }
2218 }
2219 llist = llist->next;
2220 }
2221 }
2222 }
2223 atom = atom->next;
2224 }
2225
2226 return con;
2227 }
2228
2229 int display_conflicts(res_type *res, atom_list ****grid, int xgrid, int ygrid, int zgrid)
2230 {
2231 atom_list *llist;
2232 atom_type *atom, *atom2;
2233 int i, j, k, x, y, z;
2234 int ii, jj, kk, con, iter, maxcon, merged;
2235 real dx, dy, dz, dd;
2236
2237 con = 0;
2238 atom = res->atoms;
2239 while (atom) {
2240 i = atom->gx;
2241 j = atom->gy;
2242 k = atom->gz;
2243 for (ii=i-2;ii<=i+2;ii++)
2244 for (jj=j-2;jj<=j+2;jj++)
2245 for (kk=k-2;kk<=k+2;kk++) {
2246 if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<zgrid) {
2247 llist = grid[ii][jj][kk];
2248 while (llist) {
2249 atom2 = llist->atom;
2250 if (atom && atom2 && res && atom2->res) {
2251 merged=0;
2252 if (res==atom2->res) { // self-xvol
2253 if (atom->flag & FLAG_SIDECHAIN && atom2->flag & FLAG_SIDECHAIN) merged=1;
2254 if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1;
2255 if (atom->name[0]=='C' && atom->name[1]=='A' && atom2->name[0]=='C' && atom2->name[1]=='B') merged=1;
2256 if (atom->name[0]=='C' && atom->name[1]=='B' && atom2->name[0]=='C' && atom2->name[1]=='A') merged=1;
2257 if (res->name[0]=='P') {
2258 if (atom->name[0]=='C' && atom->name[1]=='D' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1;
2259 if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]=='D') merged=1;
2260 }
2261 if (!merged) {
2262 // printf("merged: %s[%d] %s-%s %d %d\n", res->name,res->num,atom->name,atom2->name,atom->flag,atom2->flag);
2263 }
2264 } else
2265 if (res->next==atom2->res || res==atom2->res->next) {
2266 if (atom->name[0]=='C' && atom->name[1]==' ' && atom2->name[0]=='N' && atom2->name[1]==' ') merged=1;
2267 if (atom->name[0]=='N' && atom->name[1]==' ' && atom2->name[0]=='C' && atom2->name[1]==' ') merged=1;
2268 }
2269 if (atom->flag & FLAG_BACKBONE && atom2->flag & FLAG_BACKBONE) merged=1; // for now
2270 if (atom->flag & FLAG_SCM || atom2->flag & FLAG_SCM) merged=1; // for now
2271 if (!merged) {
2272 dx = atom->x-atom2->x;
2273 dx*=dx;
2274 dy = atom->y-atom2->y;
2275 dy*=dy;
2276 dz = atom->z-atom2->z;
2277 dz*=dz;
2278 dd = dx+dy+dz;
2279 if (dd<1.6*1.6) {
2280 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);
2281 con++;
2282 }
2283 }
2284 }
2285 llist = llist->next;
2286 }
2287 }
2288 }
2289 atom = atom->next;
2290 }
2291
2292 return con;
2293 }
2294
2295
2296 void allocate_grid(atom_list *****grid_, int *xgrid_, int *ygrid_, int *zgrid_)
2297 {
2298 static int xgrid, ygrid, zgrid;
2299 static atom_list ****grid = NULL;
2300 atom_list *llist, *alist;
2301 real min[3], max[3];
2302 res_type *res, *worst;
2303 atom_type *atom, *atom2;
2304 int i, j, x, y, z;
2305
2306 if (!grid && chain->residua && chain->residua->atoms) {
2307 res = chain->residua;
2308 min[0]=max[0]=res->atoms->x;
2309 min[1]=max[1]=res->atoms->y;
2310 min[2]=max[2]=res->atoms->z;
2311 while (res) {
2312 atom = res->atoms;
2313 while (atom) {
2314 if (atom->x<min[0]) min[0]=atom->x;
2315 if (atom->y<min[1]) min[1]=atom->y;
2316 if (atom->z<min[2]) min[2]=atom->z;
2317 if (atom->x>max[0]) max[0]=atom->x;
2318 if (atom->y>max[1]) max[1]=atom->y;
2319 if (atom->z>max[2]) max[2]=atom->z;
2320 atom = atom->next;
2321 }
2322 res = res->next;
2323 }
2324
2325 xgrid = (max[0]-min[0])/GRID_RES;
2326 ygrid = (max[1]-min[1])/GRID_RES;
2327 zgrid = (max[2]-min[2])/GRID_RES;
2328
2329 if (_VERBOSE) printf("Allocating grid (%d %d %d)...\n", xgrid, ygrid, zgrid);
2330
2331 grid = (atom_list****)calloc(sizeof(atom_list***)*(xgrid+1),1);
2332 for (i=0;i<xgrid+1;i++) {
2333 grid[i] = (atom_list***)calloc(sizeof(atom_list**)*(ygrid+1),1);
2334 for (j=0;j<ygrid+1;j++) {
2335 grid[i][j] = (atom_list**)calloc(sizeof(atom_list*)*(zgrid+1),1);
2336 }
2337 }
2338
2339 res = chain->residua;
2340 while (res) {
2341 atom = res->atoms;
2342 while (atom) {
2343 x = xgrid*(atom->x-min[0])/(max[0]-min[0]);
2344 y = ygrid*(atom->y-min[1])/(max[1]-min[1]);
2345 z = zgrid*(atom->z-min[2])/(max[2]-min[2]);
2346 alist = (atom_list*)calloc(sizeof(atom_list),1);
2347 alist->atom = atom;
2348 atom->gx = x;
2349 atom->gy = y;
2350 atom->gz = z;
2351 if (grid[x][y][z]!=NULL) {
2352 llist = grid[x][y][z];
2353 while (llist->next) llist=llist->next;
2354 llist->next = alist;
2355 } else {
2356 grid[x][y][z]=alist;
2357 }
2358 atom = atom->next;
2359 }
2360 res = res->next;
2361 }
2362 } else {
2363 if (_VERBOSE) printf("Grid already allocated (%d %d %d)\n", xgrid, ygrid, zgrid);
2364 }
2365
2366 *grid_ = grid;
2367 *xgrid_ = xgrid;
2368 *ygrid_ = ygrid;
2369 *zgrid_ = zgrid;
2370 }
2371
2372 void optimize_exvol(void)
2373 {
2374 real min[3], max[3];
2375 res_type *res, *worst;
2376 atom_type *atom, *atom2;
2377 int xgrid, ygrid, zgrid;
2378 atom_list ****grid, *llist, *alist;
2379 int i, j, k, l, m, x, y, z;
2380 int ii, jj, kk, con, iter, maxcon, totcon;
2381 int cpos, bestpos, pos, con0;
2382 real v1[3], v2a[3], v2b[3], v2[3], v3[3];
2383 int nsc, nca;
2384 real cax, cay, caz;
2385 real **lsys, **vv, **sc;
2386 real x1, y1, z1;
2387 real x2, y2, z2;
2388 real x3, y3, z3;
2389 real x4, y4, z4;
2390
2391 min[0]=1e5;
2392 min[1]=1e5;
2393 min[2]=1e5;
2394 max[0]=-1e5;
2395 max[1]=-1e5;
2396 max[2]=-1e5;
2397
2398 lsys = (real**)calloc(sizeof(real*)*3,1);
2399 vv = (real**)calloc(sizeof(real*)*3,1);
2400 sc = (real**)calloc(sizeof(real*)*12,1);
2401 for (i=0;i<12;i++)
2402 sc[i] = (real*)calloc(sizeof(real)*3,1);
2403 for (i=0;i<3;i++) {
2404 lsys[i] = (real*)calloc(sizeof(real)*3,1);
2405 vv[i] = (real*)calloc(sizeof(real)*3,1);
2406 }
2407
2408 allocate_grid(&grid, &xgrid, &ygrid, &zgrid);
2409
2410 if (_VERBOSE) printf("Finding excluded volume conflicts...\n");
2411
2412 iter = 0;
2413
2414 do {
2415 //printf("ITER: %d\n", iter);
2416
2417 maxcon = 0;
2418 totcon=0;
2419
2420 res = chain->residua;
2421 while (res) {
2422 if (res->protein) {
2423 con = get_conflicts(res, grid, xgrid, ygrid, zgrid);
2424 if (con>0) {
2425 totcon+=con;
2426 if (con>maxcon) {
2427 maxcon = con;
2428 worst = res;
2429 }
2430 }
2431 }
2432 res = res->next;
2433 }
2434
2435 if (_VERBOSE && iter==0) {
2436 printf("Total number of conflicts: %d\n", totcon);
2437 }
2438
2439 if (totcon==0) break;
2440
2441 if (_VERBOSE && iter==0) {
2442 printf("Maximum number of conflicts: %s[%d] : %d\n", worst->name, worst->num, maxcon);
2443 }
2444
2445 totcon=0;
2446
2447 if (maxcon>0) {
2448
2449 // try to fix...
2450
2451 res = chain->residua;
2452 for (i=0;i<chain_length;i++) {
2453 if (!strncmp(res->name,"GLY",3) || !res->protein) {
2454 if (res->next) res = res->next;
2455 continue;
2456 }
2457
2458 nsc = nheavy[res->type]+1;
2459
2460 x1 = C_ALPHA[i-2][0];
2461 y1 = C_ALPHA[i-2][1];
2462 z1 = C_ALPHA[i-2][2];
2463 x2 = C_ALPHA[i-1][0];
2464 y2 = C_ALPHA[i-1][1];
2465 z2 = C_ALPHA[i-1][2];
2466 x3 = C_ALPHA[i][0];
2467 y3 = C_ALPHA[i][1];
2468 z3 = C_ALPHA[i][2];
2469 x4 = C_ALPHA[i+1][0];
2470 y4 = C_ALPHA[i+1][1];
2471 z4 = C_ALPHA[i+1][2];
2472
2473 v1[0] = x4-x2;
2474 v1[1] = y4-y2;
2475 v1[2] = z4-z2;
2476
2477 v2a[0] = x4-x3;
2478 v2a[1] = y4-y3;
2479 v2a[2] = z4-z3;
2480
2481 v2b[0] = x3-x2;
2482 v2b[1] = y3-y2;
2483 v2b[2] = z3-z2;
2484
2485 cross(v2a, v2b, v2);
2486 cross(v1, v2, v3);
2487
2488 norm(v1);
2489 norm(v2);
2490 norm(v3);
2491
2492 con = get_conflicts(res, grid, xgrid, ygrid, zgrid);
2493
2494 if (con>0) {
2495
2496 bestpos=0;
2497 con0 = 100;
2498 for (l=0;l<10;l++) { // check two closest conformations
2499 cpos = SORTED_ROTAMERS[i][l][1];
2500 for (m=0;m<rot_stat_idx[cpos][4];m++) {
2501 for (j=0;j<3;j++) {
2502 vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j];
2503 for (k=0;k<3;k++) {
2504 if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.;
2505 }
2506 }
2507 pos = rot_stat_idx[cpos][5]+nsc*m;
2508 for (j=0;j<nsc;j++) {
2509 for (k=0;k<3;k++) {
2510 sc[j][k] = rot_stat_coords[pos+j][k];
2511 }
2512 }
2513 superimpose2(vv,lsys,3,sc,nsc);
2514 for (j=0;j<nsc;j++) {
2515 sc[j][0] += x3;
2516 sc[j][1] += y3;
2517 sc[j][2] += z3;
2518 }
2519 for (j=1;j<nsc;j++) {
2520 add_replace(res, heavy_atoms[10*res->type+j-1], sc[j][0], sc[j][1], sc[j][2], FLAG_SIDECHAIN);
2521 }
2522 con = get_conflicts(res, grid, xgrid, ygrid, zgrid);
2523 //printf("test: %d\n", con);
2524
2525 if (con<con0) {
2526 con0 = con;
2527 bestpos = pos;
2528 }
2529 if (con==0) break;
2530 }
2531 if (con==0) break;
2532 }
2533
2534 totcon += con0;
2535
2536 for (j=0;j<3;j++) {
2537 vv[0][j] = v1[j]; vv[1][j] = v2[j]; vv[2][j] = v3[j];
2538 for (k=0;k<3;k++) {
2539 if (j==k) lsys[j][k]=1.; else lsys[j][k]=0.;
2540 }
2541 }
2542 pos = bestpos;
2543 for (j=0;j<nsc;j++) {
2544 for (k=0;k<3;k++) {
2545 sc[j][k] = rot_stat_coords[pos+j][k];
2546 }
2547 }
2548 superimpose2(vv,lsys,3,sc,nsc);
2549 for (j=0;j<nsc;j++) {
2550 sc[j][0] += x3;
2551 sc[j][1] += y3;
2552 sc[j][2] += z3;
2553 }
2554 for (j=1;j<nsc;j++) {
2555 add_replace(res, heavy_atoms[10*res->type+j-1], sc[j][0], sc[j][1], sc[j][2], FLAG_SIDECHAIN);
2556 }
2557 }
2558
2559
2560
2561 res=res->next;
2562
2563 } // i
2564 }
2565
2566 iter++;
2567
2568 } while (iter<_XVOL_ITER);
2569
2570
2571 if (_VERBOSE) {
2572 if (totcon>0)
2573 printf("WARNING: %d steric conflict(s) are still there.\n", totcon);
2574 else
2575 printf("All steric conflicts removed.\n");
2576 }
2577
2578 for (i=0;i<12;i++)
2579 free(sc[i]);
2580 for (i=0;i<3;i++) {
2581 free(lsys[i]);
2582 free(vv[i]);
2583 }
2584 free(sc);
2585 free(lsys);
2586 free(vv);
2587
2588
2589 }
2590
2591
2592 void vcross(real ax,real ay,real az,real bx,real by,real bz,real *cx,real *cy,real *cz)
2593 {
2594 *cx = ay * bz - by * az;
2595 *cy = az * bx - bz * ax;
2596 *cz = ax * by - bx * ay;
2597 }
2598
2599 real vdot(real ax,real ay,real az,real bx,real by,real bz)
2600 {
2601 return ax*bx+ay*by+az*bz;
2602 }
2603
2604 real calc_torsion(atom_type *a1, atom_type *a2, atom_type *a3, atom_type *a4)
2605 {
2606 real v12x, v12y, v12z;
2607 real v43x, v43y, v43z;
2608 real zx, zy, zz;
2609 real px, py, pz;
2610 real xx, xy, xz;
2611 real yx, yy, yz;
2612 real u, v, angle;
2613
2614 v12x = a1->x-a2->x;
2615 v12y = a1->y-a2->y;
2616 v12z = a1->z-a2->z;
2617
2618 v43x = a4->x-a3->x;
2619 v43y = a4->y-a3->y;
2620 v43z = a4->z-a3->z;
2621
2622 zx = a2->x-a3->x;
2623 zy = a2->y-a3->y;
2624 zz = a2->z-a3->z;
2625
2626 vcross(zx,zy,zz,v12x,v12y,v12z,&px,&py,&pz);
2627 vcross(zx,zy,zz,v43x,v43y,v43z,&xx,&xy,&xz);
2628 vcross(zx,zy,zz,xx,xy,xz,&yx,&yy,&yz);
2629
2630 u = vdot(xx,xy,xz,xx,xy,xz);
2631 v = vdot(yx,yy,yz,yx,yy,yz);
2632
2633 angle = 360.;
2634
2635 if (u<0. || v<0.) return angle;
2636
2637 u = vdot(px,py,pz,xx,xy,xz) / sqrt(u);
2638 v = vdot(px,py,pz,yx,yy,yz) / sqrt(v);
2639
2640 if (u != 0.0 || v != 0.0) angle = atan2(v, u) * RADDEG;
2641
2642
2643 return angle;
2644
2645 }
2646
2647
2648 // Ca-N-C-Cb angle should be close to 34 deg
2649
2650 void chirality_check(void)
2651 {
2652 int i;
2653 atom_type *a_ca, *a_n, *a_c, *a_cb;
2654 atom_type *atom;
2655 res_type *res;
2656 real angle;
2657 real nx, ny, nz;
2658 real px, py, pz;
2659 real qx, qy, qz;
2660 real rx, ry, rz;
2661 real xx, xy, xz;
2662 real yx, yy, yz;
2663 real dd, costheta, sintheta;
2664
2665 if (_VERBOSE) printf("Checking chirality...\n");
2666 res = chain->residua;
2667 while (res) {
2668 a_ca = a_n = a_c = a_cb = NULL;
2669 a_ca = find_atom(res,"CA ");
2670 a_n = find_atom(res,"N ");
2671 a_c = find_atom(res,"C ");
2672 a_cb = find_atom(res,"CB ");
2673 if (a_ca && a_n && a_c && a_cb) {
2674 angle = calc_torsion(a_ca, a_n, a_c, a_cb);
2675 if (angle<0.) {
2676 if (_VERBOSE) printf("WARNING: D-aa detected at %s %3d : %5.2f", res->name, res->num, angle);
2677 xx = a_ca->x-a_n->x;
2678 xy = a_ca->y-a_n->y;
2679 xz = a_ca->z-a_n->z;
2680 yx = a_c->x-a_ca->x;
2681 yy = a_c->y-a_ca->y;
2682 yz = a_c->z-a_ca->z;
2683 vcross(xx,xy,xz,yx,yy,yz,&nx,&ny,&nz);
2684 dd = sqrt(nx*nx+ny*ny+nz*nz);
2685 nx /= dd;
2686 ny /= dd;
2687 nz /= dd;
2688 // nx, ny, nz = reflection plane normal
2689 rx = xx-yx;
2690 ry = xy-yy;
2691 rz = xz-yz;
2692 dd = sqrt(rx*rx+ry*ry+rz*rz);
2693 rx /= dd;
2694 ry /= dd;
2695 rz /= dd;
2696 costheta = -1.;
2697 sintheta = 0.;
2698 atom = res->atoms;
2699 while (atom) {
2700 if (atom->flag & FLAG_SIDECHAIN) {
2701 px = atom->x-a_ca->x;
2702 py = atom->y-a_ca->y;
2703 pz = atom->z-a_ca->z;
2704 qx = qy = qz = 0.;
2705 qx += (costheta + (1 - costheta) * rx * rx) * px;
2706 qx += ((1 - costheta) * rx * ry - rz * sintheta) * py;
2707 qx += ((1 - costheta) * rx * rz + ry * sintheta) * pz;
2708 qy += ((1 - costheta) * rx * ry + rz * sintheta) * px;
2709 qy += (costheta + (1 - costheta) * ry * ry) * py;
2710 qy += ((1 - costheta) * ry * rz - rx * sintheta) * pz;
2711 qz += ((1 - costheta) * rx * rz - ry * sintheta) * px;
2712 qz += ((1 - costheta) * ry * rz + rx * sintheta) * py;
2713 qz += (costheta + (1 - costheta) * rz * rz) * pz;
2714 qx += a_ca->x;
2715 qy += a_ca->y;
2716 qz += a_ca->z;
2717 atom->x = qx;
2718 atom->y = qy;
2719 atom->z = qz;
2720 }
2721 atom = atom->next;
2722 }
2723 angle = calc_torsion(a_ca, a_n, a_c, a_cb);
2724 if (_VERBOSE) printf(", fixed : %5.2f\n", angle);
2725 }
2726 }
2727 res = res->next;
2728 }
2729 }
2730
2731
2732 #endif
2733
2734 real hb_energy(res_type *res, atom_list ****grid, int xgrid, int ygrid, int zgrid)
2735 {
2736 atom_type *atom, *c_atom1, *o_atom1, *n_atom1, *c_atom2, *o_atom2, *n_atom2, *tmp_atom;
2737 atom_type h_atom;
2738 int i, j, k, ii, jj, kk;
2739 atom_list *llist, *alist;
2740 real dx, dy, dz, dist, min_dist1, min_dist2;
2741 real hx1, hy1, hz1, dd;
2742 real dno, dnc, dho, dhc;
2743 real ene, Q;
2744
2745 ene = 1e3;
2746
2747 if (!res || !res->prev) return ene;
2748
2749 Q = -27888.0; // DSSP h-bond energy constant
2750
2751 c_atom1 = o_atom1 = n_atom1 = NULL;
2752
2753 atom = res->prev->atoms;
2754 while (atom) {
2755 if (atom->name[0]=='C' && atom->name[1]==' ') c_atom1 = atom;
2756 if (atom->name[0]=='O' && atom->name[1]==' ') o_atom1 = atom;
2757 atom = atom->next;
2758 }
2759
2760 atom = res->atoms;
2761 while (atom) {
2762 if (atom->name[0]=='N' && atom->name[1]==' ') { n_atom1 = atom; break; }
2763 atom = atom->next;
2764 }
2765
2766 // first bond
2767
2768 min_dist2 = 1e10;
2769 o_atom2 = c_atom2 = NULL;
2770 if (n_atom1) {
2771 i = n_atom1->gx;
2772 j = n_atom1->gy;
2773 k = n_atom1->gz;
2774 for (ii=i-1;ii<=i+1;ii++) {
2775 for (jj=j-1;jj<=j+1;jj++) {
2776 for (kk=k-1;kk<=k+1;kk++) {
2777 if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<=zgrid) {
2778 llist = grid[ii][jj][kk];
2779 while (llist) {
2780 if (llist->atom->name[0]=='O' && llist->atom->name[1]==' ' && abs(llist->atom->res->locnum-n_atom1->res->locnum)>2) {
2781 tmp_atom = llist->atom;
2782 dx = n_atom1->x-tmp_atom->x;
2783 dy = n_atom1->y-tmp_atom->y;
2784 dz = n_atom1->z-tmp_atom->z;
2785 dist = dx*dx+dy*dy+dz*dz;
2786 if (dist<min_dist2 && dist<25.0) {
2787 o_atom2=tmp_atom;
2788 min_dist2 = dist;
2789 }
2790 }
2791 llist = llist->next;
2792 }
2793 }
2794 }
2795 }
2796 }
2797 }
2798
2799 if (o_atom2) {
2800 atom = o_atom2->res->atoms;
2801 while (atom) {
2802 if (atom->name[0]=='C' && atom->name[1]==' ') { c_atom2 = atom; break; }
2803 atom = atom->next;
2804 }
2805 if (c_atom2) {
2806 hx1 = o_atom1->x-c_atom1->x;
2807 hy1 = o_atom1->y-c_atom1->y;
2808 hz1 = o_atom1->z-c_atom1->z;
2809 dd = -1.081f/sqrt(hx1*hx1+hy1*hy1+hz1*hz1);
2810 hx1 *= dd;
2811 hy1 *= dd;
2812 hz1 *= dd;
2813
2814 hx1 += n_atom1->x;
2815 hy1 += n_atom1->y;
2816 hz1 += n_atom1->z;
2817
2818 add_replace(n_atom1->res, "H ", hx1, hy1, hz1, FLAG_BACKBONE);
2819
2820 // dno
2821 dx = n_atom1->x-o_atom2->x;
2822 dy = n_atom1->y-o_atom2->y;
2823 dz = n_atom1->z-o_atom2->z;
2824 dno = sqrt(dx*dx+dy*dy+dz*dz);
2825
2826 // dnc
2827 dx = n_atom1->x-c_atom2->x;
2828 dy = n_atom1->y-c_atom2->y;
2829 dz = n_atom1->z-c_atom2->z;
2830 dnc = sqrt(dx*dx+dy*dy+dz*dz);
2831
2832 // dho
2833 dx = hx1-o_atom2->x;
2834 dy = hy1-o_atom2->y;
2835 dz = hz1-o_atom2->z;
2836 dho = sqrt(dx*dx+dy*dy+dz*dz);
2837
2838 // dhc
2839 dx = hx1-c_atom2->x;
2840 dy = hy1-c_atom2->y;
2841 dz = hz1-c_atom2->z;
2842 dhc = sqrt(dx*dx+dy*dy+dz*dz);
2843
2844 if (dho<0.01F || dhc<0.01F || dnc<0.01F || dno<0.01F) {
2845 ene = -10.0;
2846 } else {
2847 ene = 0.001*(Q/dho - Q/dhc + Q/dnc - Q/dno);
2848 }
2849 }
2850 }
2851
2852 /******
2853 // second bond
2854
2855 min_dist2 = 1e10;
2856 n_atom2 = NULL;
2857 if (n_atom1) {
2858 i = o_atom1->gx;
2859 j = o_atom1->gy;
2860 k = o_atom1->gz;
2861 for (ii=i-1;ii<=i+1;ii++) {
2862 for (jj=j-1;jj<=j+1;jj++) {
2863 for (kk=k-1;kk<=k+1;kk++) {
2864 if (ii>=0 && ii<xgrid && jj>=0 && jj<ygrid && kk>=0 && kk<=zgrid) {
2865 llist = grid[ii][jj][kk];
2866 while (llist) {
2867 if (llist->atom->name[0]=='N' && llist->atom->name[1]==' ' && (abs(llist->atom->res->locnum-n_atom1->res->locnum)>2)) {
2868 tmp_atom = llist->atom;
2869 if (tmp_atom->res!=c_atom2->res) {
2870 dx = o_atom1->x-tmp_atom->x;
2871 dy = o_atom1->y-tmp_atom->y;
2872 dz = o_atom1->z-tmp_atom->z;
2873 dist = dx*dx+dy*dy+dz*dz;
2874 if (dist<min_dist2 && dist<25.0) {
2875 n_atom2=tmp_atom;
2876 min_dist2 = dist;
2877 }
2878 }
2879 }
2880 llist = llist->next;
2881 }
2882 }
2883 }
2884 }
2885 }
2886 }
2887
2888 if (n_atom2) {
2889 c_atom2 = o_atom2 = NULL;
2890 atom = n_atom2->res->atoms;
2891 while (atom) {
2892 if (atom->name[0]=='C' && atom->name[1]==' ') { c_atom2 = atom; }
2893 if (atom->name[0]=='O' && atom->name[1]==' ') { c_atom2 = atom; }
2894 atom = atom->next;
2895 }
2896
2897 if (c_atom2) {
2898 hx1 = o_atom1->x-c_atom1->x;
2899 hy1 = o_atom1->y-c_atom1->y;
2900 hz1 = o_atom1->z-c_atom1->z;
2901 dd = -1.081f/sqrt(hx1*hx1+hy1*hy1+hz1*hz1);
2902 hx1 *= dd;
2903 hy1 *= dd;
2904 hz1 *= dd;
2905
2906 hx1 += n_atom1->x;
2907 hy1 += n_atom1->y;
2908 hz1 += n_atom1->z;
2909
2910 }
2911 *******/
2912
2913
2914 return ene;
2915 }
2916
2917 // rotates a point around a vector
2918 void rot_point_vector(real *x, real *y, real *z, real u, real v, real w, real angle)
2919 {
2920 real ux, uy, uz, vx, vy, vz, wx, wy, wz, sa, ca;
2921
2922 sa = sinf(10.0*M_PI*angle/180.0);
2923 ca = cosf(10.0*M_PI*angle/180.0);
2924
2925 ux = u**x;
2926 uy = u**y;
2927 uz = u**z;
2928 vx = v**x;
2929 vy = v**y;
2930 vz = v**z;
2931 wx = w**x;
2932 wy = w**y;
2933 wz = w**z;
2934
2935 *x = u*(ux+vy+wz)+(*x*(v*v+w*w)-u*(vy+wz))*ca+(-wy+vz)*sa;
2936 *y = v*(ux+vy+wz)+(*y*(u*u+w*w)-v*(ux+wz))*ca+( wx-uz)*sa;
2937 *z = w*(ux+vy+wz)+(*z*(u*u+v*v)-w*(ux+vy))*ca+(-vx+uy)*sa;
2938 }
2939
2940
2941 // rotates a peptide plate
2942
2943 void rot_peptide(res_type *res, real angle)
2944 {
2945 atom_type *atom, *c_atom, *o_atom, *n_atom, *ca_atom1, *ca_atom2;
2946 real u, v, w, x, y, z, dd;
2947
2948 if (!res || !res->prev) return;
2949
2950 c_atom = o_atom = n_atom = ca_atom1 = ca_atom2 = NULL;
2951
2952 atom = res->prev->atoms;
2953 while (atom) {
2954 if (atom->name[0]=='C' && atom->name[1]=='A') ca_atom1 = atom;
2955 if (atom->name[0]=='C' && atom->name[1]==' ') c_atom = atom;
2956 if (atom->name[0]=='O' && atom->name[1]==' ') o_atom = atom;
2957 atom = atom->next;
2958 }
2959
2960 atom = res->atoms;
2961 while (atom) {
2962 if (atom->name[0]=='C' && atom->name[1]=='A') ca_atom2 = atom;
2963 if (atom->name[0]=='N' && atom->name[1]==' ') n_atom = atom;
2964 atom = atom->next;
2965 }
2966
2967 if (c_atom && o_atom && n_atom && ca_atom1 && ca_atom2) {
2968 u = ca_atom2->x-ca_atom1->x;
2969 v = ca_atom2->y-ca_atom1->y;
2970 w = ca_atom2->z-ca_atom1->z;
2971 dd = 1.0f/sqrt(u*u+v*v+w*w);
2972 u*=dd; v*=dd; w*=dd; // normalize ca-ca vector
2973 x = n_atom->x-ca_atom1->x;
2974 y = n_atom->y-ca_atom1->y;
2975 z = n_atom->z-ca_atom1->z;
2976 rot_point_vector(&x, &y, &z, u, v, w, angle);
2977 n_atom->x = x+ca_atom1->x;
2978 n_atom->y = y+ca_atom1->y;
2979 n_atom->z = z+ca_atom1->z;
2980 x = c_atom->x-ca_atom1->x;
2981 y = c_atom->y-ca_atom1->y;
2982 z = c_atom->z-ca_atom1->z;
2983 rot_point_vector(&x, &y, &z, u, v, w, angle);
2984 c_atom->x = x+ca_atom1->x;
2985 c_atom->y = y+ca_atom1->y;
2986 c_atom->z = z+ca_atom1->z;
2987 x = o_atom->x-ca_atom1->x;
2988 y = o_atom->y-ca_atom1->y;
2989 z = o_atom->z-ca_atom1->z;
2990 rot_point_vector(&x, &y, &z, u, v, w, angle);
2991 o_atom->x = x+ca_atom1->x;
2992 o_atom->y = y+ca_atom1->y;
2993 o_atom->z = z+ca_atom1->z;
2994 }
2995
2996 }
2997
2998 void optimize_backbone(mol_type *chain)
2999 {
3000 int xgrid, ygrid, zgrid;
3001 atom_list ****grid;
3002 atom_type *atom;
3003 res_type *res;
3004 real ene, min_ene, tot1, tot2;
3005 int i, k, best;
3006 FILE *out;
3007
3008 if (_VERBOSE) printf("Optimizing backbone...\n");
3009
3010 allocate_grid(&grid, &xgrid, &ygrid, &zgrid);
3011
3012 tot1 = tot2 = 0.0;
3013
3014 res = chain->residua;
3015 while (res) {
3016 ene = hb_energy(res, grid, xgrid, ygrid, zgrid);
3017 if (ene<-0.5) tot1 += ene;
3018 res = res->next;
3019 }
3020
3021 res = chain->residua;
3022 while (res) {
3023 if (res->type!=7) {
3024 ene = hb_energy(res, grid, xgrid, ygrid, zgrid);
3025 if (ene<1.0) { // try to optimize
3026 min_ene = ene;
3027 rot_peptide(res, -1.1);
3028 best = 0;
3029 for (i=-10;i<10;i++) {
3030 rot_peptide(res, 0.1);
3031 ene = hb_energy(res, grid, xgrid, ygrid, zgrid);
3032 if (ene<min_ene) {
3033 best = i;
3034 min_ene = ene;
3035 }
3036 }
3037 rot_peptide(res,-0.9);
3038 ene = hb_energy(res, grid, xgrid, ygrid, zgrid);
3039 if (min_ene<ene) {
3040 rot_peptide(res,0.1*best);
3041 ene = hb_energy(res, grid, xgrid, ygrid, zgrid);
3042 }
3043 }
3044 }
3045 res = res->next;
3046 }
3047
3048 res = chain->residua;
3049 while (res) {
3050 ene = hb_energy(res, grid, xgrid, ygrid, zgrid);
3051 if (ene<-0.5) tot2 += ene;
3052 res = res->next;
3053 }
3054
3055 if (_VERBOSE) printf("Backbone HB energy: before %g, after: %g, difference: %g\n", tot1, tot2, tot2-tot1);
3056
3057 }
3058
3059 int main(int argc, char **argv)
3060 {
3061 int i, j, next;
3062 char buf[100];
3063 char *name=NULL, *ini_name=NULL;
3064 char *ptr, out_name[1000];
3065 real f;
3066 mol_type *mol;
3067 struct timeb time0, time1;
3068
3069 for (i=1; i<argc; i++) {
3070 if (argv[i][0]=='-') {
3071 next=0;
3072 for (j=1; j<(int)strlen(argv[i]); j++) {
3073 switch(argv[i][j]) {
3074 case 'v': _VERBOSE=1; break;
3075 case 'c': _CA_OPTIMIZE=0; break;
3076 case 'e': _BB_REARRANGE=1; break;
3077 case 'r': _CA_RANDOM=1; break;
3078 case 'z': _CHIRAL=0; break;
3079 case 't': _CA_TRAJECTORY=1; break;
3080 case 'n': _CENTER_CHAIN=1; break;
3081 case 'b': _REBUILD_BB=0; break;
3082 case 's': _REBUILD_SC=0; break;
3083 case 'i': ini_name = argv[++i]; next=1; break;
3084 case 'g': _PDB_SG=1; break;
3085 case 'x': _TIME_SEED=1; break;
3086 case 'o': _XVOLUME=0; break;
3087 case 'h': _REBUILD_H=0; break;
3088 case 'q': _BB_OPTIMIZE=1; break;
3089 case 'p': _CISPRO=1; break;
3090 case 'u':
3091 if (sscanf(argv[++i],"%lf",&f)==1) {
3092 _CA_START_DIST = f;
3093 }
3094 next=1;
3095 break;
3096 default: {
3097 printf("Unknown option: %c\n", argv[i][j]);
3098 return -1;
3099 }
3100 }
3101 if (next) break;
3102 }
3103 } else {
3104 if (!name) name=argv[i];
3105 }
3106 }
3107
3108 if (!name) {
3109 printf("PULCHRA Protein Chain Restoration Algorithm version %4.2f\n", PULCHRA_VERSION);
3110 printf("Usage: %s [options] <pdb_file>\n", argv[0]);
3111 printf("The program default input is a PDB file.\n");
3112 printf("Output file <pdb_file.rebuild.pdb> will be created as a result.\n");
3113 printf("Valid options are:\n\n");
3114 printf(" -v : verbose output (default: off)\n");
3115 printf(" -n : center chain (default: off)\n");
3116 printf(" -x : time-seed random number generator (default: off)\n");
3117 printf(" -g : use PDBSG as an input format (CA=C-alpha, SC or CM=side chain c.m.)\n\n");
3118 printf(" -c : skip C-alpha positions optimization (default: on)\n");
3119 printf(" -p : detect cis-prolins (default: off)\n");
3120 printf(" -r : start from a random chain (default: off)\n");
3121 printf(" -i pdbfile : read the initial C-alpha coordinates from a PDB file\n");
3122 printf(" -t : save chain optimization trajectory to file <pdb_file.pdb.trajectory>\n");
3123 printf(" -u value : maximum shift from the restraint coordinates (default: 0.5A)\n\n");
3124 printf(" -e : rearrange backbone atoms (C, O are output after side chain) (default: off)\n");
3125
3126 #ifdef COMPILE_BB
3127 printf(" -b : skip backbone reconstruction (default: on)\n");
3128 printf(" -q : optimize backbone hydrogen bonds pattern (default: off)\n");
3129 printf(" -h : outputs hydrogen atoms (default: off)\n");
3130 #endif
3131
3132 #ifdef COMPILE_ROT
3133 printf(" -s : skip side chains reconstruction (default: on)\n");
3134 printf(" -o : don't attempt to fix excluded volume conflicts (default: on)\n");
3135 printf(" -z : don't check amino acid chirality (default: on)\n");
3136 #endif
3137 printf("\n");
3138 return -1;
3139 }
3140
3141 for (i=0; i<255; i++) /* prepare hash table*/
3142 AA_NUMS[i] = 20; /* dummy aa code */
3143 for (i=0; i<20; i++)
3144 AA_NUMS[(int)SHORT_AA_NAMES[i]] = i;
3145
3146 setbuf(stdout,0);
3147
3148 if (_TIME_SEED) srand(time(NULL)); else srand(1234);
3149
3150 if (_VERBOSE) printf("PULCHRA Protein Chain Restoration Algorithm version %4.2f\n", PULCHRA_VERSION);
3151
3152 ftime(&time0);
3153
3154 chain = new_mol();
3155
3156 if (read_pdb_file(name,chain,"chain")==FILE_NOT_FOUND) {
3157 if (_VERBOSE) printf("Can't read the input file!\n");
3158 return -1;
3159 }
3160
3161 if (_VERBOSE) printf("%d residua read.\n", chain->nres);
3162
3163 chain_length = chain->nres;
3164
3165 if (_CA_OPTIMIZE) {
3166 snprintf(out_name,1000,"%s.tra",name);
3167 ca_optimize(out_name, ini_name);
3168 }
3169
3170 #ifdef COMPILE_BB
3171 if (_REBUILD_BB) {
3172 rebuild_backbone();
3173 if (_BB_OPTIMIZE) {
3174 optimize_backbone(chain);
3175 }
3176 }
3177 #endif
3178
3179 #ifdef COMPILE_ROT
3180 if (_REBUILD_BB && _REBUILD_SC) {
3181 rebuild_sidechains();
3182 if (_XVOLUME)
3183 optimize_exvol();
3184 if (_CHIRAL)
3185 chirality_check();
3186 }
3187 #endif
3188
3189 if (_CENTER_CHAIN) {
3190 center_chain(chain);
3191 }
3192
3193
3194 if (_BB_REARRANGE) {
3195 if (_VERBOSE) printf("Rearranging backbone atoms...\n");
3196 }
3197
3198 ptr = strstr(name,".pdb");
3199 if (ptr) ptr[0]=0;
3200
3201 snprintf(out_name,1000,"%s.rebuilt.pdb",name);
3202
3203 if (_VERBOSE) printf("Writing output file %s...\n", out_name);
3204 write_pdb(out_name, chain);
3205
3206 ftime(&time1);
3207
3208 if (_VERBOSE) printf("Done. Reconstruction finished in %.3f s.\n", (real)0.001*(1000.*(time1.time-time0.time)+(time1.millitm-time0.millitm)));
3209
3210 return 0;
3211 }
3212