?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / cdist.c
blob09317717279d41039bb2331fbffec632bf97a4a3
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_cdist_c = "$Id$";
30
31 #include <stdlib.h>
32 #include <ctype.h>
33 #include "assert.h"
34 #include "macros.h"
35 #include "vec.h"
36 #include "txtdump.h"
37 #include "cdist.h"
38 #include "futil.h"
39 #include "tpxio.h"
40 #include "rdgroup.h"
41
42 #define NINDEX(ai,aj,natom) (natom*(ai)-((ai)*(ai+1))/2+(aj)-(ai))
43 #define INDEX(ai,aj,natom) ((ai) < (aj))?NINDEX((ai),(aj),natom):NINDEX(aj,ai,natom)
44 #define CHECK(ai,natom) if (((ai)<0) || ((ai)>=natom)) fatal_error(0,"Invalid atom number %d",ai+1)
45
46 static real vdwlen(t_atoms *atoms,int i,int j)
47 {
48 #define NAT 5
49 char anm[NAT] = "HCNOS";
50 /* For time being I give S the same vdw-parameters as O */
51 real dist[NAT][NAT] = {
52 { 2.0, 2.4, 2.4, 2.3, 2.3 },
53 { 2.4, 3.0, 2.9, 2.75, 2.75 },
54 { 2.4, 2.9, 2.7, 2.7, 2.7 },
55 { 2.3, 2.75, 2.7, 2.8, 2.8 },
56 { 2.3, 2.75, 2.7, 2.8, 2.8 }
57 };
58 char ati,atj;
59 int ai,aj;
60
61 ati=toupper((*atoms->atomname[i])[0]);
62 atj=toupper((*atoms->atomname[j])[0]);
63
64 for(ai=0; (ai<NAT) && (ati != anm[ai]); ai++)
65 ;
66 for(aj=0; (aj<NAT) && (atj != anm[aj]); aj++)
67 ;
68 if ((ai < NAT) && (aj < NAT))
69 return dist[ai][aj];
70 else {
71 fprintf(stderr,"No combination for nbs (%c %c) using 2.0 A\n",ati,atj);
72 return 2.0;
73 }
74 }
75
76 void set_dist(t_dist *d,int natoms,int ai,int aj,real lb,
77 real ub,real len)
78 {
79 int index;
80
81 if (ub <= lb)
82 fprintf(stderr,"set_dist: lb = %f, ub = %f, len = %f, atoms %d,%d\n",
83 lb,ub,len,ai,aj);
84 else {
85 CHECK(ai,natoms);
86 CHECK(aj,natoms);
87 index = INDEX(ai,aj,natoms);
88 d[index].lb = lb;
89 d[index].ub = ub;
90 d[index].len = len;
91 }
92 }
93
94 void set_ideal(t_dist *d,int natoms,int ai,int aj,real len)
95 {
96 int index;
97
98 CHECK(ai,natoms);
99 CHECK(aj,natoms);
100 index = INDEX(ai,aj,natoms);
101 d[index].len = len;
102 }
103
104 bool dist_set(t_dist *d,int natoms,int ai,int aj)
105 {
106 int index;
107
108 index = INDEX(ai,aj,natoms);
109
110 return (d[index].lb != 0.0);
111 }
112
113 static t_dist *new_dist(int natom)
114 {
115 t_dist *d;
116
117 snew(d,(natom*(natom+1))/2);
118
119 return d;
120 }
121
122 real d_lb(t_dist *d,int natoms,int ai,int aj)
123 {
124 int index;
125
126 index = INDEX(ai,aj,natoms);
127
128 return d[index].lb;
129 }
130
131 real d_ub(t_dist *d,int natoms,int ai,int aj)
132 {
133 int index;
134
135 index = INDEX(ai,aj,natoms);
136
137 return d[index].ub;
138 }
139
140 real d_len(t_dist *d,int natoms,int ai,int aj)
141 {
142 int index;
143
144 index = INDEX(ai,aj,natoms);
145
146 return d[index].len;
147 }
148
149 static t_dist *read_dist(FILE *log,char *fn,int natom,real weight[])
150 {
151 #define BLEN 255
152 FILE *fp;
153 char buf[BLEN+1];
154 int i,ai,aj,www,ndist;
155 real lb,ub,len;
156 t_dist *d;
157
158 fprintf(log,"Going to read %s\n",fn);
159 d = new_dist(natom);
160 fp = ffopen(fn,"r");
161 ndist = 0;
162 www = 0;
163 while (fgets2(buf,BLEN,fp) != NULL) {
164 if (buf[0] != '#') {
165 if (sscanf(buf,"%d%d%lf%lf%lf",&ai,&aj,&len,&lb,&ub) != 5)
166 if (sscanf(buf,"%d%d%lf%lf",&ai,&aj,&lb,&ub) != 4)
167 fatal_error(0,"Invalid dist format in %s",fn);
168 ai--;
169 aj--;
170 if ((weight[ai] != 0) || (weight[aj] != 0)) {
171 if (!dist_set(d,natom,ai,aj)) {
172 if (debug)
173 fprintf(debug,"\r%d %d %d",ndist,ai,aj);
174 set_dist(d,natom,ai,aj,lb,ub,-1);
175 ndist++;
176 }
177 }
178 else
179 www++;
180 }
181 }
182 fclose(fp);
183
184 fprintf(stderr,
185 "Read %d distances from %s (discarded %d due to zero weight)\n",
186 ndist,fn,www);
187
188 return d;
189 }
190
191 static void measure_report(FILE *fp,int nm,int natom,int nover,int nnotideal,
192 int nhb_rep,int nhb,real cutoff,real rmsd)
193 {
194 double ndtot=(natom*(natom-1.0))/2.0;
195
196 fprintf(fp,"Determined %d distances out of a possible %.0f "
197 "(%.2f %%)\nby measuring from the structure within %f A\n",
198 nm,ndtot,(nm*100.0)/ndtot,cutoff);
199 fprintf(fp,"%d distances from the database were replaced by measured ones\n",
200 nover);
201 fprintf(fp,"Of these, %d were not within the bounds of the database.\n"
202 "RMS deviation from bound limits for these was %e A.\n",
203 nnotideal,rmsd);
204 fprintf(fp,"For %d hydrogen bonds out of %d in the index file the margins "
205 "were reduced\n",nhb_rep/2,nhb/3);
206 }
207
208 static void measure_dist(FILE *log,t_dist *d,t_atoms *atoms,rvec x[],
209 real cutoff,real margin,real hbmargin,
210 int nhb,atom_id hb[])
211 {
212 int i,j,natom,dd,aa,ai,aj,nm,nover,nnotideal,nhb_rep;
213 real msd,rmsd;
214 rvec dx;
215 real ideal,lbfac,ubfac,vdw,lb,ub,nmargin;
216
217 nm = 0;
218 nover = 0;
219 nnotideal = 0;
220 msd = 0;
221 natom = atoms->nr;
222
223 lbfac = 1-margin;
224 ubfac = 1+margin;
225
226 for(ai=0; (ai<natom); ai++)
227 for(aj=ai+1; (aj<natom); aj++) {
228 rvec_sub(x[ai],x[aj],dx);
229 ideal = 10*norm(dx);
230 if (ideal == 0.0) {
231 fprintf(stderr,"Warning distance between atoms %s and %s is zero\n",
232 atomname(atoms,ai),atomname(atoms,aj));
233 }
234 else {
235 if (!dist_set(d,natom,ai,aj)) {
236 /* This distance is not determined by the database. */
237 vdw = 0; /*vdwlen(atoms,ai,aj);*/
238 if ((ideal < cutoff) || (cutoff == 0)) {
239 set_dist(d,natom,ai,aj,max(vdw,ideal*lbfac),ideal*ubfac,ideal);
240 nm++;
241 }
242 }
243 else {
244 /* These distances are already set by the database routines.
245 * However, we override the distances with the measured ones
246 * while keeping the original margins.
247 */
248 lb = d_lb(d,natom,ai,aj);
249 ub = d_ub(d,natom,ai,aj);
250 if ((ideal < lb) || (ideal > ub)) {
251 if (debug)
252 fprintf(debug,"Warning: d(%s,%s) = %8.4f. According to E&H"
253 " it should be %8.4f (dev. %.1f%%)\n",
254 atomname(atoms,ai),
255 atomname(atoms,aj),
256 ideal,(lb+ub)*0.5,50*fabs(ideal*2-lb-ub));
257 msd += (ideal < lb) ? sqr(ideal-lb) : sqr(ideal-ub);
258 nnotideal++;
259 }
260 nmargin = (ub-lb)/(ub+lb);
261 set_dist(d,natom,ai,aj,ideal*(1-nmargin),ideal*(1+nmargin),ideal);
262 nm++;
263 nover++;
264 }
265 }
266 }
267 /* Now we have set all the distances. We have to verify though
268 * that h bonds are maintained, we therefore reduce their margins.
269 */
270 nhb_rep = 0;
271 for(i=0; (i<nhb); i+= 3) {
272 /* Acceptor atom */
273 aa = hb[i+2];
274 assert(aa < natom);
275 for(j=0; (j<2); j++) {
276 /* Donor atom */
277 dd = hb[i+j];
278 assert(dd < natom);
279 if (dist_set(d,natom,dd,aa)) {
280 lb = d_lb(d,natom,dd,aa);
281 ub = d_ub(d,natom,dd,aa);
282 ideal = d_len(d,natom,dd,aa);
283 nmargin = (ub-lb)/(ub+lb);
284 if (hbmargin < nmargin) {
285 set_dist(d,natom,dd,aa,ideal*(1-hbmargin),ideal*(1+hbmargin),ideal);
286 nhb_rep++;
287 }
288 }
289 else
290 fatal_error(0,"Distance between %s and %s not set, while they do make"
291 " a hbond:\nincrease radius for measuring (now %f A)\n",
292 atomname(atoms,aa),atomname(atoms,dd),cutoff);
293 }
294 }
295
296 if (nnotideal > 0)
297 rmsd = sqrt(msd/nnotideal);
298 else
299 rmsd = 0;
300 measure_report(log,nm,natom,nover,nnotideal,nhb_rep,nhb,cutoff,rmsd);
301 measure_report(stderr,nm,natom,nover,nnotideal,nhb_rep,nhb,cutoff,rmsd);
302 }
303
304 static void dump_dist(char *fn,t_dist *d,int natom,bool bAddR)
305 {
306 FILE *fp;
307 int i,j;
308 real lb,ub,len;
309
310 fp = ffopen(fn,"w");
311 fprintf(fp,"# distance file generated by cdist\n");
312
313 for(i=0; (i<natom); i++) {
314 for(j=i+1; (j<natom); j++) {
315 if (dist_set(d,natom,i,j)) {
316 lb = d_lb(d,natom,i,j);
317 ub = d_ub(d,natom,i,j);
318 len = d_len(d,natom,i,j);
319 if (bAddR)
320 fprintf(fp,"%5d%5d%10.5f%10.5f\n",i+1,j+1,lb,ub);
321 else
322 fprintf(fp,"%5d%5d%10.5f%10.5f%10.5f\n",i+1,j+1,len,lb,ub);
323 /* Output for real-precision disco: */
324 /* fprintf(fp,"%5d%5d%15.10f%15.10f%15.10f\n",i+1,j+1,len,lb,ub); */
325 }
326 }
327 }
328 fclose(fp);
329 }
330
331 static real *read_weights(char *fn,int natom)
332 {
333 FILE *in;
334 int i,n;
335 char title[STRLEN];
336 real *w;
337 matrix box;
338 t_atoms atoms;
339 rvec *x;
340
341 /* Open file */
342 in = ffopen(fn,"r");
343
344 /* Check the number of atoms */
345 get_pdb_coordnum(in,&n);
346 if (n != natom)
347 fatal_error(0,"Number of atoms in pdb file (%d) does not match tpx (%d)",
348 n,natom);
349
350 /* Allocate space */
351 snew(w,natom);
352 snew(x,natom);
353 init_t_atoms(&atoms,natom,TRUE);
354 clear_mat(box);
355
356 /* Now read it all */
357 rewind(in);
358 read_pdbfile(in,title,&atoms,x,box,FALSE);
359 fclose(in);
360 fprintf(stderr,"%s\n",title);
361
362 /* Now copy the weights */
363 for(i=0; (i<natom); i++)
364 w[i] = atoms.pdbinfo[i].occup;
365
366 /* Now free temps */
367 sfree(x);
368 free_t_atoms(&atoms);
369
370 /* And go back */
371 return w;
372 }
373
374
375
376
377 void init_rot_matrix(real matrix[3][3],real theta, real omega)
378 {
379 matrix[0][0] = -(cos(theta)*cos(omega));
380 matrix[1][0] = -(cos(theta)*sin(omega));
381 matrix[2][0] = sin(theta);
382
383 matrix[0][1] = sin(omega);
384 matrix[1][1] = -(cos(omega));
385 matrix[2][1] = 0.0;
386
387 matrix[0][2] = sin(theta)*cos(omega);
388 matrix[1][2] = sin(theta)*sin(omega);
389 matrix[2][2] = cos(theta);
390 }
391
392 void vect_matrix(real vect[3], real matrix[3][3])
393 {
394
395 real tmp[3];
396 int i,j;
397
398 tmp[0] = 0.0;
399 tmp[1] = 0.0;
400 tmp[2] = 0.0;
401
402 for ( i=0 ; i <= 2 ; i++ ) {
403 for ( j=0 ; j <=2 ; j++ ) {
404 tmp[i]=tmp[i]+matrix[i][j]*vect[j];
405 }
406 }
407 for ( i=0 ; i <=2 ; i++ ) {
408 vect[i]=tmp[i];
409 }
410 }
411
412 void gauche(int ai,int aj,int ak,int al,t_ilist ilist[],
413 t_iparams iparams[],real *lb,t_atoms *atoms)
414 {
415
416 /* Matrix based approach */
417 int i,j;
418 real matrix1[3][3],matrix2[3][3];
419 real vect1[3],vect2[3],vect3[3];
420 real dist = 0.0;
421 real pi = M_PI;
422 real half_pi = M_PI*0.5;
423 real rij,rjk,rkl;
424 real thijk,thjkl,theta1,theta2,omega1=pi,omega2 = pi+pi/3.0;
425
426 rij = lookup_bondlength(ai,aj,ilist,iparams,TRUE,atoms);
427 rjk = lookup_bondlength(aj,ak,ilist,iparams,TRUE,atoms);
428 rkl = lookup_bondlength(ak,al,ilist,iparams,TRUE,atoms);
429 thijk = lookup_angle(ai,aj,ak,ilist,iparams,atoms);
430 thjkl = lookup_angle(aj,ak,al,ilist,iparams,atoms);
431 theta1 = pi-thijk;
432 theta2 = pi-thjkl;
433
434 /*Initialise vectors*/
435 vect1[0]=0.0;
436 vect1[1]=0.0;
437 vect1[2]=rij;
438 vect2[0]=0.0;
439 vect2[1]=0.0;
440 vect2[2]=rjk;
441 vect3[0]=0.0;
442 vect3[1]=0.0;
443 vect3[2]=rkl;
444
445 /*Initialize matrices*/
446 init_rot_matrix(matrix1,theta1,omega1);
447 init_rot_matrix(matrix2,theta2,omega2);
448
449 /* Express vect2 and vect3 in coord. system of vect1 */
450 vect_matrix(vect2,matrix1);
451 vect_matrix(vect3,matrix2);
452 vect_matrix(vect3,matrix1);
453
454 /* Add vectors */
455 for ( i=0 ; i <=2 ; i++) {
456 vect3[i] = vect1[i]+vect2[i]+vect3[i];
457 }
458
459 /* Calculate distance */
460 *lb = sqrt(vect3[0]*vect3[0]+vect3[1]*vect3[1]+vect3[2]*vect3[2]);
461 }
462
463
464 void gauche15(int ai,int aj,int ak,int al,int am,real omega1,real omega2,
465 real omega3,t_ilist ilist[],t_iparams iparams[],real *lb,
466 t_atoms *atoms)
467 {
468
469 /* Matrix based approach */
470 int i,j;
471 real matrix1[3][3],matrix2[3][3],matrix3[3][3];
472 real vect1[3],vect2[3],vect3[3],vect4[3];
473 real pi = M_PI;
474 real half_pi = M_PI*0.5;
475 real rij,rjk,rkl,rlm;
476 real thijk,thjkl,thklm,theta1,theta2,theta3;
477
478 rij = lookup_bondlength(ai,aj,ilist,iparams,TRUE,atoms);
479 rjk = lookup_bondlength(aj,ak,ilist,iparams,TRUE,atoms);
480 rkl = lookup_bondlength(ak,al,ilist,iparams,TRUE,atoms);
481 rlm = lookup_bondlength(al,am,ilist,iparams,TRUE,atoms);
482 thijk = lookup_angle(ai,aj,ak,ilist,iparams,atoms);
483 thjkl = lookup_angle(aj,ak,al,ilist,iparams,atoms);
484 thklm = lookup_angle(ak,al,am,ilist,iparams,atoms);
485 theta1 = pi-thijk;
486 theta2 = pi-thjkl;
487 theta3 = pi-thklm;
488
489 /*Initialise vectors*/
490 vect1[0]=0.0;
491 vect1[1]=0.0;
492 vect1[2]=rij;
493 vect2[0]=0.0;
494 vect2[1]=0.0;
495 vect2[2]=rjk;
496 vect3[0]=0.0;
497 vect3[1]=0.0;
498 vect3[2]=rkl;
499 vect4[0]=0.0;
500 vect4[1]=0.0;
501 vect4[2]=rlm;
502
503 /*Initialize matrices*/
504 init_rot_matrix(matrix1,theta1,omega1);
505 init_rot_matrix(matrix2,theta2,omega2);
506 init_rot_matrix(matrix3,theta3,omega3);
507
508 /* Express vect2 and vect3 in coord. system of vect1 */
509 vect_matrix(vect2,matrix1);
510 vect_matrix(vect3,matrix2);
511 vect_matrix(vect3,matrix1);
512 vect_matrix(vect4,matrix3);
513 vect_matrix(vect4,matrix2);
514 vect_matrix(vect4,matrix1);
515
516 /* Add vectors */
517 for ( i=0 ; i <=2 ; i++) {
518 vect4[i] = vect1[i]+vect2[i]+vect3[i]+vect4[i];
519 }
520
521 /* Calculate distance */
522 *lb = sqrt(vect4[0]*vect4[0]+vect4[1]*vect4[1]+vect4[2]*vect4[2]);
523 }
524
525
526 static void dump_bonds(t_atoms *atoms,t_dist *d,
527 t_ilist ilist[],t_functype functype[],
528 t_iparams ip[],real bond_margin,real angle_margin,
529 real dih_margin,real idih_margin,
530 real weight[],bool bAddR)
531 {
532 int dodist[] = { F_BONDS, F_MORSE, F_SHAKE, F_G96BONDS,
533 F_G96ANGLES, F_ANGLES, F_IDIHS /*,F_RBDIHS, F_PDIHS*/ };
534 int i,j,atom1,atom2,ai,aj,ak,al,type,ftype,nra,ndist,odist,natoms;
535 real blen,rij,rjk,c6,c12,lb=-1,ub,theta;
536
537 natoms= atoms->nr;
538 ndist = 0;
539 odist = 0;
540
541 for(j=0; (j<asize(dodist)); j++) {
542 ftype = dodist[j];
543 nra = interaction_function[ftype].nratoms;
544 for(i=0; (i<ilist[ftype].nr); i+=nra+1) {
545 type = ilist[ftype].iatoms[i];
546 ai = ilist[ftype].iatoms[i+1];
547 aj = ilist[ftype].iatoms[i+2];
548 atom1 = ai;
549 atom2 = aj;
550 blen = 0;
551 switch (ftype) {
552 case F_BONDS:
553 case F_MORSE:
554 case F_SHAKE:
555 blen = lookup_bondlength(ai,aj,ilist,ip,TRUE,atoms);
556 break;
557 case F_G96ANGLES:
558 case F_ANGLES:
559 ak = ilist[ftype].iatoms[i+3];
560 theta = lookup_angle(ai,aj,ak,ilist,ip,atoms);
561 blen = angle_length(ai,aj,ak,RAD2DEG*theta,ilist,ip,atoms);
562 atom2 = ak;
563 break;
564 case F_PDIHS:
565 case F_RBDIHS:
566 ak = ilist[ftype].iatoms[i+3];
567 al = ilist[ftype].iatoms[i+4];
568 pdih_lengths(ai,aj,ak,al,ilist,ip,&blen,&ub,atoms);
569 atom2 = al;
570 break;
571 case F_IDIHS:
572 ak = ilist[ftype].iatoms[i+3];
573 al = ilist[ftype].iatoms[i+4];
574 blen = idih_lengths(ai,aj,ak,al,ilist,ip,atoms);
575 atom2 = al;
576 break;
577 default:
578 break;
579 }
580 if ((blen != 0) && ((weight[atom1] != 0.0) || (weight[atom2] != 0.0))) {
581 switch (ftype) {
582 case F_BONDS:
583 case F_MORSE:
584 case F_SHAKE:
585 case F_G96BONDS:
586 lb = (1.0-bond_margin)*blen;
587 ub = (1.0+bond_margin)*blen;
588 break;
589 case F_ANGLES:
590 case F_G96ANGLES:
591 lb = (1.0-angle_margin)*blen;
592 ub = (1.0+angle_margin)*blen;
593 break;
594 case F_IDIHS:
595 lb = (1-idih_margin)*blen;
596 ub = (1+idih_margin)*blen;
597 break;
598 case F_PDIHS:
599 case F_RBDIHS:
600 lb = (1.0-dih_margin)*blen;
601 ub = (1.0+dih_margin)*ub;
602 break;
603 }
604 if (!dist_set(d,natoms,atom1,atom2)) {
605 set_dist(d,natoms,atom1,atom2,lb,ub,blen);
606 ndist ++;
607 }
608 else
609 odist ++;
610 }
611 }
612 }
613 fprintf(stderr,"There are %d new bonded distances + %d old\n",ndist,odist);
614 }
615
616 static bool notexcl(t_block *excl,int ai,int aj)
617 {
618 int i;
619
620 for(i=excl->index[ai]; (i<excl->index[ai+1]); i++) {
621 if (aj == excl->a[i])
622 return FALSE;
623 }
624 return TRUE;
625 }
626
627 static void dump_nonbonds(t_dist *d,t_idef *idef,t_atoms *atoms,
628 real hblen,real weight[],
629 real nb_margin,bool bAddR)
630 {
631 int i,j,k,natoms,ntype,tpi,tpj,ndist,odist;
632 real **dmat,len,lb,maxdist;
633 bool *bDA;
634 char element,ati,atj;
635
636 ntype = idef->atnr;
637 natoms = atoms->nr;
638 maxdist = (atoms->nres+3)*3.5; /* Ca-Ca distance + some buffer */
639
640 /* Determine which atoms are capable of H bonding or not */
641 snew(bDA,ntype);
642 for(j=0; (j<ntype); j++)
643 for(i=0; (i<natoms); i++) {
644 if (atoms->atom[i].type == j) {
645 element = (*atoms->atomname[i])[0];
646 bDA[j] = ((element == 'N') || (element == 'O'));
647 break;
648 }
649 }
650
651 /* Make a distance matrix for VDWaals distances */
652 /*snew(dmat,ntype);
653 for(i=0; (i<ntype); i++) {
654 snew(dmat[i],ntype);
655 for(j=0; (j<ntype); j++) {
656 if (bDA[i] && bDA[j])
657 dmat[i][j] = hblen;
658 else
659 dmat[i][j] = vdwlen(idef,i,j,hblen);
660 }
661 } */
662 sfree(bDA);
663
664 /* Now loop over atom pairs */
665 ndist = 0;
666 odist = 0;
667 for(i=0; (i<natoms); i++) {
668 tpi = atoms->atom[i].type;
669 for(j=i+1; (j<natoms); j++) {
670 if (((weight[i] != 0.0) || (weight[j] != 0.0))
671 /*&& (notexcl(&atoms->excl,i,j))*/) {
672 if (!dist_set(d,natoms,i,j)) {
673 /* We don't care about virtual particles */
674 ati=(*atoms->atomname[i])[0];
675 atj=(*atoms->atomname[j])[0];
676 if ( !(ati == 'V' || atj == 'V') ) {
677 /* *** */
678 tpj = atoms->atom[j].type;
679 len = vdwlen(atoms,i,j);
680 lb = (1-nb_margin)*len;
681 if (len > 0) {
682 set_dist(d,natoms,i,j,lb,maxdist,0.0);
683 /* set_dist(d,natoms,i,j,lb,0.0,lb); */
684 ndist++;
685 }
686 }
687 }
688 else
689 odist++;
690 }
691 }
692 }
693 fprintf(stderr,"There are %d new non-bonded distances + %d old ones\n",
694 ndist,odist);
695 }
696
697
698 int main(int argc,char *argv[])
699 {
700 static char *desc[] = {
701 "cdist read a [BB]tpx[bb] file and dumps an input file for disco.",
702 "Bond lengths etc. are read from the topology. Pairs of atoms that can",
703 "form hydrogen bonds are given a lowest possible distance of",
704 "[BB]hblen[bb] (can be specified by the user). Other nonbonded pairs",
705 "take their minimum distance from the Lennard Jones parameters",
706 "(at the combined sigma).[PAR]",
707 "The program uses proper dihedrals to give a distance too, as minimum",
708 "respectively maximum the [IT]cis[it] and [IT]trans[it] configurations",
709 "are taken. It is therefore beneficial to use the [BB]-alldih[bb] option",
710 "of [TT]pdb2gmx[tt] to generate a topology with all dihedrals in there.",
711 "If the optional pdb file is given, weights are read from the occupancy",
712 "field, so that",
713 "not all atoms are part of the disco run, only those of which one of the",
714 "weights is non-zero.[PAR]",
715 "If the option -engh is on (default) bond lengths and angles etc. are",
716 "read from another database, which is basically the Engh-Huber data",
717 "but refined to be completely self consistent. The database name is",
718 "refi_aa.dat and it resides in the $GMXLIB directory, or in the current",
719 "directory.[PAR]",
720 "The program can read a file with distances from NMR distance restraints",
721 "(-d option). Note that these distance are treated slightly different",
722 "in the disco program, and therefore these distance should be NMR",
723 "derived distance restraints only.[PAR]",
724 "Furthermore, the program can read an index file with hydrogen bond",
725 "information as generated by [TT]g_hbond[tt]. This is then used to set",
726 "tighter restraints on the hydrogen bonded atoms than on the other",
727 "non bonded atom pairs, in order to maintain secondary structure.",
728 "This option is useful only in combination with the [TT]-measure[tt]",
729 "option, when a sensible structure is known."
730 };
731 t_filenm fnm[] = {
732 { efTPS, "-s", NULL, ffREAD },
733 { efLOG, "-g", "cdist", ffWRITE },
734 { efPDB, "-q", NULL, ffOPTRD },
735 { efDAT, "-d", NULL, ffOPTRD },
736 { efDAT, "-o", "cdist", ffWRITE },
737 { efNDX, "-n", "hbond", ffOPTRD }
738 };
739 #define NFILE asize(fnm)
740
741 FILE *fp;
742 t_topology *top;
743 t_dist *dist;
744 real *weight;
745 rvec *x;
746 matrix box;
747 char *topfn,title[256];
748 int i,nhb;
749 atom_id *hb;
750 char *grpname;
751
752 /* Tolerance used in smoothing functions (real precision)*/
753
754 /* Hacked 990609 by Adam */
755
756 /* Command line options */
757 static real tol=1e-6;
758 static real bond_margin = 0.01;
759 static real angle_margin = 0.01;
760 /* static real pep_margin = 0.005; */
761 static real pep_margin = 0.01;
762 /* static real ring_margin = 0.002; */
763 static real ring_margin = 0.01;
764 /* static real arg_margin = 0.002; */
765 static real arg_margin = 0.01;
766 /* Use end_margin for asn and gln */
767 /* static real end_margin = 0.004; */
768 static real end_margin = 0.01;
769 static real val_margin = 0.01;
770 static real leu_margin = 0.01;
771 static real ile_margin = 0.03;
772 static real dih_margin = 0.01;
773 static real idih_margin = 0.01;
774 static real nb_margin = 0.05;
775 static real hb_margin = 0.02;
776 static real hblen = 2.3;
777 static real rcut = 0.0;
778 static bool bNB=TRUE,bBON=TRUE,bAddR=FALSE;
779 static bool bVir=FALSE,bEnghHuber=TRUE;
780 static char *smth_str[] = { NULL, "none", "tri", "tetra", NULL };
781 t_pargs pa[] = {
782 { "-engh",FALSE,etBOOL, &bEnghHuber,
783 "Use the Engh&Huber parameters for bond-lengths etc." },
784 { "-tol", FALSE, etREAL, &tol,
785 "HIDDENTolerance for smoothing" },
786 { "-bm", FALSE, etREAL, &bond_margin,
787 "Relative margin for bond lengths" },
788 { "-am", FALSE, etREAL, &angle_margin,
789 "Relative margin for bond angle lengths" },
790 { "-pm", FALSE, etREAL, &pep_margin,
791 "Relative margin for peptidebond dihedrals" },
792 { "-rr", FALSE, etREAL, &ring_margin,
793 "Relative margin to keep rings flat (trp,tyr,phe,hisb)" },
794 { "-ar", FALSE, etREAL, &arg_margin,
795 "Relative margin for arginine" },
796 { "-er", FALSE, etREAL, &end_margin,
797 "Relative margin for asn and gln" },
798 { "-vm", FALSE, etREAL, &val_margin,
799 "Relative margin for valine (0 disables)" },
800 { "-lm", FALSE, etREAL, &leu_margin,
801 "Relative margin for leucine (0 disables)" },
802 { "-il", FALSE, etREAL, &ile_margin,
803 "Relative margin for isoleucine (0 disables)" },
804 { "-dm", FALSE, etREAL, &dih_margin,
805 "!inactive! Relative margin for dihedral lengths" },
806 { "-im", FALSE, etREAL, &idih_margin,
807 "Relative margin for improper dihedral lengths" },
808 { "-nm", FALSE, etREAL, &nb_margin,
809 "Relative margin for nonbonded lower bounds" },
810 { "-hm", FALSE, etREAL, &hb_margin,
811 "Relative margin for hydrogen bonded atoms, which must be specified in an index file, as generated by g_hbond" },
812 { "-hb", FALSE, etREAL, &hblen,
813 "Shortest possible distance for a hydrogen bond (in Angstrom!)" },
814 { "-bon", FALSE, etBOOL, &bBON,
815 "Make bonded distance constraints" },
816 { "-nb", FALSE, etBOOL, &bNB,
817 "Make nonbonded distance constraints (lower bound only) " },
818 { "-measure", FALSE, etREAL, &rcut,
819 "Add (nonbonded) distances by examining all atoms within the distance given (in Angstrom), and using the margin given by the -nm option." },
820 { "-add",FALSE, etBOOL, &bAddR,
821 "Write restraints in format of additional restraints for disco" },
822 { "-vir",FALSE, etBOOL, &bVir,
823 "Use virtual particles"},
824 { "-sm", FALSE, etENUM, smth_str,
825 "Smoothing: none, tri (Using triangle inequality), or tetra (Partial tetrangle inequaliy)" }
826 };
827
828 CopyRight(stdout,argv[0]);
829 parse_common_args(&argc,argv,0,FALSE,NFILE,fnm,asize(pa),pa,
830 asize(desc),desc,0,NULL);
831
832 #ifndef DOUBLE
833 fprintf(stderr,"WARNING: running %s in single precision may produce bad"
834 " distances\n",argv[0]);
835 #endif
836
837 fp=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
838
839 if ((!bBON && !bNB) && (rcut==0.0))
840 fprintf(stderr,"That was interesting... (nothing done)\n");
841 else {
842 snew(top,1);
843 topfn = ftp2fn(efTPS,NFILE,fnm);
844 if (!read_tps_conf(topfn,title,top,&x,NULL,box,FALSE))
845 fatal_error(0,"No topology in %s",topfn);
846 fprintf(stderr,"Successfully read %s (%s)\n",topfn,title);
847
848 if (opt2bSet("-q",NFILE,fnm)) {
849 weight = read_weights(opt2fn("-q",NFILE,fnm),top->atoms.nr);
850 }
851 else {
852 snew(weight,top->atoms.nr);
853 for(i=0; (i<top->atoms.nr); i++)
854 weight[i] = 1.0;
855 }
856 if (opt2bSet("-d",NFILE,fnm)) {
857 dist = read_dist(fp,opt2fn("-d",NFILE,fnm),top->atoms.nr,weight);
858 }
859 else
860 dist = new_dist(top->atoms.nr);
861
862 if (bEnghHuber)
863 read_O_dist();
864
865 if (bBON) {
866 peptide_bonds(fp,dist,&top->idef,&top->atoms,weight,pep_margin,
867 top->idef.il,top->idef.iparams,bVir);
868 arg(fp,dist,&top->idef,&top->atoms,weight,arg_margin,
869 top->idef.il,top->idef.iparams,bVir);
870 asn(fp,dist,&top->idef,&top->atoms,weight,end_margin,
871 top->idef.il,top->idef.iparams,bVir);
872 gln(fp,dist,&top->idef,&top->atoms,weight,end_margin,
873 top->idef.il,top->idef.iparams,bVir);
874 phe(fp,dist,&top->idef,&top->atoms,weight,ring_margin,
875 top->idef.il,top->idef.iparams,bVir);
876 tyr(fp,dist,&top->idef,&top->atoms,weight,ring_margin,
877 top->idef.il,top->idef.iparams,bVir);
878 trp(fp,dist,&top->idef,&top->atoms,weight,ring_margin,
879 top->idef.il,top->idef.iparams,bVir);
880 hisb(fp,dist,&top->idef,&top->atoms,weight,ring_margin,
881 top->idef.il,top->idef.iparams,bVir);
882 if ( val_margin != 0 ) {
883 val(fp,dist,&top->idef,&top->atoms,weight,val_margin,
884 top->idef.il,top->idef.iparams);
885 }
886 if ( leu_margin != 0 ) {
887 leu(fp,dist,&top->idef,&top->atoms,weight,leu_margin,
888 top->idef.il,top->idef.iparams);
889 }
890 if ( ile_margin != 0 ) {
891 ile(fp,dist,&top->idef,&top->atoms,weight,ile_margin,
892 top->idef.il,top->idef.iparams);
893 }
894 fprintf(stderr,"Done residues...\n");
895 dump_bonds(&top->atoms,dist,
896 top->idef.il,top->idef.functype,top->idef.iparams,
897 bond_margin,angle_margin,dih_margin,
898 idih_margin,weight,bAddR);
899 }
900 fprintf(stderr,"Done bondeds\n");
901
902 if (rcut > 0) {
903 nhb = 0;
904 hb = NULL;
905 grpname = NULL;
906 if (ftp2bSet(efNDX,NFILE,fnm)) {
907 rd_index(ftp2fn(efNDX,NFILE,fnm),1,&nhb,&hb,&grpname);
908 fprintf(stderr,"There are %d hydrogen bonds\n",nhb/3);
909 }
910 measure_dist(fp,dist,&top->atoms,x,rcut,nb_margin,hb_margin,nhb,hb);
911 }
912 if (bNB)
913 dump_nonbonds(dist,&top->idef,&top->atoms,hblen,weight,
914 nb_margin,bAddR);
915
916 if (strcmp(smth_str[0],smth_str[1]) == 0)
917 fprintf(stderr,"No smoothing\n");
918 else if (strcmp(smth_str[0],smth_str[2]) == 0) {
919 fprintf(stderr,"Triangle smoothing only\n");
920 do_triangle (dist,&top->atoms,tol);
921 }
922 else if (strcmp(smth_str[0],smth_str[3]) == 0) {
923 fprintf(stderr,"Partial tetrangle + triangle smoothing\n");
924 do_smooth(dist,&top->atoms,tol);
925 }
926 else
927 fatal_error(0,"Uh-oh, smth_str = %s, %s, line %d",
928 smth_str[0],__FILE__,__LINE__);
929
930 dump_dist(opt2fn("-o",NFILE,fnm),dist,top->atoms.nr,bAddR);
931 }
932
933 ffclose(fp);
934
935 thanx(stdout);
936
937 return 0;
938 }
got merge a long time ago