?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / tools / gmx_chi.c
blobbf18bdba3f31fba695b2ec3e0e06394c9042175f
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 gmx_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 (gmx_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 gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,gmx_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, gmx_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 gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,gmx_bool bOmega, gmx_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,gmx_residuetype_t rt,
389 int nf,int maxchi,real **dih,
390 int nlist,t_dlist dlist[],
391 atom_id index[],
392 gmx_bool bPhi,gmx_bool bPsi,gmx_bool bOmega,gmx_bool bChi,
393 gmx_bool bNormalize,gmx_bool bSSHisto,const char *ssdump,
394 real bfac_max,t_atoms *atoms,
395 gmx_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 gmx_bool bBfac,bOccup;
427 char hisfile[256],hhisfile[256],sshisfile[256],title[256],*ss_str=NULL;
428 char **leg;
429 const char *residue_name;
430 int rt_size;
431
432 rt_size = gmx_residuetype_get_size(rt);
433 if (bSSHisto) {
434 fp = ffopen(ssdump,"r");
435 if(1 != fscanf(fp,"%d",&nres))
436 {
437 gmx_fatal(FARGS,"Error reading from file %s",ssdump);
438 }
439
440 snew(ss_str,nres+1);
441 if( 1 != fscanf(fp,"%s",ss_str))
442 {
443 gmx_fatal(FARGS,"Error reading from file %s",ssdump);
444 }
445
446 ffclose(fp);
447 /* Four dimensional array... Very cool */
448 snew(his_aa_ss,3);
449 for(i=0; (i<3); i++) {
450 snew(his_aa_ss[i],rt_size+1);
451 for(j=0; (j<=rt_size); j++) {
452 snew(his_aa_ss[i][j],edMax);
453 for(Dih=0; (Dih<edMax); Dih++)
454 snew(his_aa_ss[i][j][Dih],nbin+1);
455 }
456 }
457 }
458 snew(his_aa,edMax);
459 for(Dih=0; (Dih<edMax); Dih++) {
460 snew(his_aa[Dih],rt_size+1);
461 for(i=0; (i<=rt_size); i++) {
462 snew(his_aa[Dih][i],nbin+1);
463 }
464 }
465 snew(histmp,nbin);
466
467 snew(Jc,nlist);
468 snew(Jcsig,nlist);
469 for(i=0; (i<nlist); i++) {
470 snew(Jc[i],NJC);
471 snew(Jcsig[i],NJC);
472 }
473
474 j=0;
475 n=0;
476 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
477 for(i=0; (i<nlist); i++) {
478 if (((Dih < edOmega) ) ||
479 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
480 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
481 make_histo(log,nf,dih[j],nbin,histmp,-M_PI,M_PI);
482
483 if (bSSHisto) {
484 /* Assume there is only one structure, the first.
485 * Compute index in histogram.
486 */
487 /* Check the atoms to see whether their B-factors are low enough
488 * Check atoms to see their occupancy is 1.
489 */
490 bBfac = bOccup = TRUE;
491 for(nn=0; (nn<4); nn++,n++) {
492 bBfac = bBfac && (atoms->pdbinfo[index[n]].bfac <= bfac_max);
493 bOccup = bOccup && (atoms->pdbinfo[index[n]].occup == 1);
494 }
495 if (bOccup && ((bfac_max <= 0) || ((bfac_max > 0) && bBfac))) {
496 hindex = ((dih[j][0]+M_PI)*nbin)/(2*M_PI);
497 range_check(hindex,0,nbin);
498
499 /* Assign dihedral to either of the structure determined
500 * histograms
501 */
502 switch(ss_str[dlist[i].resnr]) {
503 case 'E':
504 his_aa_ss[0][dlist[i].index][Dih][hindex]++;
505 break;
506 case 'H':
507 his_aa_ss[1][dlist[i].index][Dih][hindex]++;
508 break;
509 default:
510 his_aa_ss[2][dlist[i].index][Dih][hindex]++;
511 break;
512 }
513 }
514 else if (debug)
515 fprintf(debug,"Res. %d has imcomplete occupancy or bfacs > %g\n",
516 dlist[i].resnr,bfac_max);
517 }
518 else
519 n += 4;
520
521 switch (Dih) {
522 case edPhi:
523 calc_distribution_props(nbin,histmp,-M_PI,NKKKPHI,kkkphi,&S2);
524
525 for(m=0; (m<NKKKPHI); m++) {
526 Jc[i][m] = kkkphi[m].Jc;
527 Jcsig[i][m] = kkkphi[m].Jcsig;
528 }
529 break;
530 case edPsi:
531 calc_distribution_props(nbin,histmp,-M_PI,NKKKPSI,kkkpsi,&S2);
532
533 for(m=0; (m<NKKKPSI); m++) {
534 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
535 Jcsig[i][NKKKPHI+m] = kkkpsi[m].Jcsig;
536 }
537 break;
538 case edChi1:
539 calc_distribution_props(nbin,histmp,-M_PI,NKKKCHI,kkkchi1,&S2);
540 for(m=0; (m<NKKKCHI); m++) {
541 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
542 Jcsig[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jcsig;
543 }
544 break;
545 default: /* covers edOmega and higher Chis than Chi1 */
546 calc_distribution_props(nbin,histmp,-M_PI,0,NULL,&S2);
547 break;
548 }
549 dlist[i].S2[Dih] = S2;
550
551 /* Sum distribution per amino acid type as well */
552 for(k=0; (k<nbin); k++) {
553 his_aa[Dih][dlist[i].index][k] += histmp[k];
554 histmp[k] = 0;
555 }
556 j++;
557 } else { /* dihed not defined */
558 dlist[i].S2[Dih] = 0.0 ;
559 }
560 }
561 }
562 sfree(histmp);
563
564 /* Print out Jcouplings */
565 fprintf(log,"\n *** J-Couplings from simulation (plus std. dev.) ***\n\n");
566 fprintf(log,"Residue ");
567 for(i=0; (i<NKKKPHI); i++)
568 fprintf(log,"%7s SD",kkkphi[i].name);
569 for(i=0; (i<NKKKPSI); i++)
570 fprintf(log,"%7s SD",kkkpsi[i].name);
571 for(i=0; (i<NKKKCHI); i++)
572 fprintf(log,"%7s SD",kkkchi1[i].name);
573 fprintf(log,"\n");
574 for(i=0; (i<NJC+1); i++)
575 fprintf(log,"------------");
576 fprintf(log,"\n");
577 for(i=0; (i<nlist); i++) {
578 fprintf(log,"%-10s",dlist[i].name);
579 for(j=0; (j<NJC); j++)
580 fprintf(log," %5.2f %4.2f",Jc[i][j],Jcsig[i][j]);
581 fprintf(log,"\n");
582 }
583 fprintf(log,"\n");
584
585 /* and to -jc file... */
586 if (bDo_jc) {
587 fp=xvgropen(fn,"\\S3\\NJ-Couplings from Karplus Equation","Residue",
588 "Coupling",oenv);
589 snew(leg,NJC);
590 for(i=0; (i<NKKKPHI); i++){
591 leg[i] = strdup(kkkphi[i].name);
592 }
593 for(i=0; (i<NKKKPSI); i++){
594 leg[i+NKKKPHI]=strdup(kkkpsi[i].name);
595 }
596 for(i=0; (i<NKKKCHI); i++){
597 leg[i+NKKKPHI+NKKKPSI]=strdup(kkkchi1[i].name);
598 }
599 xvgr_legend(fp,NJC,(const char**)leg,oenv);
600 fprintf(fp,"%5s ","#Res.");
601 for(i=0; (i<NJC); i++)
602 fprintf(fp,"%10s ",leg[i]);
603 fprintf(fp,"\n");
604 for(i=0; (i<nlist); i++) {
605 fprintf(fp,"%5d ",dlist[i].resnr);
606 for(j=0; (j<NJC); j++)
607 fprintf(fp," %8.3f",Jc[i][j]);
608 fprintf(fp,"\n");
609 }
610 ffclose(fp);
611 for(i=0; (i<NJC); i++)
612 sfree(leg[i]);
613 }
614 /* finished -jc stuff */
615
616 snew(normhisto,nbin);
617 for(i=0; (i<rt_size); i++) {
618 for(Dih=0; (Dih<edMax); Dih++){
619 /* First check whether something is in there */
620 for(j=0; (j<nbin); j++)
621 if (his_aa[Dih][i][j] != 0)
622 break;
623 if ((j < nbin) &&
624 ((bPhi && (Dih==edPhi)) ||
625 (bPsi && (Dih==edPsi)) ||
626 (bOmega &&(Dih==edOmega)) ||
627 (bChi && (Dih>=edChi1)))) {
628 if (bNormalize)
629 normalize_histo(nbin,his_aa[Dih][i],(360.0/nbin),normhisto);
630
631 residue_name = gmx_residuetype_get_name(rt,i);
632 switch (Dih) {
633 case edPhi:
634 sprintf(hisfile,"histo-phi%s",residue_name);
635 sprintf(title,"\\xf\\f{} Distribution for %s",residue_name);
636 break;
637 case edPsi:
638 sprintf(hisfile,"histo-psi%s",residue_name);
639 sprintf(title,"\\xy\\f{} Distribution for %s",residue_name);
640 break;
641 case edOmega:
642 sprintf(hisfile,"histo-omega%s",residue_name);
643 sprintf(title,"\\xw\\f{} Distribution for %s",residue_name);
644 break;
645 default:
646 sprintf(hisfile,"histo-chi%d%s",Dih-NONCHI+1,residue_name);
647 sprintf(title,"\\xc\\f{}\\s%d\\N Distribution for %s",
648 Dih-NONCHI+1,residue_name);
649 }
650 strcpy(hhisfile,hisfile);
651 strcat(hhisfile,".xvg");
652 fp=xvgropen(hhisfile,title,"Degrees","",oenv);
653 fprintf(fp,"@ with g0\n");
654 xvgr_world(fp,-180,0,180,0.1,oenv);
655 fprintf(fp,"# this effort to set graph size fails unless you run with -autoscale none or -autoscale y flags\n");
656 fprintf(fp,"@ xaxis tick on\n");
657 fprintf(fp,"@ xaxis tick major 90\n");
658 fprintf(fp,"@ xaxis tick minor 30\n");
659 fprintf(fp,"@ xaxis ticklabel prec 0\n");
660 fprintf(fp,"@ yaxis tick off\n");
661 fprintf(fp,"@ yaxis ticklabel off\n");
662 fprintf(fp,"@ type xy\n");
663 if (bSSHisto) {
664 for(k=0; (k<3); k++) {
665 sprintf(sshisfile,"%s-%s.xvg",hisfile,sss[k]);
666 ssfp[k] = ffopen(sshisfile,"w");
667 }
668 }
669 for(j=0; (j<nbin); j++) {
670 angle = -180 + (360/nbin)*j ;
671 if (bNormalize)
672 fprintf(fp,"%5d %10g\n",angle,normhisto[j]);
673 else
674 fprintf(fp,"%5d %10d\n",angle,his_aa[Dih][i][j]);
675 if (bSSHisto)
676 for(k=0; (k<3); k++)
677 fprintf(ssfp[k],"%5d %10d\n",angle,
678 his_aa_ss[k][i][Dih][j]);
679 }
680 fprintf(fp,"&\n");
681 ffclose(fp);
682 if (bSSHisto) {
683 for(k=0; (k<3); k++) {
684 fprintf(ssfp[k],"&\n");
685 ffclose(ssfp[k]);
686 }
687 }
688 }
689 }
690 }
691 sfree(normhisto);
692
693 if (bSSHisto) {
694 /* Four dimensional array... Very cool */
695 for(i=0; (i<3); i++) {
696 for(j=0; (j<=rt_size); j++) {
697 for(Dih=0; (Dih<edMax); Dih++)
698 sfree(his_aa_ss[i][j][Dih]);
699 sfree(his_aa_ss[i][j]);
700 }
701 sfree(his_aa_ss[i]);
702 }
703 sfree(his_aa_ss);
704 sfree(ss_str);
705 }
706 }
707
708 static FILE *rama_file(const char *fn,const char *title,const char *xaxis,
709 const char *yaxis,const output_env_t oenv)
710 {
711 FILE *fp;
712
713 fp = xvgropen(fn,title,xaxis,yaxis,oenv);
714 fprintf(fp,"@ with g0\n");
715 xvgr_world(fp,-180,-180,180,180,oenv);
716 fprintf(fp,"@ xaxis tick on\n");
717 fprintf(fp,"@ xaxis tick major 90\n");
718 fprintf(fp,"@ xaxis tick minor 30\n");
719 fprintf(fp,"@ xaxis ticklabel prec 0\n");
720 fprintf(fp,"@ yaxis tick on\n");
721 fprintf(fp,"@ yaxis tick major 90\n");
722 fprintf(fp,"@ yaxis tick minor 30\n");
723 fprintf(fp,"@ yaxis ticklabel prec 0\n");
724 fprintf(fp,"@ s0 type xy\n");
725 fprintf(fp,"@ s0 symbol 2\n");
726 fprintf(fp,"@ s0 symbol size 0.410000\n");
727 fprintf(fp,"@ s0 symbol fill 1\n");
728 fprintf(fp,"@ s0 symbol color 1\n");
729 fprintf(fp,"@ s0 symbol linewidth 1\n");
730 fprintf(fp,"@ s0 symbol linestyle 1\n");
731 fprintf(fp,"@ s0 symbol center false\n");
732 fprintf(fp,"@ s0 symbol char 0\n");
733 fprintf(fp,"@ s0 skip 0\n");
734 fprintf(fp,"@ s0 linestyle 0\n");
735 fprintf(fp,"@ s0 linewidth 1\n");
736 fprintf(fp,"@ type xy\n");
737
738 return fp;
739 }
740
741 static void do_rama(int nf,int nlist,t_dlist dlist[],real **dih,
742 gmx_bool bViol,gmx_bool bRamOmega,const output_env_t oenv)
743 {
744 FILE *fp,*gp=NULL;
745 gmx_bool bOm;
746 char fn[256];
747 int i,j,k,Xi1,Xi2,Phi,Psi,Om=0,nlevels;
748 #define NMAT 120
749 real **mat=NULL,phi,psi,omega,axis[NMAT],lo,hi;
750 t_rgb rlo = { 1.0, 0.0, 0.0 };
751 t_rgb rmid= { 1.0, 1.0, 1.0 };
752 t_rgb rhi = { 0.0, 0.0, 1.0 };
753
754 for(i=0; (i<nlist); i++) {
755 if ((has_dihedral(edPhi,&(dlist[i]))) &&
756 (has_dihedral(edPsi,&(dlist[i])))) {
757 sprintf(fn,"ramaPhiPsi%s.xvg",dlist[i].name);
758 fp = rama_file(fn,"Ramachandran Plot",
759 "\\8f\\4 (deg)","\\8y\\4 (deg)",oenv);
760 bOm = bRamOmega && has_dihedral(edOmega,&(dlist[i]));
761 if (bOm) {
762 Om = dlist[i].j0[edOmega];
763 snew(mat,NMAT);
764 for(j=0; (j<NMAT); j++) {
765 snew(mat[j],NMAT);
766 axis[j] = -180+(360*j)/NMAT;
767 }
768 }
769 if (bViol) {
770 sprintf(fn,"violPhiPsi%s.xvg",dlist[i].name);
771 gp = ffopen(fn,"w");
772 }
773 Phi = dlist[i].j0[edPhi];
774 Psi = dlist[i].j0[edPsi];
775 for(j=0; (j<nf); j++) {
776 phi = RAD2DEG*dih[Phi][j];
777 psi = RAD2DEG*dih[Psi][j];
778 fprintf(fp,"%10g %10g\n",phi,psi);
779 if (bViol)
780 fprintf(gp,"%d\n",!bAllowed(dih[Phi][j],RAD2DEG*dih[Psi][j]));
781 if (bOm) {
782 omega = RAD2DEG*dih[Om][j];
783 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
784 += omega;
785 }
786 }
787 if (bViol)
788 ffclose(gp);
789 ffclose(fp);
790 if (bOm) {
791 sprintf(fn,"ramomega%s.xpm",dlist[i].name);
792 fp = ffopen(fn,"w");
793 lo = hi = 0;
794 for(j=0; (j<NMAT); j++)
795 for(k=0; (k<NMAT); k++) {
796 mat[j][k] /= nf;
797 lo=min(mat[j][k],lo);
798 hi=max(mat[j][k],hi);
799 }
800 /* Symmetrise */
801 if (fabs(lo) > fabs(hi))
802 hi = -lo;
803 else
804 lo = -hi;
805 /* Add 180 */
806 for(j=0; (j<NMAT); j++)
807 for(k=0; (k<NMAT); k++)
808 mat[j][k] += 180;
809 lo += 180;
810 hi += 180;
811 nlevels = 20;
812 write_xpm3(fp,0,"Omega/Ramachandran Plot","Deg","Phi","Psi",
813 NMAT,NMAT,axis,axis,mat,lo,180.0,hi,rlo,rmid,rhi,&nlevels);
814 ffclose(fp);
815 for(j=0; (j<NMAT); j++)
816 sfree(mat[j]);
817 sfree(mat);
818 }
819 }
820 if ((has_dihedral(edChi1,&(dlist[i]))) &&
821 (has_dihedral(edChi2,&(dlist[i])))) {
822 sprintf(fn,"ramaX1X2%s.xvg",dlist[i].name);
823 fp = rama_file(fn,"\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
824 "\\8c\\4\\s1\\N (deg)","\\8c\\4\\s2\\N (deg)",oenv);
825 Xi1 = dlist[i].j0[edChi1];
826 Xi2 = dlist[i].j0[edChi2];
827 for(j=0; (j<nf); j++)
828 fprintf(fp,"%10g %10g\n",RAD2DEG*dih[Xi1][j],RAD2DEG*dih[Xi2][j]);
829 ffclose(fp);
830 }
831 else
832 fprintf(stderr,"No chi1 & chi2 angle for %s\n",dlist[i].name);
833 }
834 }
835
836
837 static void print_transitions(const char *fn,int maxchi,int nlist,
838 t_dlist dlist[], t_atoms *atoms,rvec x[],
839 matrix box, gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,real dt,
840 const output_env_t oenv)
841 {
842 /* based on order_params below */
843 FILE *fp;
844 int nh[edMax];
845 int i,Dih,Xi;
846
847 /* must correspond with enum in pp2shift.h:38 */
848 char *leg[edMax];
849 #define NLEG asize(leg)
850
851 leg[0] = strdup("Phi");
852 leg[1] = strdup("Psi");
853 leg[2] = strdup("Omega");
854 leg[3] = strdup("Chi1");
855 leg[4] = strdup("Chi2");
856 leg[5] = strdup("Chi3");
857 leg[6] = strdup("Chi4");
858 leg[7] = strdup("Chi5");
859 leg[8] = strdup("Chi6");
860
861 /* Print order parameters */
862 fp=xvgropen(fn,"Dihedral Rotamer Transitions","Residue","Transitions/ns",
863 oenv);
864 xvgr_legend(fp,NONCHI+maxchi,(const char**)leg,oenv);
865
866 for (Dih=0; (Dih<edMax); Dih++)
867 nh[Dih]=0;
868
869 fprintf(fp,"%5s ","#Res.");
870 fprintf(fp,"%10s %10s %10s ",leg[edPhi],leg[edPsi],leg[edOmega]);
871 for (Xi=0; Xi<maxchi; Xi++)
872 fprintf(fp,"%10s ",leg[NONCHI+Xi]);
873 fprintf(fp,"\n");
874
875 for(i=0; (i<nlist); i++) {
876 fprintf(fp,"%5d ",dlist[i].resnr);
877 for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
878 fprintf(fp,"%10.3f ",dlist[i].ntr[Dih]/dt);
879 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
880 fprintf(fp,"\n");
881 }
882 ffclose(fp);
883 }
884
885 static void order_params(FILE *log,
886 const char *fn,int maxchi,int nlist,t_dlist dlist[],
887 const char *pdbfn,real bfac_init,
888 t_atoms *atoms,rvec x[],int ePBC,matrix box,
889 gmx_bool bPhi,gmx_bool bPsi,gmx_bool bChi,const output_env_t oenv)
890 {
891 FILE *fp;
892 int nh[edMax];
893 char buf[STRLEN];
894 int i,Dih,Xi;
895 real S2Max, S2Min;
896
897 /* except for S2Min/Max, must correspond with enum in pp2shift.h:38 */
898 const char *const_leg[2+edMax]= { "S2Min","S2Max","Phi","Psi","Omega",
899 "Chi1", "Chi2", "Chi3", "Chi4", "Chi5",
900 "Chi6" };
901 #define NLEG asize(leg)
902
903 char *leg[2+edMax];
904
905 for(i=0;i<NLEG;i++)
906 leg[i]=strdup(const_leg[i]);
907
908 /* Print order parameters */
909 fp=xvgropen(fn,"Dihedral Order Parameters","Residue","S2",oenv);
910 xvgr_legend(fp,2+NONCHI+maxchi,const_leg,oenv);
911
912 for (Dih=0; (Dih<edMax); Dih++)
913 nh[Dih]=0;
914
915 fprintf(fp,"%5s ","#Res.");
916 fprintf(fp,"%10s %10s ",leg[0],leg[1]);
917 fprintf(fp,"%10s %10s %10s ",leg[2+edPhi],leg[2+edPsi],leg[2+edOmega]);
918 for (Xi=0; Xi<maxchi; Xi++)
919 fprintf(fp,"%10s ",leg[2+NONCHI+Xi]);
920 fprintf(fp,"\n");
921
922 for(i=0; (i<nlist); i++) {
923 S2Max=-10;
924 S2Min=10;
925 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
926 if (dlist[i].S2[Dih]!=0) {
927 if (dlist[i].S2[Dih] > S2Max)
928 S2Max=dlist[i].S2[Dih];
929 if (dlist[i].S2[Dih] < S2Min)
930 S2Min=dlist[i].S2[Dih];
931 }
932 if (dlist[i].S2[Dih] > 0.8)
933 nh[Dih]++;
934 }
935 fprintf(fp,"%5d ",dlist[i].resnr);
936 fprintf(fp,"%10.3f %10.3f ",S2Min,S2Max);
937 for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
938 fprintf(fp,"%10.3f ",dlist[i].S2[Dih]);
939 fprintf(fp,"\n");
940 /* fprintf(fp,"%12s\n",dlist[i].name); this confuses xmgrace */
941 }
942 ffclose(fp);
943
944 if (NULL != pdbfn) {
945 real x0,y0,z0;
946
947 if (NULL == atoms->pdbinfo)
948 snew(atoms->pdbinfo,atoms->nr);
949 for(i=0; (i<atoms->nr); i++)
950 atoms->pdbinfo[i].bfac=bfac_init;
951
952 for(i=0; (i<nlist); i++) {
953 atoms->pdbinfo[dlist[i].atm.N].bfac=-dlist[i].S2[0];/* Phi */
954 atoms->pdbinfo[dlist[i].atm.H].bfac=-dlist[i].S2[0];/* Phi */
955 atoms->pdbinfo[dlist[i].atm.C].bfac=-dlist[i].S2[1];/* Psi */
956 atoms->pdbinfo[dlist[i].atm.O].bfac=-dlist[i].S2[1];/* Psi */
957 for (Xi=0; (Xi<maxchi); Xi++) { /* Chi's */
958 if (dlist[i].atm.Cn[Xi+3]!=-1) {
959 atoms->pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac=-dlist[i].S2[NONCHI+Xi];
960 }
961 }
962 }
963
964 fp=ffopen(pdbfn,"w");
965 fprintf(fp,"REMARK generated by g_chi\n");
966 fprintf(fp,"REMARK "
967 "B-factor field contains negative of dihedral order parameters\n");
968 write_pdbfile(fp,NULL,atoms,x,ePBC,box,' ',0,NULL,TRUE);
969 x0=y0=z0=1000.0;
970 for (i=0; (i<atoms->nr); i++) {
971 x0 = min(x0, x[i][XX]);
972 y0 = min(y0, x[i][YY]);
973 z0 = min(z0, x[i][ZZ]);
974 }
975 x0*=10.0;/* nm -> angstrom */
976 y0*=10.0;/* nm -> angstrom */
977 z0*=10.0;/* nm -> angstrom */
978 sprintf(buf,"%s%%6.f%%6.2f\n",pdbformat);
979 for (i=0; (i<10); i++) {
980 fprintf(fp,buf,"ATOM ", atoms->nr+1+i, "CA", "LEG",' ',
981 atoms->nres+1, ' ',x0, y0, z0+(1.2*i), 0.0, -0.1*i);
982 }
983 ffclose(fp);
984 }
985
986 fprintf(log,"Dihedrals with S2 > 0.8\n");
987 fprintf(log,"Dihedral: ");
988 if (bPhi) fprintf(log," Phi ");
989 if (bPsi) fprintf(log," Psi ");
990 if (bChi)
991 for(Xi=0; (Xi<maxchi); Xi++)
992 fprintf(log," %s ", leg[2+NONCHI+Xi]);
993 fprintf(log,"\nNumber: ");
994 if (bPhi) fprintf(log,"%4d ",nh[0]);
995 if (bPsi) fprintf(log,"%4d ",nh[1]);
996 if (bChi)
997 for(Xi=0; (Xi<maxchi); Xi++)
998 fprintf(log,"%4d ",nh[NONCHI+Xi]);
999 fprintf(log,"\n");
1000
1001 for(i=0;i<NLEG;i++)
1002 sfree(leg[i]);
1003
1004 }
1005
1006 int gmx_chi(int argc,char *argv[])
1007 {
1008 const char *desc[] = {
1009 "g_chi computes phi, psi, omega and chi dihedrals for all your ",
1010 "amino acid backbone and sidechains.",
1011 "It can compute dihedral angle as a function of time, and as",
1012 "histogram distributions.",
1013 "The distributions (histo-(dihedral)(RESIDUE).xvg) are cumulative over all residues of each type.[PAR]",
1014 "If option [TT]-corr[tt] is given, the program will",
1015 "calculate dihedral autocorrelation functions. The function used",
1016 "is C(t) = < cos(chi(tau)) cos(chi(tau+t)) >. The use of cosines",
1017 "rather than angles themselves, resolves the problem of periodicity.",
1018 "(Van der Spoel & Berendsen (1997), [BB]Biophys. J. 72[bb], 2032-2041).",
1019 "Separate files for each dihedral of each residue",
1020 "(corr(dihedral)(RESIDUE)(nresnr).xvg) are output, as well as a",
1021 "file containing the information for all residues (argument of [TT]-corr[tt]).[PAR]",
1022 "With option [TT]-all[tt], the angles themselves as a function of time for",
1023 "each residue are printed to separate files (dihedral)(RESIDUE)(nresnr).xvg.",
1024 "These can be in radians or degrees.[PAR]",
1025 "A log file (argument [TT]-g[tt]) is also written. This contains [BR]",
1026 "(a) information about the number of residues of each type.[BR]",
1027 "(b) The NMR 3J coupling constants from the Karplus equation.[BR]",
1028 "(c) a table for each residue of the number of transitions between ",
1029 "rotamers per nanosecond, and the order parameter S2 of each dihedral.[BR]",
1030 "(d) a table for each residue of the rotamer occupancy.[BR]",
1031 "All rotamers are taken as 3-fold, except for omegas and chi-dihedrals",
1032 "to planar groups (i.e. chi2 of aromatics asp and asn, chi3 of glu",
1033 "and gln, and chi4 of arg), which are 2-fold. \"rotamer 0\" means ",
1034 "that the dihedral was not in the core region of each rotamer. ",
1035 "The width of the core region can be set with [TT]-core_rotamer[tt][PAR]",
1036
1037 "The S2 order parameters are also output to an xvg file",
1038 "(argument [TT]-o[tt] ) and optionally as a pdb file with",
1039 "the S2 values as B-factor (argument [TT]-p[tt]). ",
1040 "The total number of rotamer transitions per timestep",
1041 "(argument [TT]-ot[tt]), the number of transitions per rotamer",
1042 "(argument [TT]-rt[tt]), and the 3J couplings (argument [TT]-jc[tt]), ",
1043 "can also be written to .xvg files.[PAR]",
1044
1045 "If [TT]-chi_prod[tt] is set (and maxchi > 0), cumulative rotamers, e.g.",
1046 "1+9(chi1-1)+3(chi2-1)+(chi3-1) (if the residue has three 3-fold ",
1047 "dihedrals and maxchi >= 3)",
1048 "are calculated. As before, if any dihedral is not in the core region,",
1049 "the rotamer is taken to be 0. The occupancies of these cumulative ",
1050 "rotamers (starting with rotamer 0) are written to the file",
1051 "that is the argument of [TT]-cp[tt], and if the [TT]-all[tt] flag",
1052 "is given, the rotamers as functions of time",
1053 "are written to chiproduct(RESIDUE)(nresnr).xvg ",
1054 "and their occupancies to histo-chiproduct(RESIDUE)(nresnr).xvg.[PAR]",
1055
1056 "The option [TT]-r[tt] generates a contour plot of the average omega angle",
1057 "as a function of the phi and psi angles, that is, in a Ramachandran plot",
1058 "the average omega angle is plotted using color coding.",
1059
1060 };
1061
1062 const char *bugs[] = {
1063 "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.",
1064 "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.",
1065 "-r0 option does not work properly",
1066 "Rotamers with multiplicity 2 are printed in chi.log as if they had multiplicity 3, with the 3rd (g(+)) always having probability 0"
1067 };
1068
1069 /* defaults */
1070 static int r0=1,ndeg=1,maxchi=2;
1071 static gmx_bool bAll=FALSE;
1072 static gmx_bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
1073 static real bfac_init=-1.0,bfac_max=0;
1074 static const char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
1075 static gmx_bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
1076 static gmx_bool bNormHisto=TRUE,bChiProduct=FALSE,bHChi=FALSE,bRAD=FALSE,bPBC=TRUE;
1077 static real core_frac=0.5 ;
1078 t_pargs pa[] = {
1079 { "-r0", FALSE, etINT, {&r0},
1080 "starting residue" },
1081 { "-phi", FALSE, etBOOL, {&bPhi},
1082 "Output for Phi dihedral angles" },
1083 { "-psi", FALSE, etBOOL, {&bPsi},
1084 "Output for Psi dihedral angles" },
1085 { "-omega",FALSE, etBOOL, {&bOmega},
1086 "Output for Omega dihedrals (peptide bonds)" },
1087 { "-rama", FALSE, etBOOL, {&bRama},
1088 "Generate Phi/Psi and Chi1/Chi2 ramachandran plots" },
1089 { "-viol", FALSE, etBOOL, {&bViol},
1090 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
1091 { "-periodic", FALSE, etBOOL, {&bPBC},
1092 "Print dihedral angles modulo 360 degrees" },
1093 { "-all", FALSE, etBOOL, {&bAll},
1094 "Output separate files for every dihedral." },
1095 { "-rad", FALSE, etBOOL, {&bRAD},
1096 "in angle vs time files, use radians rather than degrees."},
1097 { "-shift", FALSE, etBOOL, {&bShift},
1098 "Compute chemical shifts from Phi/Psi angles" },
1099 { "-binwidth", FALSE, etINT, {&ndeg},
1100 "bin width for histograms (degrees)" },
1101 { "-core_rotamer", FALSE, etREAL, {&core_frac},
1102 "only the central -core_rotamer*(360/multiplicity) belongs to each rotamer (the rest is assigned to rotamer 0)" },
1103 { "-maxchi", FALSE, etENUM, {maxchistr},
1104 "calculate first ndih Chi dihedrals" },
1105 { "-normhisto", FALSE, etBOOL, {&bNormHisto},
1106 "Normalize histograms" },
1107 { "-ramomega",FALSE,etBOOL, {&bRamOmega},
1108 "compute average omega as a function of phi/psi and plot it in an xpm plot" },
1109 { "-bfact", FALSE, etREAL, {&bfac_init},
1110 "B-factor value for pdb file for atoms with no calculated dihedral order parameter"},
1111 { "-chi_prod",FALSE,etBOOL, {&bChiProduct},
1112 "compute a single cumulative rotamer for each residue"},
1113 { "-HChi",FALSE,etBOOL, {&bHChi},
1114 "Include dihedrals to sidechain hydrogens"},
1115 { "-bmax", FALSE, etREAL, {&bfac_max},
1116 "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." }
1117 };
1118
1119 FILE *log;
1120 int natoms,nlist,idum,nbin;
1121 t_atoms atoms;
1122 rvec *x;
1123 int ePBC;
1124 matrix box;
1125 char title[256],grpname[256];
1126 t_dlist *dlist;
1127 gmx_bool bChi,bCorr,bSSHisto;
1128 gmx_bool bDo_rt, bDo_oh, bDo_ot, bDo_jc ;
1129 real dt=0, traj_t_ns;
1130 output_env_t oenv;
1131 gmx_residuetype_t rt;
1132
1133 atom_id isize,*index;
1134 int ndih,nactdih,nf;
1135 real **dih,*trans_frac,*aver_angle,*time;
1136 int i,j,**chi_lookup,*xity;
1137
1138 t_filenm fnm[] = {
1139 { efSTX, "-s", NULL, ffREAD },
1140 { efTRX, "-f", NULL, ffREAD },
1141 { efXVG, "-o", "order", ffWRITE },
1142 { efPDB, "-p", "order", ffOPTWR },
1143 { efDAT, "-ss", "ssdump", ffOPTRD },
1144 { efXVG, "-jc", "Jcoupling", ffWRITE },
1145 { efXVG, "-corr", "dihcorr",ffOPTWR },
1146 { efLOG, "-g", "chi", ffWRITE },
1147 /* add two more arguments copying from g_angle */
1148 { efXVG, "-ot", "dihtrans", ffOPTWR },
1149 { efXVG, "-oh", "trhisto", ffOPTWR },
1150 { efXVG, "-rt", "restrans", ffOPTWR },
1151 { efXVG, "-cp", "chiprodhisto", ffOPTWR }
1152 };
1153 #define NFILE asize(fnm)
1154 int npargs;
1155 t_pargs *ppa;
1156
1157 CopyRight(stderr,argv[0]);
1158 npargs = asize(pa);
1159 ppa = add_acf_pargs(&npargs,pa);
1160 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
1161 NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs,
1162 &oenv);
1163
1164 /* Handle result from enumerated type */
1165 sscanf(maxchistr[0],"%d",&maxchi);
1166 bChi = (maxchi > 0);
1167
1168 log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
1169
1170 if (bRamOmega) {
1171 bOmega = TRUE;
1172 bPhi = TRUE;
1173 bPsi = TRUE;
1174 }
1175
1176 /* set some options */
1177 bDo_rt=(opt2bSet("-rt",NFILE,fnm));
1178 bDo_oh=(opt2bSet("-oh",NFILE,fnm));
1179 bDo_ot=(opt2bSet("-ot",NFILE,fnm));
1180 bDo_jc=(opt2bSet("-jc",NFILE,fnm));
1181 bCorr=(opt2bSet("-corr",NFILE,fnm));
1182 if (bCorr)
1183 fprintf(stderr,"Will calculate autocorrelation\n");
1184
1185 if (core_frac > 1.0 ) {
1186 fprintf(stderr, "core_rotamer fraction > 1.0 ; will use 1.0\n");
1187 core_frac=1.0 ;
1188 }
1189 if (core_frac < 0.0 ) {
1190 fprintf(stderr, "core_rotamer fraction < 0.0 ; will use 0.0\n");
1191 core_frac=0.0 ;
1192 }
1193
1194 if (maxchi > MAXCHI) {
1195 fprintf(stderr,
1196 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
1197 MAXCHI, maxchi);
1198 maxchi=MAXCHI;
1199 }
1200 bSSHisto = ftp2bSet(efDAT,NFILE,fnm);
1201 nbin = 360/ndeg ;
1202
1203 /* Find the chi angles using atoms struct and a list of amino acids */
1204 get_stx_coordnum(ftp2fn(efSTX,NFILE,fnm),&natoms);
1205 init_t_atoms(&atoms,natoms,TRUE);
1206 snew(x,natoms);
1207 read_stx_conf(ftp2fn(efSTX,NFILE,fnm),title,&atoms,x,NULL,&ePBC,box);
1208 fprintf(log,"Title: %s\n",title);
1209
1210 gmx_residuetype_init(&rt);
1211 dlist=mk_dlist(log,&atoms,&nlist,bPhi,bPsi,bChi,bHChi,maxchi,r0,rt);
1212 fprintf(stderr,"%d residues with dihedrals found\n", nlist);
1213
1214 if (nlist == 0)
1215 gmx_fatal(FARGS,"No dihedrals in your structure!\n");
1216
1217 /* Make a linear index for reading all. */
1218 index=make_chi_ind(nlist,dlist,&ndih);
1219 isize=4*ndih;
1220 fprintf(stderr,"%d dihedrals found\n", ndih);
1221
1222 snew(dih,ndih);
1223
1224 /* COMPUTE ALL DIHEDRALS! */
1225 read_ang_dih(ftp2fn(efTRX,NFILE,fnm),FALSE,TRUE,FALSE,bPBC,1,&idum,
1226 &nf,&time,isize,index,&trans_frac,&aver_angle,dih,oenv);
1227
1228 dt=(time[nf-1]-time[0])/(nf-1); /* might want this for corr or n. transit*/
1229 if (bCorr)
1230 {
1231 if (nf < 2)
1232 {
1233 gmx_fatal(FARGS,"Need at least 2 frames for correlation");
1234 }
1235 }
1236
1237 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi
1238 * pass nactdih instead of ndih to low_ana_dih_trans and get_chi_product_traj
1239 * to prevent accessing off end of arrays when maxchi < 5 or 6. */
1240 nactdih = reset_em_all(nlist,dlist,nf,dih,maxchi);
1241
1242 if (bAll)
1243 dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega,bRAD,oenv);
1244
1245 /* Histogramming & J coupling constants & calc of S2 order params */
1246 histogramming(log,nbin,rt,nf,maxchi,dih,nlist,dlist,index,
1247 bPhi,bPsi,bOmega,bChi,
1248 bNormHisto,bSSHisto,ftp2fn(efDAT,NFILE,fnm),bfac_max,&atoms,
1249 bDo_jc,opt2fn("-jc",NFILE,fnm),oenv);
1250
1251 /* transitions
1252 *
1253 * added multiplicity */
1254
1255 snew(xity,ndih) ;
1256 mk_multiplicity_lookup(xity, maxchi, dih, nlist, dlist,ndih);
1257
1258 strcpy(grpname, "All residues, ");
1259 if(bPhi)
1260 strcat(grpname, "Phi ");
1261 if(bPsi)
1262 strcat(grpname, "Psi ");
1263 if(bOmega)
1264 strcat(grpname, "Omega ");
1265 if(bChi){
1266 strcat(grpname, "Chi 1-") ;
1267 sprintf(grpname + strlen(grpname), "%i", maxchi);
1268 }
1269
1270
1271 low_ana_dih_trans(bDo_ot, opt2fn("-ot",NFILE,fnm),
1272 bDo_oh, opt2fn("-oh",NFILE,fnm),maxchi,
1273 dih, nlist, dlist, nf, nactdih, grpname, xity,
1274 *time, dt, FALSE, core_frac,oenv) ;
1275
1276 /* Order parameters */
1277 order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
1278 ftp2fn_null(efPDB,NFILE,fnm),bfac_init,
1279 &atoms,x,ePBC,box,bPhi,bPsi,bChi,oenv);
1280
1281 /* Print ramachandran maps! */
1282 if (bRama)
1283 do_rama(nf,nlist,dlist,dih,bViol,bRamOmega,oenv);
1284
1285 if (bShift)
1286 do_pp2shifts(log,nf,nlist,dlist,dih);
1287
1288 /* rprint S^2, transitions, and rotamer occupancies to log */
1289 traj_t_ns = 0.001 * (time[nf-1]-time[0]) ;
1290 pr_dlist(log,nlist,dlist,traj_t_ns,edPrintST,bPhi,bPsi,bChi,bOmega,maxchi);
1291 pr_dlist(log,nlist,dlist,traj_t_ns,edPrintRO,bPhi,bPsi,bChi,bOmega,maxchi);
1292 ffclose(log);
1293 /* transitions to xvg */
1294 if (bDo_rt)
1295 print_transitions(opt2fn("-rt",NFILE,fnm),maxchi,nlist,dlist,
1296 &atoms,x,box,bPhi,bPsi,bChi,traj_t_ns,oenv);
1297
1298 /* chi_product trajectories (ie one "rotamer number" for each residue) */
1299 if (bChiProduct && bChi ) {
1300 snew(chi_lookup,nlist) ;
1301 for (i=0;i<nlist;i++)
1302 snew(chi_lookup[i], maxchi) ;
1303 mk_chi_lookup(chi_lookup, maxchi, dih, nlist, dlist);
1304
1305 get_chi_product_traj(dih,nf,nactdih,nlist,
1306 maxchi,dlist,time,chi_lookup,xity,
1307 FALSE,bNormHisto, core_frac,bAll,
1308 opt2fn("-cp",NFILE,fnm),oenv);
1309
1310 for (i=0;i<nlist;i++)
1311 sfree(chi_lookup[i]);
1312 }
1313
1314 /* Correlation comes last because it fucks up the angles */
1315 if (bCorr)
1316 do_dihcorr(opt2fn("-corr",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
1317 maxchi,bPhi,bPsi,bChi,bOmega,oenv);
1318
1319
1320 do_view(oenv,opt2fn("-o",NFILE,fnm),"-nxy");
1321 do_view(oenv,opt2fn("-jc",NFILE,fnm),"-nxy");
1322 if (bCorr)
1323 do_view(oenv,opt2fn("-corr",NFILE,fnm),"-nxy");
1324
1325 gmx_residuetype_destroy(rt);
1326
1327 thanx(stderr);
1328
1329 return 0;
1330 }
got merge a long time ago