?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Parallel vs. sequentiual code: I get binary identical result (apart from comment...
[gromacs/qmmm-gamess-us.git] / src / tools / gmx_chi.c
blob35c40a253ebd187f9c9d4e617886b2e83c53b8c2
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
26 *
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38 #include <stdio.h>
39 #include <math.h>
40
41 #include "confio.h"
42 #include "pdbio.h"
43 #include "copyrite.h"
44 #include "gmx_fatal.h"
45 #include "futil.h"
46 #include "gstat.h"
47 #include "macros.h"
48 #include "maths.h"
49 #include "physics.h"
50 #include "index.h"
51 #include "smalloc.h"
52 #include "statutil.h"
53 #include "tpxio.h"
54 #include "string.h"
55 #include "sysstuff.h"
56 #include "txtdump.h"
57 #include "typedefs.h"
58 #include "vec.h"
59 #include "strdb.h"
60 #include "xvgr.h"
61 #include "matio.h"
62 #include "gmx_ana.h"
63
64 static bool bAllowed(real phi,real psi)
65 {
66 static const char *map[] = {
67 "1100000000000000001111111000000000001111111111111111111111111",
68 "1100000000000000001111110000000000011111111111111111111111111",
69 "1100000000000000001111110000000000011111111111111111111111111",
70 "1100000000000000001111100000000000111111111111111111111111111",
71 "1100000000000000001111100000000000111111111111111111111111111",
72 "1100000000000000001111100000000001111111111111111111111111111",
73 "1100000000000000001111100000000001111111111111111111111111111",
74 "1100000000000000001111100000000011111111111111111111111111111",
75 "1110000000000000001111110000000111111111111111111111111111111",
76 "1110000000000000001111110000001111111111111111111111111111111",
77 "1110000000000000001111111000011111111111111111111111111111111",
78 "1110000000000000001111111100111111111111111111111111111111111",
79 "1110000000000000001111111111111111111111111111111111111111111",
80 "1110000000000000001111111111111111111111111111111111111111111",
81 "1110000000000000001111111111111111111111111111111111111111111",
82 "1110000000000000001111111111111111111111111111111111111111111",
83 "1110000000000000001111111111111110011111111111111111111111111",
84 "1110000000000000001111111111111100000111111111111111111111111",
85 "1110000000000000001111111111111000000000001111111111111111111",
86 "1100000000000000001111111111110000000000000011111111111111111",
87 "1100000000000000001111111111100000000000000011111111111111111",
88 "1000000000000000001111111111000000000000000001111111111111110",
89 "0000000000000000001111111110000000000000000000111111111111100",
90 "0000000000000000000000000000000000000000000000000000000000000",
91 "0000000000000000000000000000000000000000000000000000000000000",
92 "0000000000000000000000000000000000000000000000000000000000000",
93 "0000000000000000000000000000000000000000000000000000000000000",
94 "0000000000000000000000000000000000000000000000000000000000000",
95 "0000000000000000000000000000000000000000000000000000000000000",
96 "0000000000000000000000000000000000000000000000000000000000000",
97 "0000000000000000000000000000000000000000000000000000000000000",
98 "0000000000000000000000000000000000000000000000000000000000000",
99 "0000000000000000000000000000000000000000000000000000000000000",
100 "0000000000000000000000000000000000000000000000000000000000000",
101 "0000000000000000000000000000000000000000000000000000000000000",
102 "0000000000000000000000000000000000000000000000000000000000000",
103 "0000000000000000000000000000000000000000000000000000000000000",
104 "0000000000000000000000000000000000000000000000000000000000000",
105 "0000000000000000000000000000000000111111111111000000000000000",
106 "1100000000000000000000000000000001111111111111100000000000111",
107 "1100000000000000000000000000000001111111111111110000000000111",
108 "0000000000000000000000000000000000000000000000000000000000000",
109 "0000000000000000000000000000000000000000000000000000000000000",
110 "0000000000000000000000000000000000000000000000000000000000000",
111 "0000000000000000000000000000000000000000000000000000000000000",
112 "0000000000000000000000000000000000000000000000000000000000000",
113 "0000000000000000000000000000000000000000000000000000000000000",
114 "0000000000000000000000000000000000000000000000000000000000000",
115 "0000000000000000000000000000000000000000000000000000000000000",
116 "0000000000000000000000000000000000000000000000000000000000000",
117 "0000000000000000000000000000000000000000000000000000000000000",
118 "0000000000000000000000000000000000000000000000000000000000000",
119 "0000000000000000000000000000000000000000000000000000000000000",
120 "0000000000000000000000000000000000000000000000000000000000000",
121 "0000000000000000000000000000000000000000000000000000000000000",
122 "0000000000000000000000000000000000000000000000000000000000000",
123 "0000000000000000000000000000000000000000000000000000000000000",
124 "0000000000000000000000000000000000000000000000000000000000000",
125 "0000000000000000000000000000000000000000000000000000000000000",
126 "0000000000000000000000000000000000000000000000000000000000000",
127 "0000000000000000000000000000000000000000000000000000000000000"
128 };
129 #define NPP asize(map)
130 int x,y;
131
132 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
133 x = INDEX(phi);
134 y = INDEX(psi);
135 #undef INDEX
136 return (bool) map[x][y];
137 }
138
139 atom_id *make_chi_ind(int nl,t_dlist dl[],int *ndih)
140 {
141 atom_id *id;
142 int i,Xi,n;
143
144 /* There are nl residues with max edMax dihedrals with 4 atoms each */
145 snew(id,nl*edMax*4);
146
147 n=0;
148 for(i=0; (i<nl); i++)
149 {
150 /* Phi, fake the first one */
151 dl[i].j0[edPhi] = n/4;
152 if(dl[i].atm.minC >= 0)
153 id[n++]=dl[i].atm.minC;
154 else
155 id[n++]=dl[i].atm.H;
156 id[n++]=dl[i].atm.N;
157 id[n++]=dl[i].atm.Cn[1];
158 id[n++]=dl[i].atm.C;
159 }
160 for(i=0; (i<nl); i++)
161 {
162 /* Psi, fake the last one */
163 dl[i].j0[edPsi] = n/4;
164 id[n++]=dl[i].atm.N;
165 id[n++]=dl[i].atm.Cn[1];
166 id[n++]=dl[i].atm.C;
167 if ( i< (nl-1) )
168 id[n++]=dl[i+1].atm.N;
169 else
170 id[n++]=dl[i].atm.O;
171 }
172 for(i=0; (i<nl); i++)
173 {
174 /* Omega */
175 if (has_dihedral(edOmega,&(dl[i])))
176 {
177 dl[i].j0[edOmega] = n/4;
178 id[n++]=dl[i].atm.minO;
179 id[n++]=dl[i].atm.minC;
180 id[n++]=dl[i].atm.N;
181 id[n++]=dl[i].atm.H;
182 }
183 }
184 for(Xi=0; (Xi<MAXCHI); Xi++)
185 {
186 /* Chi# */
187 for(i=0; (i<nl); i++)
188 {
189 if (dl[i].atm.Cn[Xi+3] != -1)
190 {
191 dl[i].j0[edChi1+Xi] = n/4;
192 id[n++]=dl[i].atm.Cn[Xi];
193 id[n++]=dl[i].atm.Cn[Xi+1];
194 id[n++]=dl[i].atm.Cn[Xi+2];
195 id[n++]=dl[i].atm.Cn[Xi+3];
196 }
197 }
198 }
199 *ndih=n/4;
200
201 return id;
202 }
203
204 int bin(real chi,int mult)
205 {
206 mult=3;
207
208 return (int) (chi*mult/360.0);
209 }
210
211
212 static void do_dihcorr(const char *fn,int nf,int ndih,real **dih,real dt,
213 int nlist,t_dlist dlist[],real time[],int maxchi,
214 bool bPhi,bool bPsi,bool bChi,bool bOmega,
215 const output_env_t oenv)
216 {
217 char name1[256],name2[256];
218 int i,j,Xi;
219
220 do_autocorr(fn,oenv,"Dihedral Autocorrelation Function",
221 nf,ndih,dih,dt,eacCos,FALSE);
222 /* Dump em all */
223 j=0;
224 for(i=0; (i<nlist); i++) {
225 if (bPhi)
226 print_one(oenv,"corrphi",dlist[i].name,"Phi ACF for", "C(t)", nf/2,time,
227 dih[j]);
228 j++;
229 }
230 for(i=0; (i<nlist); i++) {
231 if (bPsi)
232 print_one(oenv,"corrpsi",dlist[i].name,"Psi ACF for","C(t)",nf/2,time,
233 dih[j]);
234 j++;
235 }
236 for(i=0; (i<nlist); i++) {
237 if (has_dihedral(edOmega,&dlist[i])) {
238 if (bOmega)
239 print_one(oenv,"corromega",dlist[i].name,"Omega ACF for","C(t)",
240 nf/2,time,dih[j]);
241 j++;
242 }
243 }
244 for(Xi=0; (Xi<maxchi); Xi++) {
245 sprintf(name1, "corrchi%d", Xi+1);
246 sprintf(name2, "Chi%d ACF for", Xi+1);
247 for(i=0; (i<nlist); i++) {
248 if (dlist[i].atm.Cn[Xi+3] != -1) {
249 if (bChi)
250 print_one(oenv,name1,dlist[i].name,name2,"C(t)",nf/2,time,dih[j]);
251 j++;
252 }
253 }
254 }
255 fprintf(stderr,"\n");
256 }
257
258 static void copy_dih_data(real in[], real out[], int nf, bool bLEAVE)
259 {
260 /* if bLEAVE, do nothing to data in copying to out
261 * otherwise multiply by 180/pi to convert rad to deg */
262 int i ;
263 real mult ;
264 if (bLEAVE)
265 mult = 1 ;
266 else
267 mult = (180.0/M_PI);
268 for (i=0;(i<nf);i++){
269 out[i]=in[i]*mult ;
270 }
271 }
272
273 static void dump_em_all(int nlist,t_dlist dlist[],int nf,real time[],
274 real **dih,int maxchi,
275 bool bPhi,bool bPsi,bool bChi,bool bOmega, bool bRAD,
276 const output_env_t oenv)
277 {
278 char name[256], titlestr[256], ystr[256];
279 real *data ;
280 int i,j,Xi;
281
282 snew(data,nf);
283 if (bRAD)
284 strcpy(ystr,"Angle (rad)");
285 else
286 strcpy(ystr,"Angle (degrees)");
287
288 /* Dump em all */
289 j = 0;
290 for(i=0; (i<nlist); i++) {
291 /* grs debug printf("OK i %d j %d\n", i, j) ; */
292 if (bPhi) {
293 copy_dih_data(dih[j],data,nf,bRAD);
294 print_one(oenv,"phi",dlist[i].name,"\\xf\\f{}",ystr, nf,time,data);
295 }
296 j++;
297 }
298 for(i=0; (i<nlist); i++) {
299 if (bPsi) {
300 copy_dih_data(dih[j],data,nf,bRAD);
301 print_one(oenv,"psi",dlist[i].name,"\\xy\\f{}",ystr, nf,time,data);
302 }
303 j++;
304 }
305 for(i=0; (i<nlist); i++)
306 if (has_dihedral(edOmega,&(dlist[i]))) {
307 if (bOmega){
308 copy_dih_data(dih[j],data,nf,bRAD);
309 print_one(oenv,"omega",dlist[i].name,"\\xw\\f{}",ystr,nf,time,data);
310 }
311 j++;
312 }
313
314 for(Xi=0; (Xi<maxchi); Xi++)
315 for(i=0; (i<nlist); i++)
316 if (dlist[i].atm.Cn[Xi+3] != -1) {
317 if (bChi) {
318 sprintf(name,"chi%d",Xi+1);
319 sprintf(titlestr,"\\xc\\f{}\\s%d\\N",Xi+1);
320 copy_dih_data(dih[j],data,nf,bRAD);
321 print_one(oenv,name,dlist[i].name,titlestr,ystr, nf,time,data);
322 }
323 j++;
324 }
325 fprintf(stderr,"\n");
326 }
327
328 static void reset_one(real dih[],int nf,real phase)
329 {
330 int j;
331
332 for(j=0; (j<nf); j++) {
333 dih[j] += phase;
334 while(dih[j] < -M_PI)
335 dih[j] += 2*M_PI;
336 while(dih[j] >= M_PI)
337 dih[j] -= 2*M_PI;
338 }
339 }
340
341 static int reset_em_all(int nlist,t_dlist dlist[],int nf,
342 real **dih,int maxchi)
343 {
344 int i,j,Xi;
345
346 /* Reset em all */
347 j=0;
348 /* Phi */
349 for(i=0; (i<nlist); i++)
350 {
351 if (dlist[i].atm.minC == -1)
352 {
353 reset_one(dih[j++],nf,M_PI);
354 }
355 else
356 {
357 reset_one(dih[j++],nf,0);
358 }
359 }
360 /* Psi */
361 for(i=0; (i<nlist-1); i++)
362 {
363 reset_one(dih[j++],nf,0);
364 }
365 /* last Psi is faked from O */
366 reset_one(dih[j++],nf,M_PI);
367
368 /* Omega */
369 for(i=0; (i<nlist); i++)
370 if (has_dihedral(edOmega,&dlist[i]))
371 reset_one(dih[j++],nf,0);
372 /* Chi 1 thru maxchi */
373 for(Xi=0; (Xi<maxchi); Xi++)
374 {
375 for(i=0; (i<nlist); i++)
376 {
377 if (dlist[i].atm.Cn[Xi+3] != -1)
378 {
379 reset_one(dih[j],nf,0);
380 j++;
381 }
382 }
383 }
384 fprintf(stderr,"j after resetting (nr. active dihedrals) = %d\n",j);
385 return j ;
386 }
387
388 static void histogramming(FILE *log,int nbin, int naa,char **aa,
389 int nf,int maxchi,real **dih,
390 int nlist,t_dlist dlist[],
391 atom_id index[],
392 bool bPhi,bool bPsi,bool bOmega,bool bChi,
393 bool bNormalize,bool bSSHisto,const char *ssdump,
394 real bfac_max,t_atoms *atoms,
395 bool bDo_jc, const char *fn,
396 const output_env_t oenv)
397 {
398 /* also gets 3J couplings and order parameters S2 */
399 t_karplus kkkphi[] = {
400 { "J_NHa1", 6.51, -1.76, 1.6, -M_PI/3, 0.0, 0.0 },
401 { "J_NHa2", 6.51, -1.76, 1.6, M_PI/3, 0.0, 0.0 },
402 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0, 0.0 },
403 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0, 0.0 },
404 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0, 0.0 }
405 };
406 t_karplus kkkpsi[] = {
407 { "J_HaN", -0.88, -0.61,-0.27,M_PI/3, 0.0, 0.0 }
408 };
409 t_karplus kkkchi1[] = {
410 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0, 0.0 },
411 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0, 0.0 }
412 };
413 #define NKKKPHI asize(kkkphi)
414 #define NKKKPSI asize(kkkpsi)
415 #define NKKKCHI asize(kkkchi1)
416 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
417
418 FILE *fp,*ssfp[3]={NULL,NULL,NULL};
419 const char *sss[3] = { "sheet", "helix", "coil" };
420 real S2;
421 real *normhisto;
422 real **Jc,**Jcsig;
423 int ****his_aa_ss=NULL;
424 int ***his_aa,**his_aa1,*histmp;
425 int i,j,k,m,n,nn,Dih,nres,hindex,angle;
426 bool bBfac,bOccup;
427 char hisfile[256],hhisfile[256],sshisfile[256],title[256],*ss_str=NULL;
428 char **leg;
429
430 if (bSSHisto) {
431 fp = ffopen(ssdump,"r");
432 if(1 != fscanf(fp,"%d",&nres))
433 {
434 gmx_fatal(FARGS,"Error reading from file %s",ssdump);
435 }
436
437 snew(ss_str,nres+1);
438 if( 1 != fscanf(fp,"%s",ss_str))
439 {
440 gmx_fatal(FARGS,"Error reading from file %s",ssdump);
441 }
442
443 ffclose(fp);
444 /* Four dimensional array... Very cool */
445 snew(his_aa_ss,3);
446 for(i=0; (i<3); i++) {
447 snew(his_aa_ss[i],naa+1);
448 for(j=0; (j<=naa); j++) {
449 snew(his_aa_ss[i][j],edMax);
450 for(Dih=0; (Dih<edMax); Dih++)
451 snew(his_aa_ss[i][j][Dih],nbin+1);
452 }
453 }
454 }
455 snew(his_aa,edMax);
456 for(Dih=0; (Dih<edMax); Dih++) {
457 snew(his_aa[Dih],naa+1);
458 for(i=0; (i<=naa); i++) {
459 snew(his_aa[Dih][i],nbin+1);
460 }
461 }
462 snew(histmp,nbin);
463
464 snew(Jc,nlist);
465 snew(Jcsig,nlist);
466 for(i=0; (i<nlist); i++) {
467 snew(Jc[i],NJC);
468 snew(Jcsig[i],NJC);
469 }
470
471 j=0;
472 n=0;
473 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
474 for(i=0; (i<nlist); i++) {
475 if (((Dih < edOmega) ) ||
476 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
477 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
478 make_histo(log,nf,dih[j],nbin,histmp,-M_PI,M_PI);
479
480 if (bSSHisto) {
481 /* Assume there is only one structure, the first.
482 * Compute index in histogram.
483 */
484 /* Check the atoms to see whether their B-factors are low enough
485 * Check atoms to see their occupancy is 1.
486 */
487 bBfac = bOccup = TRUE;
488 for(nn=0; (nn<4); nn++,n++) {
489 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
490 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
491 }
492 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac))) {
493 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
494 range_check(hindex,0,nbin);
495
496 /* Assign dihedral to either of the structure determined
497 * histograms
498 */
499 switch(ss_str[dlist[i].resnr]) {
500 case 'E':
501 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
502 break;
503 case 'H':
504 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
505 break;
506 default:
507 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
508 break;
509 }
510 }
511 else if (debug)
512 fprintf(debug,"Res. %d has imcomplete occupancy or bfacs > %g\n",
513 dlist[i].resnr,bfac_max);
514 }
515 else
516 n += 4;
517
518 switch (Dih) {
519 case edPhi:
520 calc_distribution_props(nbin,histmp,-M_PI,NKKKPHI,kkkphi,&S2);
521
522 for(m=0; (m<NKKKPHI); m++) {
523 Jc[i][m] = kkkphi[m].Jc;
524 Jcsig[i][m] = kkkphi[m].Jcsig;
525 }
526 break;
527 case edPsi:
528 calc_distribution_props(nbin,histmp,-M_PI,NKKKPSI,kkkpsi,&S2);
529
530 for(m=0; (m<NKKKPSI); m++) {
531 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
532 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
533 }
534 break;
535 case edChi1:
536 calc_distribution_props(nbin,histmp,-M_PI,NKKKCHI,kkkchi1,&S2);
537 for(m=0; (m<NKKKCHI); m++) {
538 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
539 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
540 }
541 break;
542 default: /* covers edOmega and higher Chis than Chi1 */
543 calc_distribution_props(nbin,histmp,-M_PI,0,NULL,&S2);
544 break;
545 }
546 dlist[i].S2[Dih] = S2;
547
548 /* Sum distribution per amino acid type as well */
549 for(k=0; (k<nbin); k++) {
550 his_aa[Dih][dlist[i].index][k] += histmp[k];
551 histmp[k] = 0;
552 }
553 j++;
554 } else { /* dihed not defined */
555 dlist[i].S2[Dih] = 0.0 ;
556 }
557 }
558 }
559 sfree(histmp);
560
561 /* Print out Jcouplings */
562 fprintf(log,"\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
563 fprintf(log,"Residue ");
564 for(i=0; (i<NKKKPHI); i++)
565 fprintf(log,"%7s SD",kkkphi[i].name);
566 for(i=0; (i<NKKKPSI); i++)
567 fprintf(log,"%7s SD",kkkpsi[i].name);
568 for(i=0; (i<NKKKCHI); i++)
569 fprintf(log,"%7s SD",kkkchi1[i].name);
570 fprintf(log,"\n");
571 for(i=0; (i<NJC+1); i++)
572 fprintf(log,"------------");
573 fprintf(log,"\n");
574 for(i=0; (i<nlist); i++) {
575 fprintf(log,"%-10s",dlist[i].name);
576 for(j=0; (j<NJC); j++)
577 fprintf(log," %5.2f %4.2f",Jc[i][j],Jcsig[i][j]);
578 fprintf(log,"\n");
579 }
580 fprintf(log,"\n");
581
582 /* and to -jc file... */
583 if (bDo_jc) {
584 fp=xvgropen(fn,"\\S3\\NJ-Couplings from Karplus Equation","Residue",
585 "Coupling",oenv);
586 snew(leg,NJC);
587 for(i=0; (i<NKKKPHI); i++){
588 leg[i] = strdup(kkkphi[i].name);
589 }
590 for(i=0; (i<NKKKPSI); i++){
591 leg[i+NKKKPHI]=strdup(kkkpsi[i].name);
592 }
593 for(i=0; (i<NKKKCHI); i++){
594 leg[i+NKKKPHI+NKKKPSI]=strdup(kkkchi1[i].name);
595 }
596 xvgr_legend(fp,NJC,leg,oenv);
597 fprintf(fp,"%5s ","#Res.");
598 for(i=0; (i<NJC); i++)
599 fprintf(fp,"%10s ",leg[i]);
600 fprintf(fp,"\n");
601 for(i=0; (i<nlist); i++) {
602 fprintf(fp,"%5d ",dlist[i].resnr);
603 for(j=0; (j<NJC); j++)
604 fprintf(fp," %8.3f",Jc[i][j]);
605 fprintf(fp,"\n");
606 }
607 ffclose(fp);
608 for(i=0; (i<NJC); i++)
609 sfree(leg[i]);
610 }
611 /* finished -jc stuff */
612
613 snew(normhisto,nbin);
614 for(i=0; (i<naa); i++) {
615 for(Dih=0; (Dih<edMax); Dih++){
616 /* First check whether something is in there */
617 for(j=0; (j<nbin); j++)
618 if (his_aa[Dih][i][j] != 0)
619 break;
620 if ((j < nbin) &&
621 ((bPhi && (Dih==edPhi)) ||
622 (bPsi && (Dih==edPsi)) ||
623 (bOmega &&(Dih==edOmega)) ||
624 (bChi && (Dih>=edChi1)))) {
625 if (bNormalize)
626 normalize_histo(nbin,his_aa[Dih][i],(360.0/nbin),normhisto);
627
628 switch (Dih) {
629 case edPhi:
630 sprintf(hisfile,"histo-phi%s",aa[i]);
631 sprintf(title,"\\xf\\f{} Distribution for %s",aa[i]);
632 break;
633 case edPsi:
634 sprintf(hisfile,"histo-psi%s",aa[i]);
635 sprintf(title,"\\xy\\f{} Distribution for %s",aa[i]);
636 break;
637 case edOmega:
638 sprintf(hisfile,"histo-omega%s",aa[i]);
639 sprintf(title,"\\xw\\f{} Distribution for %s",aa[i]);
640 break;
641 default:
642 sprintf(hisfile,"histo-chi%d%s",Dih-NONCHI+1,aa[i]);
643 sprintf(title,"\\xc\\f{}\\s%d\\N Distribution for %s",
644 Dih-NONCHI+1,aa[i]);
645 }
646 strcpy(hhisfile,hisfile);
647 strcat(hhisfile,".xvg");
648 fp=xvgropen(hhisfile,title,"Degrees","",oenv);
649 fprintf(fp,"@ with g0\n");
650 xvgr_world(fp,-180,0,180,0.1,oenv);
651 fprintf(fp,"# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
652 fprintf(fp,"@ xaxis tick on\n");
653 fprintf(fp,"@ xaxis tick major 90\n");
654 fprintf(fp,"@ xaxis tick minor 30\n");
655 fprintf(fp,"@ xaxis ticklabel prec 0\n");
656 fprintf(fp,"@ yaxis tick off\n");
657 fprintf(fp,"@ yaxis ticklabel off\n");
658 fprintf(fp,"@ type xy\n");
659 if (bSSHisto) {
660 for(k=0; (k<3); k++) {
661 sprintf(sshisfile,"%s-%s.xvg",hisfile,sss[k]);
662 ssfp[k] = ffopen(sshisfile,"w");
663 }
664 }
665 for(j=0; (j<nbin); j++) {
666 angle = -180 + (360/nbin)*j ;
667 if (bNormalize)
668 fprintf(fp,"%5d %10g\n",angle,normhisto[j]);
669 else
670 fprintf(fp,"%5d %10d\n",angle,his_aa[Dih][i][j]);
671 if (bSSHisto)
672 for(k=0; (k<3); k++)
673 fprintf(ssfp[k],"%5d %10d\n",angle,
674 his_aa_ss[k][i][Dih][j]);
675 }
676 fprintf(fp,"&\n");
677 ffclose(fp);
678 if (bSSHisto) {
679 for(k=0; (k<3); k++) {
680 fprintf(ssfp[k],"&\n");
681 ffclose(ssfp[k]);
682 }
683 }
684 }
685 }
686 }
687 sfree(normhisto);
688
689 if (bSSHisto) {
690 /* Four dimensional array... Very cool */
691 for(i=0; (i<3); i++) {
692 for(j=0; (j<=naa); j++) {
693 for(Dih=0; (Dih<edMax); Dih++)
694 sfree(his_aa_ss[i][j][Dih]);
695 sfree(his_aa_ss[i][j]);
696 }
697 sfree(his_aa_ss[i]);
698 }
699 sfree(his_aa_ss);
700 sfree(ss_str);
701 }
702 }
703
704 static FILE *rama_file(const char *fn,const char *title,const char *xaxis,
705 const char *yaxis,const output_env_t oenv)
706 {
707 FILE *fp;
708
709 fp = xvgropen(fn,title,xaxis,yaxis,oenv);
710 fprintf(fp,"@ with g0\n");
711 xvgr_world(fp,-180,-180,180,180,oenv);
712 fprintf(fp,"@ xaxis tick on\n");
713 fprintf(fp,"@ xaxis tick major 90\n");
714 fprintf(fp,"@ xaxis tick minor 30\n");
715 fprintf(fp,"@ xaxis ticklabel prec 0\n");
716 fprintf(fp,"@ yaxis tick on\n");
717 fprintf(fp,"@ yaxis tick major 90\n");
718 fprintf(fp,"@ yaxis tick minor 30\n");
719 fprintf(fp,"@ yaxis ticklabel prec 0\n");
720 fprintf(fp,"@ s0 type xy\n");
721 fprintf(fp,"@ s0 symbol 2\n");
722 fprintf(fp,"@ s0 symbol size 0.410000\n");
723 fprintf(fp,"@ s0 symbol fill 1\n");
724 fprintf(fp,"@ s0 symbol color 1\n");
725 fprintf(fp,"@ s0 symbol linewidth 1\n");
726 fprintf(fp,"@ s0 symbol linestyle 1\n");
727 fprintf(fp,"@ s0 symbol center false\n");
728 fprintf(fp,"@ s0 symbol char 0\n");
729 fprintf(fp,"@ s0 skip 0\n");
730 fprintf(fp,"@ s0 linestyle 0\n");
731 fprintf(fp,"@ s0 linewidth 1\n");
732 fprintf(fp,"@ type xy\n");
733
734 return fp;
735 }
736
737 static void do_rama(int nf,int nlist,t_dlist dlist[],real **dih,
738 bool bViol,bool bRamOmega,const output_env_t oenv)
739 {
740 FILE *fp,*gp=NULL;
741 bool bOm;
742 char fn[256];
743 int i,j,k,Xi1,Xi2,Phi,Psi,Om=0,nlevels;
744 #define NMAT 120
745 real **mat=NULL,phi,psi,omega,axis[NMAT],lo,hi;
746 t_rgb rlo = { 1.0, 0.0, 0.0 };
747 t_rgb rmid= { 1.0, 1.0, 1.0 };
748 t_rgb rhi = { 0.0, 0.0, 1.0 };
749
750 for(i=0; (i<nlist); i++) {
751 if ((has_dihedral(edPhi,&(dlist[i]))) &&
752 (has_dihedral(edPsi,&(dlist[i])))) {
753 sprintf(fn,"ramaPhiPsi%s.xvg",dlist[i].name);
754 fp = rama_file(fn,"Ramachandran Plot",
755 "\\8f\\4 (deg)","\\8y\\4 (deg)",oenv);
756 bOm = bRamOmega && has_dihedral(edOmega,&(dlist[i]));
757 if (bOm) {
758 Om = dlist[i].j0[edOmega];
759 snew(mat,NMAT);
760 for(j=0; (j<NMAT); j++) {
761 snew(mat[j],NMAT);
762 axis[j] = -180+(360*j)/NMAT;
763 }
764 }
765 if (bViol) {
766 sprintf(fn,"violPhiPsi%s.xvg",dlist[i].name);
767 gp = ffopen(fn,"w");
768 }
769 Phi = dlist[i].j0[edPhi];
770 Psi = dlist[i].j0[edPsi];
771 for(j=0; (j<nf); j++) {
772 phi = RAD2DEG*dih[Phi][j];
773 psi = RAD2DEG*dih[Psi][j];
774 fprintf(fp,"%10g %10g\n",phi,psi);
775 if (bViol)
776 fprintf(gp,"%d\n",!bAllowed(dih[Phi][j],RAD2DEG*dih[Psi][j]));
777 if (bOm) {
778 omega = RAD2DEG*dih[Om][j];
779 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
780 += omega;
781 }
782 }
783 if (bViol)
784 ffclose(gp);
785 ffclose(fp);
786 if (bOm) {
787 sprintf(fn,"ramomega%s.xpm",dlist[i].name);
788 fp = ffopen(fn,"w");
789 lo = hi = 0;
790 for(j=0; (j<NMAT); j++)
791 for(k=0; (k<NMAT); k++) {
792 mat[j][k] /= nf;
793 lo=min(mat[j][k],lo);
794 hi=max(mat[j][k],hi);
795 }
796 /* Symmetrise */
797 if (fabs(lo) > fabs(hi))
798 hi = -lo;
799 else
800 lo = -hi;
801 /* Add 180 */
802 for(j=0; (j<NMAT); j++)
803 for(k=0; (k<NMAT); k++)
804 mat[j][k] += 180;
805 lo += 180;
806 hi += 180;
807 nlevels = 20;
808 write_xpm3(fp,0,"Omega/Ramachandran Plot","Deg","Phi","Psi",
809 NMAT,NMAT,axis,axis,mat,lo,180.0,hi,rlo,rmid,rhi,&nlevels);
810 ffclose(fp);
811 for(j=0; (j<NMAT); j++)
812 sfree(mat[j]);
813 sfree(mat);
814 }
815 }
816 if ((has_dihedral(edChi1,&(dlist[i]))) &&
817 (has_dihedral(edChi2,&(dlist[i])))) {
818 sprintf(fn,"ramaX1X2%s.xvg",dlist[i].name);
819 fp = rama_file(fn,"\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
820 "\\8c\\4\\s1\\N (deg)","\\8c\\4\\s2\\N (deg)",oenv);
821 Xi1 = dlist[i].j0[edChi1];
822 Xi2 = dlist[i].j0[edChi2];
823 for(j=0; (j<nf); j++)
824 fprintf(fp,"%10g %10g\n",RAD2DEG*dih[Xi1][j],RAD2DEG*dih[Xi2][j]);
825 ffclose(fp);
826 }
827 else
828 fprintf(stderr,"No chi1 & chi2 angle for %s\n",dlist[i].name);
829 }
830 }
831
832
833 static void print_transitions(const char *fn,int maxchi,int nlist,
834 t_dlist dlist[], t_atoms *atoms,rvec x[],
835 matrix box, bool bPhi,bool bPsi,bool bChi,real dt,
836 const output_env_t oenv)
837 {
838 /* based on order_params below */
839 FILE *fp;
840 int nh[edMax];
841 int i,Dih,Xi;
842
843 /* must correspond with enum in pp2shift.h:38 */
844 char *leg[edMax];
845 #define NLEG asize(leg)
846
847 leg[0] = strdup("Phi");
848 leg[1] = strdup("Psi");
849 leg[2] = strdup("Omega");
850 leg[3] = strdup("Chi1");
851 leg[4] = strdup("Chi2");
852 leg[5] = strdup("Chi3");
853 leg[6] = strdup("Chi4");
854 leg[7] = strdup("Chi5");
855 leg[8] = strdup("Chi6");
856
857 /* Print order parameters */
858 fp=xvgropen(fn,"Dihedral Rotamer Transitions","Residue","Transitions/ns",
859 oenv);
860 xvgr_legend(fp,NONCHI+maxchi,leg,oenv);
861
862 for (Dih=0; (Dih<edMax); Dih++)
863 nh[Dih]=0;
864
865 fprintf(fp,"%5s ","#Res.");
866 fprintf(fp,"%10s %10s %10s ",leg[edPhi],leg[edPsi],leg[edOmega]);
867 for (Xi=0; Xi<maxchi; Xi++)
868 fprintf(fp,"%10s ",leg[NONCHI+Xi]);
869 fprintf(fp,"\n");
870
871 for(i=0; (i<nlist); i++) {
872 fprintf(fp,"%5d ",dlist[i].resnr);
873 for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
874 fprintf(fp,"%10.3f ",dlist[i].ntr[Dih]/dt);
875 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
876 fprintf(fp,"\n");
877 }
878 ffclose(fp);
879 }
880
881 static void order_params(FILE *log,
882 const char *fn,int maxchi,int nlist,t_dlist dlist[],
883 const char *pdbfn,real bfac_init,
884 t_atoms *atoms,rvec x[],int ePBC,matrix box,
885 bool bPhi,bool bPsi,bool bChi,const output_env_t oenv)
886 {
887 FILE *fp;
888 int nh[edMax];
889 int i,Dih,Xi;
890 real S2Max, S2Min;
891
892 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
893 const char *const_leg[2+edMax]= { "S2Min","S2Max","Phi","Psi","Omega",
894 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
895 "Chi6" };
896 #define NLEG asize(leg)
897
898 char *leg[2+edMax];
899
900 for(i=0;i<NLEG;i++)
901 leg[i]=strdup(const_leg[i]);
902
903 /* Print order parameters */
904 fp=xvgropen(fn,"Dihedral Order Parameters","Residue","S2",oenv);
905 xvgr_legend(fp,2+NONCHI+maxchi,leg,oenv);
906
907 for (Dih=0; (Dih<edMax); Dih++)
908 nh[Dih]=0;
909
910 fprintf(fp,"%5s ","#Res.");
911 fprintf(fp,"%10s %10s ",leg[0],leg[1]);
912 fprintf(fp,"%10s %10s %10s ",leg[2+edPhi],leg[2+edPsi],leg[2+edOmega]);
913 for (Xi=0; Xi<maxchi; Xi++)
914 fprintf(fp,"%10s ",leg[2+NONCHI+Xi]);
915 fprintf(fp,"\n");
916
917 for(i=0; (i<nlist); i++) {
918 S2Max=-10;
919 S2Min=10;
920 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
921 if (dlist[i].S2[Dih]!=0) {
922 if (dlist[i].S2[Dih] > S2Max)
923 S2Max=dlist[i].S2[Dih];
924 if (dlist[i].S2[Dih] < S2Min)
925 S2Min=dlist[i].S2[Dih];
926 }
927 if (dlist[i].S2[Dih] > 0.8)
928 nh[Dih]++;
929 }
930 fprintf(fp,"%5d ",dlist[i].resnr);
931 fprintf(fp,"%10.3f %10.3f ",S2Min,S2Max);
932 for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
933 fprintf(fp,"%10.3f ",dlist[i].S2[Dih]);
934 fprintf(fp,"\n");
935 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
936 }
937 ffclose(fp);
938
939 if (pdbfn) {
940 real x0,y0,z0;
941
942 for(i=0; (i<atoms->nr); i++)
943 atoms->pdbinfo[i].bfac=bfac_init;
944
945 for(i=0; (i<nlist); i++) {
946 atoms->pdbinfo[dlist[i].atm.N].bfac=-dlist[i].S2[0];/* Phi */
947 atoms->pdbinfo[dlist[i].atm.H].bfac=-dlist[i].S2[0];/* Phi */
948 atoms->pdbinfo[dlist[i].atm.C].bfac=-dlist[i].S2[1];/* Psi */
949 atoms->pdbinfo[dlist[i].atm.O].bfac=-dlist[i].S2[1];/* Psi */
950 for (Xi=0; (Xi<maxchi); Xi++) { /* Chi's */
951 if (dlist[i].atm.Cn[Xi+3]!=-1) {
952 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac=-dlist[i].S2[NONCHI+Xi];
953 }
954 }
955 }
956
957 fp=ffopen(pdbfn,"w");
958 fprintf(fp,"REMARK generated by g_chi\n");
959 fprintf(fp,"REMARK "
960 "B-factor field contains negative of dihedral order parameters\n");
961 write_pdbfile(fp,NULL,atoms,x,ePBC,box,0,0,NULL);
962 x0=y0=z0=1000.0;
963 for (i=0; (i<atoms->nr); i++) {
964 x0 = min(x0, x[i][XX]);
965 y0 = min(y0, x[i][YY]);
966 z0 = min(z0, x[i][ZZ]);
967 }
968 x0*=10.0;/* nm -> angstrom */
969 y0*=10.0;/* nm -> angstrom */
970 z0*=10.0;/* nm -> angstrom */
971 for (i=0; (i<10); i++)
972 fprintf(fp,pdbformat,"ATOM ", atoms->nr+1+i, "CA", "LEG",' ',
973 atoms->nres+1, x0, y0, z0+(1.2*i), 0.0, -0.1*i);
974 ffclose(fp);
975 }
976
977 fprintf(log,"Dihedrals with S2 > 0.8\n");
978 fprintf(log,"Dihedral: ");
979 if (bPhi) fprintf(log," Phi ");
980 if (bPsi) fprintf(log," Psi ");
981 if (bChi)
982 for(Xi=0; (Xi<maxchi); Xi++)
983 fprintf(log," %s ", leg[2+NONCHI+Xi]);
984 fprintf(log,"\nNumber: ");
985 if (bPhi) fprintf(log,"%4d ",nh[0]);
986 if (bPsi) fprintf(log,"%4d ",nh[1]);
987 if (bChi)
988 for(Xi=0; (Xi<maxchi); Xi++)
989 fprintf(log,"%4d ",nh[NONCHI+Xi]);
990 fprintf(log,"\n");
991
992 for(i=0;i<NLEG;i++)
993 sfree(leg[i]);
994
995 }
996
997 int gmx_chi(int argc,char *argv[])
998 {
999 const char *desc[] = {
1000 "g_chi computes phi, psi, omega and chi dihedrals for all your ",
1001 "amino acid backbone and sidechains.",
1002 "It can compute dihedral angle as a function of time, and as",
1003 "histogram distributions.",
1004 "The distributions (histo-(dihedral)(RESIDUE).xvg) are cumulative over all residues of each type.[PAR]",
1005 "If option [TT]-corr[tt] is given, the program will",
1006 "calculate dihedral autocorrelation functions. The function used",
1007 "is C(t) = < cos(chi(tau)) cos(chi(tau+t)) >. The use of cosines",
1008 "rather than angles themselves, resolves the problem of periodicity.",
1009 "(Van der Spoel & Berendsen (1997), [BB]Biophys. J. 72[bb], 2032-2041).",
1010 "Separate files for each dihedral of each residue",
1011 "(corr(dihedral)(RESIDUE)(nresnr).xvg) are output, as well as a",
1012 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1013 "With option [TT]-all[tt], the angles themselves as a function of time for",
1014 "each residue are printed to separate files (dihedral)(RESIDUE)(nresnr).xvg.",
1015 "These can be in radians or degrees.[PAR]",
1016 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1017 "(a) information about the number of residues of each type.[BR]",
1018 "(b) The NMR 3J coupling constants from the Karplus equation.[BR]",
1019 "(c) a table for each residue of the number of transitions between ",
1020 "rotamers per nanosecond, and the order parameter S2 of each dihedral.[BR]",
1021 "(d) a table for each residue of the rotamer occupancy.[BR]",
1022 "All rotamers are taken as 3-fold, except for omegas and chi-dihedrals",
1023 "to planar groups (i.e. chi2 of aromatics asp and asn, chi3 of glu",
1024 "and gln, and chi4 of arg), which are 2-fold. \"rotamer 0\" means ",
1025 "that the dihedral was not in the core region of each rotamer. ",
1026 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1027
1028 "The S2 order parameters are also output to an xvg file",
1029 "(argument [TT]-o[tt] ) and optionally as a pdb file with",
1030 "the S2 values as B-factor (argument [TT]-p[tt]). ",
1031 "The total number of rotamer transitions per timestep",
1032 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1033 "(argument [TT]-rt[tt]), and the 3J couplings (argument [TT]-jc[tt]), ",
1034 "can also be written to .xvg files.[PAR]",
1035
1036 "If [TT]-chi_prod[tt] is set (and maxchi > 0), cumulative rotamers, e.g.",
1037 "1+9(chi1-1)+3(chi2-1)+(chi3-1) (if the residue has three 3-fold ",
1038 "dihedrals and maxchi >= 3)",
1039 "are calculated. As before, if any dihedral is not in the core region,",
1040 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1041 "rotamers (starting with rotamer 0) are written to the file",
1042 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1043 "is given, the rotamers as functions of time",
1044 "are written to chiproduct(RESIDUE)(nresnr).xvg ",
1045 "and their occupancies to histo-chiproduct(RESIDUE)(nresnr).xvg.[PAR]",
1046
1047 "The option [TT]-r[tt] generates a contour plot of the average omega angle",
1048 "as a function of the phi and psi angles, that is, in a Ramachandran plot",
1049 "the average omega angle is plotted using color coding.",
1050
1051 };
1052
1053 const char *bugs[] = {
1054 "Produces MANY output files (up to about 4 times the number of residues in the protein, twice that if autocorrelation functions are calculated). Typically several hundred files are output.",
1055 "Phi and psi dihedrals are calculated in a non-standard way, using H-N-CA-C for phi instead of C(-)-N-CA-C, and N-CA-C-O for psi instead of N-CA-C-N(+). This causes (usually small) discrepancies with the output of other tools like g_rama.",
1056 "-r0 option does not work properly",
1057 "Rotamers with multiplicity 2 are printed in chi.log as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1058 };
1059
1060 /* defaults */
1061 static int r0=1,ndeg=1,maxchi=2;
1062 static bool bAll=FALSE;
1063 static bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
1064 static real bfac_init=-1.0,bfac_max=0;
1065 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1066 static bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
1067 static bool bNormHisto=TRUE,bChiProduct=FALSE,bHChi=FALSE,bRAD=FALSE,bPBC=TRUE;
1068 static real core_frac=0.5 ;
1069 t_pargs pa[] = {
1070 { "-r0", FALSE, etINT, {&r0},
1071 "starting residue" },
1072 { "-phi", FALSE, etBOOL, {&bPhi},
1073 "Output for Phi dihedral angles" },
1074 { "-psi", FALSE, etBOOL, {&bPsi},
1075 "Output for Psi dihedral angles" },
1076 { "-omega",FALSE, etBOOL, {&bOmega},
1077 "Output for Omega dihedrals (peptide bonds)" },
1078 { "-rama", FALSE, etBOOL, {&bRama},
1079 "Generate Phi/Psi and Chi1/Chi2 ramachandran plots" },
1080 { "-viol", FALSE, etBOOL, {&bViol},
1081 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1082 { "-periodic", FALSE, etBOOL, {&bPBC},
1083 "Print dihedral angles modulo 360 degrees" },
1084 { "-all", FALSE, etBOOL, {&bAll},
1085 "Output separate files for every dihedral." },
1086 { "-rad", FALSE, etBOOL, {&bRAD},
1087 "in angle vs time files, use radians rather than degrees."},
1088 { "-shift", FALSE, etBOOL, {&bShift},
1089 "Compute chemical shifts from Phi/Psi angles" },
1090 { "-binwidth", FALSE, etINT, {&ndeg},
1091 "bin width for histograms (degrees)" },
1092 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1093 "only the central -core_rotamer*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1094 { "-maxchi", FALSE, etENUM, {maxchistr},
1095 "calculate first ndih Chi dihedrals" },
1096 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1097 "Normalize histograms" },
1098 { "-ramomega",FALSE,etBOOL, {&bRamOmega},
1099 "compute average omega as a function of phi/psi and plot it in an xpm plot" },
1100 { "-bfact", FALSE, etREAL, {&bfac_init},
1101 "B-factor value for pdb file for atoms with no calculated dihedral order parameter"},
1102 { "-chi_prod",FALSE,etBOOL, {&bChiProduct},
1103 "compute a single cumulative rotamer for each residue"},
1104 { "-HChi",FALSE,etBOOL, {&bHChi},
1105 "Include dihedrals to sidechain hydrogens"},
1106 { "-bmax", FALSE, etREAL, {&bfac_max},
1107 "Maximum B-factor on any of the atoms that make up a dihedral, for the dihedral angle to be considere in the statistics. Applies to database work where a number of X-Ray structures is analyzed. -bmax <= 0 means no limit." }
1108 };
1109
1110 FILE *log;
1111 int natoms,nlist,naa,idum,nbin;
1112 t_atoms atoms;
1113 rvec *x;
1114 int ePBC;
1115 matrix box;
1116 char title[256],grpname[256];
1117 t_dlist *dlist;
1118 char **aa;
1119 bool bChi,bCorr,bSSHisto;
1120 bool bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
1121 real dt=0, traj_t_ns;
1122 output_env_t oenv;
1123
1124 atom_id isize,*index;
1125 int ndih,nactdih,nf;
1126 real **dih,*trans_frac,*aver_angle,*time;
1127 int i,j,**chi_lookup,*xity;
1128
1129 t_filenm fnm[] = {
1130 { efSTX, "-s", NULL, ffREAD },
1131 { efTRX, "-f", NULL, ffREAD },
1132 { efXVG, "-o", "order", ffWRITE },
1133 { efPDB, "-p", "order", ffOPTWR },
1134 { efDAT, "-ss", "ssdump", ffOPTRD },
1135 { efXVG, "-jc", "Jcoupling", ffWRITE },
1136 { efXVG, "-corr", "dihcorr",ffOPTWR },
1137 { efLOG, "-g", "chi", ffWRITE },
1138 /* add two more arguments copying from g_angle */
1139 { efXVG, "-ot", "dihtrans", ffOPTWR },
1140 { efXVG, "-oh", "trhisto", ffOPTWR },
1141 { efXVG, "-rt", "restrans", ffOPTWR },
1142 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1143 };
1144 #define NFILE asize(fnm)
1145 int npargs;
1146 t_pargs *ppa;
1147
1148 CopyRight(stderr,argv[0]);
1149 npargs = asize(pa);
1150 ppa = add_acf_pargs(&npargs,pa);
1151 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1152 NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
1153 &oenv);
1154
1155 /* Handle result from enumerated type */
1156 sscanf(maxchistr[0],"%d",&maxchi);
1157 bChi = (maxchi > 0);
1158
1159 log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
1160
1161 if (bRamOmega) {
1162 bOmega = TRUE;
1163 bPhi = TRUE;
1164 bPsi = TRUE;
1165 }
1166
1167 /* set some options */
1168 bDo_rt=(opt2bSet("-rt",NFILE,fnm));
1169 bDo_oh=(opt2bSet("-oh",NFILE,fnm));
1170 bDo_ot=(opt2bSet("-ot",NFILE,fnm));
1171 bDo_jc=(opt2bSet("-jc",NFILE,fnm));
1172 bCorr=(opt2bSet("-corr",NFILE,fnm));
1173 if (bCorr)
1174 fprintf(stderr,"Will calculate autocorrelation\n");
1175
1176 if (core_frac > 1.0 ) {
1177 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1178 core_frac=1.0 ;
1179 }
1180 if (core_frac < 0.0 ) {
1181 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1182 core_frac=0.0 ;
1183 }
1184
1185 if (maxchi > MAXCHI) {
1186 fprintf(stderr,
1187 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1188 MAXCHI, maxchi);
1189 maxchi=MAXCHI;
1190 }
1191 bSSHisto = ftp2bSet(efDAT,NFILE,fnm);
1192 nbin = 360/ndeg ;
1193
1194 /* Find the chi angles using atoms struct and a list of amino acids */
1195 get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natoms);
1196 init_t_atoms(&atoms,natoms,TRUE);
1197 snew(x,natoms);
1198 read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,NULL,&ePBC,box);
1199 fprintf(log,"Title: %s\n",title);
1200
1201 naa=get_strings("aminoacids.dat",&aa);
1202 dlist=mk_dlist(log,&atoms,&nlist,bPhi,bPsi,bChi,bHChi,maxchi,r0,naa,aa);
1203 fprintf(stderr,"%d residues with dihedrals found\n", nlist);
1204
1205 if (nlist == 0)
1206 gmx_fatal(FARGS,"No dihedrals in your structure!\n");
1207
1208 /* Make a linear index for reading all. */
1209 index=make_chi_ind(nlist,dlist,&ndih);
1210 isize=4*ndih;
1211 fprintf(stderr,"%d dihedrals found\n", ndih);
1212
1213 snew(dih,ndih);
1214
1215 /* COMPUTE ALL DIHEDRALS! */
1216 read_ang_dih(ftp2fn(efTRX,NFILE,fnm),FALSE,TRUE,FALSE,bPBC,1,&idum,
1217 &nf,&time,isize,index,&trans_frac,&aver_angle,dih,oenv);
1218
1219 dt=(time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1220 if (bCorr)
1221 {
1222 if (nf < 2)
1223 {
1224 gmx_fatal(FARGS,"Need at least 2 frames for correlation");
1225 }
1226 }
1227
1228 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1229 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1230 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1231 nactdih = reset_em_all(nlist,dlist,nf,dih,maxchi);
1232
1233 if (bAll)
1234 dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega,bRAD,oenv);
1235
1236 /* Histogramming & J coupling constants & calc of S2 order params */
1237 histogramming(log,nbin,naa,aa,nf,maxchi,dih,nlist,dlist,index,
1238 bPhi,bPsi,bOmega,bChi,
1239 bNormHisto,bSSHisto,ftp2fn(efDAT,NFILE,fnm),bfac_max,&atoms,
1240 bDo_jc,opt2fn("-jc",NFILE,fnm),oenv);
1241
1242 /* transitions
1243 *
1244 * added multiplicity */
1245
1246 snew(xity,ndih) ;
1247 mk_multiplicity_lookup(xity, maxchi, dih, nlist, dlist,ndih);
1248
1249 strcpy(grpname, "All residues, ");
1250 if(bPhi)
1251 strcat(grpname, "Phi ");
1252 if(bPsi)
1253 strcat(grpname, "Psi ");
1254 if(bOmega)
1255 strcat(grpname, "Omega ");
1256 if(bChi){
1257 strcat(grpname, "Chi 1-") ;
1258 sprintf(grpname + strlen(grpname), "%i", maxchi);
1259 }
1260
1261
1262 low_ana_dih_trans(bDo_ot, opt2fn("-ot",NFILE,fnm),
1263 bDo_oh, opt2fn("-oh",NFILE,fnm),maxchi,
1264 dih, nlist, dlist, nf, nactdih, grpname, xity,
1265 *time, dt, FALSE, core_frac,oenv) ;
1266
1267 /* Order parameters */
1268 order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
1269 ftp2fn_null(efPDB,NFILE,fnm),bfac_init,
1270 &atoms,x,ePBC,box,bPhi,bPsi,bChi,oenv);
1271
1272 /* Print ramachandran maps! */
1273 if (bRama)
1274 do_rama(nf,nlist,dlist,dih,bViol,bRamOmega,oenv);
1275
1276 if (bShift)
1277 do_pp2shifts(log,nf,nlist,dlist,dih);
1278
1279 /* rprint S^2, transitions, and rotamer occupancies to log */
1280 traj_t_ns = 0.001 * (time[nf-1]-time[0]) ;
1281 pr_dlist(log,nlist,dlist,traj_t_ns,edPrintST,bPhi,bPsi,bChi,bOmega,maxchi);
1282 pr_dlist(log,nlist,dlist,traj_t_ns,edPrintRO,bPhi,bPsi,bChi,bOmega,maxchi);
1283 ffclose(log);
1284 /* transitions to xvg */
1285 if (bDo_rt)
1286 print_transitions(opt2fn("-rt",NFILE,fnm),maxchi,nlist,dlist,
1287 &atoms,x,box,bPhi,bPsi,bChi,traj_t_ns,oenv);
1288
1289 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1290 if (bChiProduct && bChi ) {
1291 snew(chi_lookup,nlist) ;
1292 for (i=0;i<nlist;i++)
1293 snew(chi_lookup[i], maxchi) ;
1294 mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1295
1296 get_chi_product_traj(dih,nf,nactdih,nlist,
1297 maxchi,dlist,time,chi_lookup,xity,
1298 FALSE,bNormHisto, core_frac,bAll,
1299 opt2fn("-cp",NFILE,fnm),oenv);
1300
1301 for (i=0;i<nlist;i++)
1302 sfree(chi_lookup[i]);
1303 }
1304
1305 /* Correlation comes last because it fucks up the angles */
1306 if (bCorr)
1307 do_dihcorr(opt2fn("-corr",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
1308 maxchi,bPhi,bPsi,bChi,bOmega,oenv);
1309
1310
1311 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1312 do_view(oenv,opt2fn("-jc",NFILE,fnm),"-nxy");
1313 if (bCorr)
1314 do_view(oenv,opt2fn("-corr",NFILE,fnm),"-nxy");
1315
1316 thanx(stderr);
1317
1318 return 0;
1319 }
Atempt to write QMMM interface for Gamess-US