?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fixed make_edi.c
[gromacs/rigid-bodies.git] / src / tools / gmx_chi.c
blob472735256757d9a166176e9ccd4a669a81f0ad9a
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 char buf[STRLEN];
890 int i,Dih,Xi;
891 real S2Max, S2Min;
892
893 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
894 const char *const_leg[2+edMax]= { "S2Min","S2Max","Phi","Psi","Omega",
895 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
896 "Chi6" };
897 #define NLEG asize(leg)
898
899 char *leg[2+edMax];
900
901 for(i=0;i<NLEG;i++)
902 leg[i]=strdup(const_leg[i]);
903
904 /* Print order parameters */
905 fp=xvgropen(fn,"Dihedral Order Parameters","Residue","S2",oenv);
906 xvgr_legend(fp,2+NONCHI+maxchi,leg,oenv);
907
908 for (Dih=0; (Dih<edMax); Dih++)
909 nh[Dih]=0;
910
911 fprintf(fp,"%5s ","#Res.");
912 fprintf(fp,"%10s %10s ",leg[0],leg[1]);
913 fprintf(fp,"%10s %10s %10s ",leg[2+edPhi],leg[2+edPsi],leg[2+edOmega]);
914 for (Xi=0; Xi<maxchi; Xi++)
915 fprintf(fp,"%10s ",leg[2+NONCHI+Xi]);
916 fprintf(fp,"\n");
917
918 for(i=0; (i<nlist); i++) {
919 S2Max=-10;
920 S2Min=10;
921 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
922 if (dlist[i].S2[Dih]!=0) {
923 if (dlist[i].S2[Dih] > S2Max)
924 S2Max=dlist[i].S2[Dih];
925 if (dlist[i].S2[Dih] < S2Min)
926 S2Min=dlist[i].S2[Dih];
927 }
928 if (dlist[i].S2[Dih] > 0.8)
929 nh[Dih]++;
930 }
931 fprintf(fp,"%5d ",dlist[i].resnr);
932 fprintf(fp,"%10.3f %10.3f ",S2Min,S2Max);
933 for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
934 fprintf(fp,"%10.3f ",dlist[i].S2[Dih]);
935 fprintf(fp,"\n");
936 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
937 }
938 ffclose(fp);
939
940 if (NULL != pdbfn) {
941 real x0,y0,z0;
942
943 if (NULL == atoms->pdbinfo)
944 snew(atoms->pdbinfo,atoms->nr);
945 for(i=0; (i<atoms->nr); i++)
946 atoms->pdbinfo[i].bfac=bfac_init;
947
948 for(i=0; (i<nlist); i++) {
949 atoms->pdbinfo[dlist[i].atm.N].bfac=-dlist[i].S2[0];/* Phi */
950 atoms->pdbinfo[dlist[i].atm.H].bfac=-dlist[i].S2[0];/* Phi */
951 atoms->pdbinfo[dlist[i].atm.C].bfac=-dlist[i].S2[1];/* Psi */
952 atoms->pdbinfo[dlist[i].atm.O].bfac=-dlist[i].S2[1];/* Psi */
953 for (Xi=0; (Xi<maxchi); Xi++) { /* Chi's */
954 if (dlist[i].atm.Cn[Xi+3]!=-1) {
955 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac=-dlist[i].S2[NONCHI+Xi];
956 }
957 }
958 }
959
960 fp=ffopen(pdbfn,"w");
961 fprintf(fp,"REMARK generated by g_chi\n");
962 fprintf(fp,"REMARK "
963 "B-factor field contains negative of dihedral order parameters\n");
964 write_pdbfile(fp,NULL,atoms,x,ePBC,box,0,0,NULL);
965 x0=y0=z0=1000.0;
966 for (i=0; (i<atoms->nr); i++) {
967 x0 = min(x0, x[i][XX]);
968 y0 = min(y0, x[i][YY]);
969 z0 = min(z0, x[i][ZZ]);
970 }
971 x0*=10.0;/* nm -> angstrom */
972 y0*=10.0;/* nm -> angstrom */
973 z0*=10.0;/* nm -> angstrom */
974 sprintf(buf,"%s%%6.f%%6.2f\n",pdbformat);
975 for (i=0; (i<10); i++) {
976 fprintf(fp,buf,"ATOM ", atoms->nr+1+i, "CA", "LEG",' ',
977 atoms->nres+1, ' ',x0, y0, z0+(1.2*i), 0.0, -0.1*i);
978 }
979 ffclose(fp);
980 }
981
982 fprintf(log,"Dihedrals with S2 > 0.8\n");
983 fprintf(log,"Dihedral: ");
984 if (bPhi) fprintf(log," Phi ");
985 if (bPsi) fprintf(log," Psi ");
986 if (bChi)
987 for(Xi=0; (Xi<maxchi); Xi++)
988 fprintf(log," %s ", leg[2+NONCHI+Xi]);
989 fprintf(log,"\nNumber: ");
990 if (bPhi) fprintf(log,"%4d ",nh[0]);
991 if (bPsi) fprintf(log,"%4d ",nh[1]);
992 if (bChi)
993 for(Xi=0; (Xi<maxchi); Xi++)
994 fprintf(log,"%4d ",nh[NONCHI+Xi]);
995 fprintf(log,"\n");
996
997 for(i=0;i<NLEG;i++)
998 sfree(leg[i]);
999
1000 }
1001
1002 int gmx_chi(int argc,char *argv[])
1003 {
1004 const char *desc[] = {
1005 "g_chi computes phi, psi, omega and chi dihedrals for all your ",
1006 "amino acid backbone and sidechains.",
1007 "It can compute dihedral angle as a function of time, and as",
1008 "histogram distributions.",
1009 "The distributions (histo-(dihedral)(RESIDUE).xvg) are cumulative over all residues of each type.[PAR]",
1010 "If option [TT]-corr[tt] is given, the program will",
1011 "calculate dihedral autocorrelation functions. The function used",
1012 "is C(t) = < cos(chi(tau)) cos(chi(tau+t)) >. The use of cosines",
1013 "rather than angles themselves, resolves the problem of periodicity.",
1014 "(Van der Spoel & Berendsen (1997), [BB]Biophys. J. 72[bb], 2032-2041).",
1015 "Separate files for each dihedral of each residue",
1016 "(corr(dihedral)(RESIDUE)(nresnr).xvg) are output, as well as a",
1017 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1018 "With option [TT]-all[tt], the angles themselves as a function of time for",
1019 "each residue are printed to separate files (dihedral)(RESIDUE)(nresnr).xvg.",
1020 "These can be in radians or degrees.[PAR]",
1021 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1022 "(a) information about the number of residues of each type.[BR]",
1023 "(b) The NMR 3J coupling constants from the Karplus equation.[BR]",
1024 "(c) a table for each residue of the number of transitions between ",
1025 "rotamers per nanosecond, and the order parameter S2 of each dihedral.[BR]",
1026 "(d) a table for each residue of the rotamer occupancy.[BR]",
1027 "All rotamers are taken as 3-fold, except for omegas and chi-dihedrals",
1028 "to planar groups (i.e. chi2 of aromatics asp and asn, chi3 of glu",
1029 "and gln, and chi4 of arg), which are 2-fold. \"rotamer 0\" means ",
1030 "that the dihedral was not in the core region of each rotamer. ",
1031 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1032
1033 "The S2 order parameters are also output to an xvg file",
1034 "(argument [TT]-o[tt] ) and optionally as a pdb file with",
1035 "the S2 values as B-factor (argument [TT]-p[tt]). ",
1036 "The total number of rotamer transitions per timestep",
1037 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1038 "(argument [TT]-rt[tt]), and the 3J couplings (argument [TT]-jc[tt]), ",
1039 "can also be written to .xvg files.[PAR]",
1040
1041 "If [TT]-chi_prod[tt] is set (and maxchi > 0), cumulative rotamers, e.g.",
1042 "1+9(chi1-1)+3(chi2-1)+(chi3-1) (if the residue has three 3-fold ",
1043 "dihedrals and maxchi >= 3)",
1044 "are calculated. As before, if any dihedral is not in the core region,",
1045 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1046 "rotamers (starting with rotamer 0) are written to the file",
1047 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1048 "is given, the rotamers as functions of time",
1049 "are written to chiproduct(RESIDUE)(nresnr).xvg ",
1050 "and their occupancies to histo-chiproduct(RESIDUE)(nresnr).xvg.[PAR]",
1051
1052 "The option [TT]-r[tt] generates a contour plot of the average omega angle",
1053 "as a function of the phi and psi angles, that is, in a Ramachandran plot",
1054 "the average omega angle is plotted using color coding.",
1055
1056 };
1057
1058 const char *bugs[] = {
1059 "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.",
1060 "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.",
1061 "-r0 option does not work properly",
1062 "Rotamers with multiplicity 2 are printed in chi.log as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1063 };
1064
1065 /* defaults */
1066 static int r0=1,ndeg=1,maxchi=2;
1067 static bool bAll=FALSE;
1068 static bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
1069 static real bfac_init=-1.0,bfac_max=0;
1070 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1071 static bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
1072 static bool bNormHisto=TRUE,bChiProduct=FALSE,bHChi=FALSE,bRAD=FALSE,bPBC=TRUE;
1073 static real core_frac=0.5 ;
1074 t_pargs pa[] = {
1075 { "-r0", FALSE, etINT, {&r0},
1076 "starting residue" },
1077 { "-phi", FALSE, etBOOL, {&bPhi},
1078 "Output for Phi dihedral angles" },
1079 { "-psi", FALSE, etBOOL, {&bPsi},
1080 "Output for Psi dihedral angles" },
1081 { "-omega",FALSE, etBOOL, {&bOmega},
1082 "Output for Omega dihedrals (peptide bonds)" },
1083 { "-rama", FALSE, etBOOL, {&bRama},
1084 "Generate Phi/Psi and Chi1/Chi2 ramachandran plots" },
1085 { "-viol", FALSE, etBOOL, {&bViol},
1086 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1087 { "-periodic", FALSE, etBOOL, {&bPBC},
1088 "Print dihedral angles modulo 360 degrees" },
1089 { "-all", FALSE, etBOOL, {&bAll},
1090 "Output separate files for every dihedral." },
1091 { "-rad", FALSE, etBOOL, {&bRAD},
1092 "in angle vs time files, use radians rather than degrees."},
1093 { "-shift", FALSE, etBOOL, {&bShift},
1094 "Compute chemical shifts from Phi/Psi angles" },
1095 { "-binwidth", FALSE, etINT, {&ndeg},
1096 "bin width for histograms (degrees)" },
1097 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1098 "only the central -core_rotamer*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1099 { "-maxchi", FALSE, etENUM, {maxchistr},
1100 "calculate first ndih Chi dihedrals" },
1101 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1102 "Normalize histograms" },
1103 { "-ramomega",FALSE,etBOOL, {&bRamOmega},
1104 "compute average omega as a function of phi/psi and plot it in an xpm plot" },
1105 { "-bfact", FALSE, etREAL, {&bfac_init},
1106 "B-factor value for pdb file for atoms with no calculated dihedral order parameter"},
1107 { "-chi_prod",FALSE,etBOOL, {&bChiProduct},
1108 "compute a single cumulative rotamer for each residue"},
1109 { "-HChi",FALSE,etBOOL, {&bHChi},
1110 "Include dihedrals to sidechain hydrogens"},
1111 { "-bmax", FALSE, etREAL, {&bfac_max},
1112 "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." }
1113 };
1114
1115 FILE *log;
1116 int natoms,nlist,naa,idum,nbin;
1117 t_atoms atoms;
1118 rvec *x;
1119 int ePBC;
1120 matrix box;
1121 char title[256],grpname[256];
1122 t_dlist *dlist;
1123 char **aa;
1124 bool bChi,bCorr,bSSHisto;
1125 bool bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
1126 real dt=0, traj_t_ns;
1127 output_env_t oenv;
1128
1129 atom_id isize,*index;
1130 int ndih,nactdih,nf;
1131 real **dih,*trans_frac,*aver_angle,*time;
1132 int i,j,**chi_lookup,*xity;
1133
1134 t_filenm fnm[] = {
1135 { efSTX, "-s", NULL, ffREAD },
1136 { efTRX, "-f", NULL, ffREAD },
1137 { efXVG, "-o", "order", ffWRITE },
1138 { efPDB, "-p", "order", ffOPTWR },
1139 { efDAT, "-ss", "ssdump", ffOPTRD },
1140 { efXVG, "-jc", "Jcoupling", ffWRITE },
1141 { efXVG, "-corr", "dihcorr",ffOPTWR },
1142 { efLOG, "-g", "chi", ffWRITE },
1143 /* add two more arguments copying from g_angle */
1144 { efXVG, "-ot", "dihtrans", ffOPTWR },
1145 { efXVG, "-oh", "trhisto", ffOPTWR },
1146 { efXVG, "-rt", "restrans", ffOPTWR },
1147 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1148 };
1149 #define NFILE asize(fnm)
1150 int npargs;
1151 t_pargs *ppa;
1152
1153 CopyRight(stderr,argv[0]);
1154 npargs = asize(pa);
1155 ppa = add_acf_pargs(&npargs,pa);
1156 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1157 NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
1158 &oenv);
1159
1160 /* Handle result from enumerated type */
1161 sscanf(maxchistr[0],"%d",&maxchi);
1162 bChi = (maxchi > 0);
1163
1164 log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
1165
1166 if (bRamOmega) {
1167 bOmega = TRUE;
1168 bPhi = TRUE;
1169 bPsi = TRUE;
1170 }
1171
1172 /* set some options */
1173 bDo_rt=(opt2bSet("-rt",NFILE,fnm));
1174 bDo_oh=(opt2bSet("-oh",NFILE,fnm));
1175 bDo_ot=(opt2bSet("-ot",NFILE,fnm));
1176 bDo_jc=(opt2bSet("-jc",NFILE,fnm));
1177 bCorr=(opt2bSet("-corr",NFILE,fnm));
1178 if (bCorr)
1179 fprintf(stderr,"Will calculate autocorrelation\n");
1180
1181 if (core_frac > 1.0 ) {
1182 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1183 core_frac=1.0 ;
1184 }
1185 if (core_frac < 0.0 ) {
1186 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1187 core_frac=0.0 ;
1188 }
1189
1190 if (maxchi > MAXCHI) {
1191 fprintf(stderr,
1192 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1193 MAXCHI, maxchi);
1194 maxchi=MAXCHI;
1195 }
1196 bSSHisto = ftp2bSet(efDAT,NFILE,fnm);
1197 nbin = 360/ndeg ;
1198
1199 /* Find the chi angles using atoms struct and a list of amino acids */
1200 get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natoms);
1201 init_t_atoms(&atoms,natoms,TRUE);
1202 snew(x,natoms);
1203 read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,NULL,&ePBC,box);
1204 fprintf(log,"Title: %s\n",title);
1205
1206 naa=get_strings("aminoacids.dat",&aa);
1207 dlist=mk_dlist(log,&atoms,&nlist,bPhi,bPsi,bChi,bHChi,maxchi,r0,naa,aa);
1208 fprintf(stderr,"%d residues with dihedrals found\n", nlist);
1209
1210 if (nlist == 0)
1211 gmx_fatal(FARGS,"No dihedrals in your structure!\n");
1212
1213 /* Make a linear index for reading all. */
1214 index=make_chi_ind(nlist,dlist,&ndih);
1215 isize=4*ndih;
1216 fprintf(stderr,"%d dihedrals found\n", ndih);
1217
1218 snew(dih,ndih);
1219
1220 /* COMPUTE ALL DIHEDRALS! */
1221 read_ang_dih(ftp2fn(efTRX,NFILE,fnm),FALSE,TRUE,FALSE,bPBC,1,&idum,
1222 &nf,&time,isize,index,&trans_frac,&aver_angle,dih,oenv);
1223
1224 dt=(time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1225 if (bCorr)
1226 {
1227 if (nf < 2)
1228 {
1229 gmx_fatal(FARGS,"Need at least 2 frames for correlation");
1230 }
1231 }
1232
1233 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1234 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1235 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1236 nactdih = reset_em_all(nlist,dlist,nf,dih,maxchi);
1237
1238 if (bAll)
1239 dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega,bRAD,oenv);
1240
1241 /* Histogramming & J coupling constants & calc of S2 order params */
1242 histogramming(log,nbin,naa,aa,nf,maxchi,dih,nlist,dlist,index,
1243 bPhi,bPsi,bOmega,bChi,
1244 bNormHisto,bSSHisto,ftp2fn(efDAT,NFILE,fnm),bfac_max,&atoms,
1245 bDo_jc,opt2fn("-jc",NFILE,fnm),oenv);
1246
1247 /* transitions
1248 *
1249 * added multiplicity */
1250
1251 snew(xity,ndih) ;
1252 mk_multiplicity_lookup(xity, maxchi, dih, nlist, dlist,ndih);
1253
1254 strcpy(grpname, "All residues, ");
1255 if(bPhi)
1256 strcat(grpname, "Phi ");
1257 if(bPsi)
1258 strcat(grpname, "Psi ");
1259 if(bOmega)
1260 strcat(grpname, "Omega ");
1261 if(bChi){
1262 strcat(grpname, "Chi 1-") ;
1263 sprintf(grpname + strlen(grpname), "%i", maxchi);
1264 }
1265
1266
1267 low_ana_dih_trans(bDo_ot, opt2fn("-ot",NFILE,fnm),
1268 bDo_oh, opt2fn("-oh",NFILE,fnm),maxchi,
1269 dih, nlist, dlist, nf, nactdih, grpname, xity,
1270 *time, dt, FALSE, core_frac,oenv) ;
1271
1272 /* Order parameters */
1273 order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
1274 ftp2fn_null(efPDB,NFILE,fnm),bfac_init,
1275 &atoms,x,ePBC,box,bPhi,bPsi,bChi,oenv);
1276
1277 /* Print ramachandran maps! */
1278 if (bRama)
1279 do_rama(nf,nlist,dlist,dih,bViol,bRamOmega,oenv);
1280
1281 if (bShift)
1282 do_pp2shifts(log,nf,nlist,dlist,dih);
1283
1284 /* rprint S^2, transitions, and rotamer occupancies to log */
1285 traj_t_ns = 0.001 * (time[nf-1]-time[0]) ;
1286 pr_dlist(log,nlist,dlist,traj_t_ns,edPrintST,bPhi,bPsi,bChi,bOmega,maxchi);
1287 pr_dlist(log,nlist,dlist,traj_t_ns,edPrintRO,bPhi,bPsi,bChi,bOmega,maxchi);
1288 ffclose(log);
1289 /* transitions to xvg */
1290 if (bDo_rt)
1291 print_transitions(opt2fn("-rt",NFILE,fnm),maxchi,nlist,dlist,
1292 &atoms,x,box,bPhi,bPsi,bChi,traj_t_ns,oenv);
1293
1294 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1295 if (bChiProduct && bChi ) {
1296 snew(chi_lookup,nlist) ;
1297 for (i=0;i<nlist;i++)
1298 snew(chi_lookup[i], maxchi) ;
1299 mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1300
1301 get_chi_product_traj(dih,nf,nactdih,nlist,
1302 maxchi,dlist,time,chi_lookup,xity,
1303 FALSE,bNormHisto, core_frac,bAll,
1304 opt2fn("-cp",NFILE,fnm),oenv);
1305
1306 for (i=0;i<nlist;i++)
1307 sfree(chi_lookup[i]);
1308 }
1309
1310 /* Correlation comes last because it fucks up the angles */
1311 if (bCorr)
1312 do_dihcorr(opt2fn("-corr",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
1313 maxchi,bPhi,bPsi,bChi,bOmega,oenv);
1314
1315
1316 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1317 do_view(oenv,opt2fn("-jc",NFILE,fnm),"-nxy");
1318 if (bCorr)
1319 do_view(oenv,opt2fn("-corr",NFILE,fnm),"-nxy");
1320
1321 thanx(stderr);
1322
1323 return 0;
1324 }