?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / g_chi.c
blob3e162390aca7ac02f45f11a5b54e5b31aeace358
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_g_chi_c = "$Id$";
30
31 #include <stdio.h>
32 #include <math.h>
33 #include "confio.h"
34 #include "pdbio.h"
35 #include "copyrite.h"
36 #include "fatal.h"
37 #include "futil.h"
38 #include "gstat.h"
39 #include "macros.h"
40 #include "maths.h"
41 #include "physics.h"
42 #include "rdgroup.h"
43 #include "smalloc.h"
44 #include "statutil.h"
45 #include "tpxio.h"
46 #include "string.h"
47 #include "sysstuff.h"
48 #include "txtdump.h"
49 #include "typedefs.h"
50 #include "vec.h"
51 #include "strdb.h"
52 #include "xvgr.h"
53 #include "pp2shift.h"
54 #include "matio.h"
55
56 static bool bAllowed(real phi,real psi)
57 {
58 static char *map[] = {
59 "1100000000000000001111111000000000001111111111111111111111111",
60 "1100000000000000001111110000000000011111111111111111111111111",
61 "1100000000000000001111110000000000011111111111111111111111111",
62 "1100000000000000001111100000000000111111111111111111111111111",
63 "1100000000000000001111100000000000111111111111111111111111111",
64 "1100000000000000001111100000000001111111111111111111111111111",
65 "1100000000000000001111100000000001111111111111111111111111111",
66 "1100000000000000001111100000000011111111111111111111111111111",
67 "1110000000000000001111110000000111111111111111111111111111111",
68 "1110000000000000001111110000001111111111111111111111111111111",
69 "1110000000000000001111111000011111111111111111111111111111111",
70 "1110000000000000001111111100111111111111111111111111111111111",
71 "1110000000000000001111111111111111111111111111111111111111111",
72 "1110000000000000001111111111111111111111111111111111111111111",
73 "1110000000000000001111111111111111111111111111111111111111111",
74 "1110000000000000001111111111111111111111111111111111111111111",
75 "1110000000000000001111111111111110011111111111111111111111111",
76 "1110000000000000001111111111111100000111111111111111111111111",
77 "1110000000000000001111111111111000000000001111111111111111111",
78 "1100000000000000001111111111110000000000000011111111111111111",
79 "1100000000000000001111111111100000000000000011111111111111111",
80 "1000000000000000001111111111000000000000000001111111111111110",
81 "0000000000000000001111111110000000000000000000111111111111100",
82 "0000000000000000000000000000000000000000000000000000000000000",
83 "0000000000000000000000000000000000000000000000000000000000000",
84 "0000000000000000000000000000000000000000000000000000000000000",
85 "0000000000000000000000000000000000000000000000000000000000000",
86 "0000000000000000000000000000000000000000000000000000000000000",
87 "0000000000000000000000000000000000000000000000000000000000000",
88 "0000000000000000000000000000000000000000000000000000000000000",
89 "0000000000000000000000000000000000000000000000000000000000000",
90 "0000000000000000000000000000000000000000000000000000000000000",
91 "0000000000000000000000000000000000000000000000000000000000000",
92 "0000000000000000000000000000000000000000000000000000000000000",
93 "0000000000000000000000000000000000000000000000000000000000000",
94 "0000000000000000000000000000000000000000000000000000000000000",
95 "0000000000000000000000000000000000000000000000000000000000000",
96 "0000000000000000000000000000000000000000000000000000000000000",
97 "0000000000000000000000000000000000111111111111000000000000000",
98 "1100000000000000000000000000000001111111111111100000000000111",
99 "1100000000000000000000000000000001111111111111110000000000111",
100 "0000000000000000000000000000000000000000000000000000000000000",
101 "0000000000000000000000000000000000000000000000000000000000000",
102 "0000000000000000000000000000000000000000000000000000000000000",
103 "0000000000000000000000000000000000000000000000000000000000000",
104 "0000000000000000000000000000000000000000000000000000000000000",
105 "0000000000000000000000000000000000000000000000000000000000000",
106 "0000000000000000000000000000000000000000000000000000000000000",
107 "0000000000000000000000000000000000000000000000000000000000000",
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 };
121 #define NPP asize(map)
122 int x,y;
123
124 #define INDEX(ppp) ((((int) (360+ppp*RAD2DEG)) % 360)/6)
125 x = INDEX(phi);
126 y = INDEX(psi);
127 #undef INDEX
128 return (bool) map[x][y];
129 }
130
131 atom_id *make_chi_ind(int nl,t_dlist dl[],int *ndih)
132 {
133 atom_id *id;
134 int i,Xi,n;
135
136 /* There are nl sidechains with max edMax dihedrals with 4 atoms each */
137 snew(id,nl*edMax*4);
138
139 n=0;
140 for(i=0; (i<nl); i++) { /* Phi */
141 dl[i].j0[edPhi] = n/4;
142 if (dl[i].atm.H == -1)
143 id[n++]=dl[i].atm.minC;
144 else
145 id[n++]=dl[i].atm.H;
146 id[n++]=dl[i].atm.N;
147 id[n++]=dl[i].atm.Cn[1];
148 id[n++]=dl[i].atm.C;
149 }
150 for(i=0; (i<nl); i++) { /* Psi */
151 dl[i].j0[edPsi] = n/4;
152 id[n++]=dl[i].atm.N;
153 id[n++]=dl[i].atm.Cn[1];
154 id[n++]=dl[i].atm.C;
155 id[n++]=dl[i].atm.O;
156 }
157 for(i=0; (i<nl); i++) { /* Omega */
158 if (has_dihedral(edOmega,&(dl[i]))) {
159 dl[i].j0[edOmega] = n/4;
160 id[n++]=dl[i].atm.minO;
161 id[n++]=dl[i].atm.minC;
162 id[n++]=dl[i].atm.N;
163 id[n++]=dl[i].atm.Cn[1];
164 }
165 }
166 for(Xi=0; (Xi<MAXCHI); Xi++) { /* Chi# */
167 for(i=0; (i<nl); i++) {
168 if (dl[i].atm.Cn[Xi+3] != -1) {
169 dl[i].j0[edChi1+Xi] = n/4;
170 id[n++]=dl[i].atm.Cn[Xi];
171 id[n++]=dl[i].atm.Cn[Xi+1];
172 id[n++]=dl[i].atm.Cn[Xi+2];
173 id[n++]=dl[i].atm.Cn[Xi+3];
174 }
175 }
176 }
177 *ndih=n/4;
178
179 return id;
180 }
181
182 int bin(real chi,int mult)
183 {
184 mult=3;
185
186 return (int) (chi*mult/360.0);
187 }
188
189 static void print_one(char *base,char *name,char *title,
190 int nf,real time[],real data[])
191 {
192 FILE *fp;
193 char buf[256],t2[256];
194 int k;
195
196 sprintf(buf,"%s%s.xvg",base,name);
197 fprintf(stderr,"\rPrinting %s ",buf);
198 sprintf(t2,"%s %s",title,name);
199 fp=xvgropen(buf,t2,"Time (ps)","C(t)");
200 for(k=0; (k<nf); k++)
201 fprintf(fp,"%10g %10g\n",time[k],data[k]);
202 ffclose(fp);
203 }
204
205 static void do_dihcorr(char *fn,int nf,int ndih,real **dih,real dt,
206 int nlist,t_dlist dlist[],real time[],int maxchi,
207 bool bPhi,bool bPsi,bool bChi,bool bOmega)
208 {
209 char name1[256],name2[256];
210 int i,j,Xi;
211
212 do_autocorr(fn,"Dihedral Autocorrelation Function",
213 nf,ndih,dih,dt,eacCos,FALSE);
214 /* Dump em all */
215 j=0;
216 for(i=0; (i<nlist); i++) {
217 if (bPhi)
218 print_one("corrphi",dlist[i].name,"Phi ACF for",nf/2,time,dih[j]);
219 j++;
220 }
221 for(i=0; (i<nlist); i++) {
222 if (bPsi)
223 print_one("corrpsi",dlist[i].name,"Psi ACF for",nf/2,time,dih[j]);
224 j++;
225 }
226 for(i=0; (i<nlist); i++) {
227 if (has_dihedral(edOmega,&dlist[i])) {
228 if (bOmega)
229 print_one("corromega",dlist[i].name,"Omega ACF for",nf/2,time,dih[j]);
230 j++;
231 }
232 }
233 for(Xi=0; (Xi<maxchi); Xi++) {
234 sprintf(name1, "corrchi%d", Xi+1);
235 sprintf(name2, "Chi%d ACF for", Xi+1);
236 for(i=0; (i<nlist); i++) {
237 if (dlist[i].atm.Cn[Xi+3] != -1) {
238 if (bChi)
239 print_one(name1,dlist[i].name,name2,nf/2,time,dih[j]);
240 j++;
241 }
242 }
243 }
244 fprintf(stderr,"\n");
245 }
246
247 static void dump_em_all(int nlist,t_dlist dlist[],int nf,real time[],
248 real **dih,int maxchi,
249 bool bPhi,bool bPsi,bool bChi,bool bOmega)
250 {
251 char name[256];
252 int i,j,Xi;
253
254 /* Dump em all */
255 j = 0;
256 for(i=0; (i<nlist); i++) {
257 if (bPhi)
258 print_one("phi",dlist[i].name,name,nf,time,dih[j]);
259 j++;
260 }
261 for(i=0; (i<nlist); i++) {
262 if (bPsi)
263 print_one("psi",dlist[i].name,name,nf,time,dih[j]);
264 j++;
265 }
266 for(i=0; (i<nlist); i++)
267 if (has_dihedral(edOmega,&(dlist[i]))) {
268 if (bOmega)
269 print_one("omega",dlist[i].name,name,nf,time,dih[j]);
270 j++;
271 }
272
273 for(Xi=0; (Xi<maxchi); Xi++)
274 for(i=0; (i<nlist); i++)
275 if (dlist[i].atm.Cn[Xi+3] != -1) {
276 if (bChi) {
277 sprintf(name,"chi%d",Xi+1);
278 print_one(name,dlist[i].name,name,nf,time,dih[j]);
279 }
280 j++;
281 }
282 fprintf(stderr,"\n");
283 }
284
285 static void reset_one(real dih[],int nf,real phase)
286 {
287 int j;
288
289 for(j=0; (j<nf); j++) {
290 dih[j] += phase;
291 while(dih[j] < -M_PI)
292 dih[j] += 2*M_PI;
293 while(dih[j] >= M_PI)
294 dih[j] -= 2*M_PI;
295 }
296 }
297
298 static void reset_em_all(int nlist,t_dlist dlist[],int nf,
299 real **dih,int maxchi,bool bPhi,bool bPsi,bool bChi)
300 {
301 int i,j,Xi;
302
303 /* Reset em all */
304 j=0;
305 /* Phi */
306 for(i=0; (i<nlist); i++)
307 if (dlist[i].atm.H == -1)
308 reset_one(dih[j++],nf,0);
309 else
310 reset_one(dih[j++],nf,M_PI);
311 /* Psi */
312 for(i=0; (i<nlist); i++)
313 reset_one(dih[j++],nf,M_PI);
314 /* Omega */
315 for(i=0; (i<nlist); i++)
316 if (has_dihedral(edOmega,&dlist[i]))
317 reset_one(dih[j++],nf,0);
318 /* Chi 1 thru maxchi */
319 for(Xi=0; (Xi<maxchi); Xi++)
320 for(i=0; (i<nlist); i++)
321 if (dlist[i].atm.Cn[Xi+3] != -1) {
322 reset_one(dih[j],nf,0);
323 j++;
324 }
325 fprintf(stderr,"j after resetting = %d\n",j);
326 }
327
328 static void histogramming(FILE *log,int naa,char **aa,
329 int nf,int maxchi,real **dih,
330 int nlist,t_dlist dlist[],
331 bool bPhi,bool bPsi,bool bOmega,bool bChi)
332 {
333 t_karplus kkkphi[] = {
334 { "J_NHa", 6.51, -1.76, 1.6, -M_PI/3, 0.0 },
335 { "J_HaC'", 4.0, 1.1, 0.1, 0.0, 0.0 },
336 { "J_NHCb", 4.7, -1.5, -0.2, M_PI/3, 0.0 },
337 { "J_Ci-1Hai", 4.5, -1.3, -1.2, 2*M_PI/3, 0.0 }
338 };
339 t_karplus kkkpsi[] = {
340 { "J_HaN", -0.88, -0.61,-0.27,M_PI/3, 0.0 }
341 };
342 t_karplus kkkchi1[] = {
343 { "JHaHb2", 9.5, -1.6, 1.8, -M_PI/3, 0 },
344 { "JHaHb3", 9.5, -1.6, 1.8, 0, 0 }
345 };
346 #define NKKKPHI asize(kkkphi)
347 #define NKKKPSI asize(kkkpsi)
348 #define NKKKCHI asize(kkkchi1)
349 #define NJC (NKKKPHI+NKKKPSI+NKKKCHI)
350
351 FILE *fp;
352 real S2;
353 real *normhisto;
354 real **Jc;
355 int ***his_aa,**his_aa1,*histmp;
356 int i,j,k,m,Dih;
357 char hisfile[256],title[256];
358
359 snew(his_aa,edMax);
360 for(Dih=0; (Dih<edMax); Dih++) {
361 snew(his_aa[Dih],naa+1);
362 for(i=0; (i<=naa); i++) {
363 snew(his_aa[Dih][i],NHISTO);
364 }
365 }
366 snew(histmp,NHISTO);
367
368 snew(Jc,nlist);
369 for(i=0; (i<nlist); i++)
370 snew(Jc[i],NJC);
371
372 j=0;
373 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
374 for(i=0; (i<nlist); i++) {
375 if (((Dih < edOmega) ) ||
376 ((Dih == edOmega) && (has_dihedral(edOmega,&(dlist[i])))) ||
377 ((Dih > edOmega) && (dlist[i].atm.Cn[Dih-NONCHI+3] != -1))) {
378 make_histo(log,nf,dih[j],NHISTO,histmp,-M_PI,M_PI);
379
380 switch (Dih) {
381 case edPhi:
382 calc_distribution_props(NHISTO,histmp,-M_PI,NKKKPHI,kkkphi,&S2);
383
384 for(m=0; (m<NKKKPHI); m++)
385 Jc[i][m] = kkkphi[m].Jc;
386 break;
387 case edPsi:
388 calc_distribution_props(NHISTO,histmp,-M_PI,NKKKPSI,kkkpsi,&S2);
389
390 for(m=0; (m<NKKKPSI); m++)
391 Jc[i][NKKKPHI+m] = kkkpsi[m].Jc;
392 break;
393 case edOmega:
394 calc_distribution_props(NHISTO,histmp,-M_PI,0,NULL,&S2);
395 break;
396 case edChi1:
397 calc_distribution_props(NHISTO,histmp,-M_PI,NKKKCHI,kkkchi1,&S2);
398 for(m=0; (m<NKKKCHI); m++)
399 Jc[i][NKKKPHI+NKKKPSI+m] = kkkchi1[m].Jc;
400 break;
401 }
402 dlist[i].S2[Dih] = S2;
403
404 /* Sum distribution per amino acid type as well */
405 for(k=0; (k<NHISTO); k++) {
406 his_aa[Dih][dlist[i].index][k] += histmp[k];
407 histmp[k] = 0;
408 }
409 j++;
410 }
411 }
412 }
413 sfree(histmp);
414
415 /* Print out Jcouplings */
416 fprintf(log,"\n *** J-Couplings from simulation ***\n\n");
417 fprintf(log,"Residue ");
418 for(i=0; (i<NKKKPHI); i++)
419 fprintf(log,"%10s",kkkphi[i].name);
420 for(i=0; (i<NKKKPSI); i++)
421 fprintf(log,"%10s",kkkpsi[i].name);
422 for(i=0; (i<NKKKCHI); i++)
423 fprintf(log,"%10s",kkkchi1[i].name);
424 fprintf(log,"\n");
425 for(i=0; (i<NJC+1); i++)
426 fprintf(log,"----------");
427 fprintf(log,"\n");
428 for(i=0; (i<nlist); i++) {
429 fprintf(log,"%-10s",dlist[i].name);
430 for(j=0; (j<NJC); j++)
431 fprintf(log," %8.3f",Jc[i][j]);
432 fprintf(log,"\n");
433 }
434 fprintf(log,"\n");
435
436 snew(normhisto,NHISTO);
437 for(i=0; (i<naa); i++) {
438 for(Dih=0; (Dih<edMax); Dih++){
439 /* First check whether something is in there */
440 for(j=0; (j<NHISTO); j++)
441 if (his_aa[Dih][i][j] != 0)
442 break;
443 if ((j < NHISTO) &&
444 ((bPhi && (Dih==edPhi)) ||
445 (bPsi && (Dih==edPsi)) ||
446 (bOmega &&(Dih==edOmega)) ||
447 (bChi && (Dih>=edChi1)))) {
448 normalize_histo(NHISTO,his_aa[Dih][i],(360.0/NHISTO),normhisto);
449
450 switch (Dih) {
451 case edPhi:
452 sprintf(hisfile,"histo-phi%s.xvg",aa[i]);
453 sprintf(title,"\\8f\\4 Distribution for %s",aa[i]);
454 break;
455 case edPsi:
456 sprintf(hisfile,"histo-psi%s.xvg",aa[i]);
457 sprintf(title,"\\8y\\4 Distribution for %s",aa[i]);
458 break;
459 case edOmega:
460 sprintf(hisfile,"histo-omega%s.xvg",aa[i]);
461 sprintf(title,"\\8w\\4 Distribution for %s",aa[i]);
462 break;
463 default:
464 sprintf(hisfile,"histo-chi%d%s.xvg",Dih-NONCHI+1,aa[i]);
465 sprintf(title,"\\8c\\4\\s%d\\N Distribution for %s",
466 Dih-NONCHI+1,aa[i]);
467 }
468 fp=xvgropen(hisfile,title,"Degrees","");
469 fprintf(fp,"@ with g0\n");
470 xvgr_world(fp,-180,0,180,0.1);
471 fprintf(fp,"@ xaxis tick on\n");
472 fprintf(fp,"@ xaxis tick major 90\n");
473 fprintf(fp,"@ xaxis tick minor 30\n");
474 fprintf(fp,"@ xaxis ticklabel prec 0\n");
475 fprintf(fp,"@ yaxis tick off\n");
476 fprintf(fp,"@ yaxis ticklabel off\n");
477 fprintf(fp,"@ type xy\n");
478 for(j=0; (j<NHISTO); j++)
479 fprintf(fp,"%5d %10g\n",j-180,normhisto[j]);
480 fprintf(fp,"&\n");
481 ffclose(fp);
482 }
483 }
484 }
485 sfree(normhisto);
486 }
487
488 static FILE *rama_file(char *fn,char *title,char *xaxis,char *yaxis)
489 {
490 FILE *fp;
491
492 fp = xvgropen(fn,title,xaxis,yaxis);
493 fprintf(fp,"@ with g0\n");
494 xvgr_world(fp,-180,-180,180,180);
495 fprintf(fp,"@ xaxis tick on\n");
496 fprintf(fp,"@ xaxis tick major 90\n");
497 fprintf(fp,"@ xaxis tick minor 30\n");
498 fprintf(fp,"@ xaxis ticklabel prec 0\n");
499 fprintf(fp,"@ yaxis tick on\n");
500 fprintf(fp,"@ yaxis tick major 90\n");
501 fprintf(fp,"@ yaxis tick minor 30\n");
502 fprintf(fp,"@ yaxis ticklabel prec 0\n");
503 fprintf(fp,"@ s0 type xy\n");
504 fprintf(fp,"@ s0 symbol 2\n");
505 fprintf(fp,"@ s0 symbol size 0.410000\n");
506 fprintf(fp,"@ s0 symbol fill 1\n");
507 fprintf(fp,"@ s0 symbol color 1\n");
508 fprintf(fp,"@ s0 symbol linewidth 1\n");
509 fprintf(fp,"@ s0 symbol linestyle 1\n");
510 fprintf(fp,"@ s0 symbol center false\n");
511 fprintf(fp,"@ s0 symbol char 0\n");
512 fprintf(fp,"@ s0 skip 0\n");
513 fprintf(fp,"@ s0 linestyle 0\n");
514 fprintf(fp,"@ s0 linewidth 1\n");
515 fprintf(fp,"@ type xy\n");
516
517 return fp;
518 }
519
520 static void do_rama(int nf,int nlist,t_dlist dlist[],real **dih,
521 bool bViol,bool bRamOmega)
522 {
523 FILE *fp,*gp=NULL;
524 bool bOm;
525 char fn[256];
526 int i,j,k,Xi1,Xi2,Phi,Psi,Om=0,nlevels;
527 #define NMAT 120
528 real **mat=NULL,phi,psi,omega,axis[NMAT],lo,hi;
529 t_rgb rlo = { 1.0, 0.0, 0.0 };
530 t_rgb rmid= { 1.0, 1.0, 1.0 };
531 t_rgb rhi = { 0.0, 0.0, 1.0 };
532
533 for(i=0; (i<nlist); i++) {
534 if ((has_dihedral(edPhi,&(dlist[i]))) &&
535 (has_dihedral(edPsi,&(dlist[i])))) {
536 sprintf(fn,"ramaPhiPsi%s.xvg",dlist[i].name);
537 fp = rama_file(fn,"Ramachandran Plot",
538 "\\8f\\4 (deg)","\\8y\\4 (deg)");
539 bOm = bRamOmega && has_dihedral(edOmega,&(dlist[i]));
540 if (bOm) {
541 Om = dlist[i].j0[edOmega];
542 snew(mat,NMAT);
543 for(j=0; (j<NMAT); j++) {
544 snew(mat[j],NMAT);
545 axis[j] = -180+(360*j)/NMAT;
546 }
547 }
548 if (bViol) {
549 sprintf(fn,"violPhiPsi%s.xvg",dlist[i].name);
550 gp = ffopen(fn,"w");
551 }
552 Phi = dlist[i].j0[edPhi];
553 Psi = dlist[i].j0[edPsi];
554 for(j=0; (j<nf); j++) {
555 phi = RAD2DEG*dih[Phi][j];
556 psi = RAD2DEG*dih[Psi][j];
557 fprintf(fp,"%10g %10g\n",phi,psi);
558 if (bViol)
559 fprintf(gp,"%d\n",!bAllowed(dih[Phi][j],RAD2DEG*dih[Psi][j]));
560 if (bOm) {
561 omega = RAD2DEG*dih[Om][j];
562 mat[(int)((phi*NMAT)/360)+NMAT/2][(int)((psi*NMAT)/360)+NMAT/2]
563 += omega;
564 }
565 }
566 if (bViol)
567 fclose(gp);
568 fclose(fp);
569 if (bOm) {
570 sprintf(fn,"ramomega%s.xpm",dlist[i].name);
571 fp = ffopen(fn,"w");
572 lo = hi = 0;
573 for(j=0; (j<NMAT); j++)
574 for(k=0; (k<NMAT); k++) {
575 mat[j][k] /= nf;
576 lo=min(mat[j][k],lo);
577 hi=max(mat[j][k],hi);
578 }
579 /* Symmetrise */
580 if (fabs(lo) > fabs(hi))
581 hi = -lo;
582 else
583 lo = -hi;
584 /* Add 180 */
585 for(j=0; (j<NMAT); j++)
586 for(k=0; (k<NMAT); k++)
587 mat[j][k] += 180;
588 lo += 180;
589 hi += 180;
590 nlevels = 20;
591 write_xpm3(fp,"Omega/Ramachandran Plot","Deg","Phi","Psi",
592 NMAT,NMAT,axis,axis,mat,lo,180.0,hi,rlo,rmid,rhi,&nlevels);
593 fclose(fp);
594 for(j=0; (j<NMAT); j++)
595 sfree(mat[j]);
596 sfree(mat);
597 }
598 }
599 if ((has_dihedral(edChi1,&(dlist[i]))) &&
600 (has_dihedral(edChi2,&(dlist[i])))) {
601 sprintf(fn,"ramaX1X2%s.xvg",dlist[i].name);
602 fp = rama_file(fn,"\\8c\\4\\s1\\N-\\8c\\4\\s2\\N Ramachandran Plot",
603 "\\8c\\4\\s1\\N (deg)","\\8c\\4\\s2\\N (deg)");
604 Xi1 = dlist[i].j0[edChi1];
605 Xi2 = dlist[i].j0[edChi2];
606 for(j=0; (j<nf); j++)
607 fprintf(fp,"%10g %10g\n",RAD2DEG*dih[Xi1][j],RAD2DEG*dih[Xi2][j]);
608 fclose(fp);
609 }
610 else
611 fprintf(stderr,"No chi1 & chi2 angle for %s\n",dlist[i].name);
612 }
613 }
614
615 static void order_params(FILE *log,
616 char *fn,int maxchi,int nlist,t_dlist dlist[],
617 bool bPDB,char *pdbfn,real bfac_init,
618 t_topology *top,rvec x[],
619 bool bPhi,bool bPsi,bool bChi)
620 {
621 FILE *fp;
622 int nh[edMax];
623 int i,Dih,Xi;
624 real S2Max, S2Min;
625
626 /* Print order parameters */
627 fp=xvgropen(fn,"Dihedral Order Parameters","Residue","S2");
628
629 for (Dih=0; (Dih<edMax); Dih++)
630 nh[Dih]=0;
631
632 fprintf(fp,"%5s ","#Res.");
633 fprintf(fp,"%10s %10s ","S2Min","S2Max");
634 fprintf(fp,"%10s %10s ","Phi","Psi");
635 for (Xi=1; (Xi<=maxchi); Xi++)
636 fprintf(fp,"%8s%2d ","Xi",Xi);
637 fprintf(fp,"%12s\n","Res. Name");
638
639 for(i=0; (i<nlist); i++) {
640 S2Max=-10;
641 S2Min=10;
642 for (Dih=0; (Dih<NONCHI+maxchi); Dih++) {
643 if (dlist[i].S2[Dih]!=0) {
644 if (dlist[i].S2[Dih] > S2Max)
645 S2Max=dlist[i].S2[Dih];
646 if (dlist[i].S2[Dih] < S2Min)
647 S2Min=dlist[i].S2[Dih];
648 }
649 if (dlist[i].S2[Dih] > 0.8)
650 nh[Dih]++;
651 }
652 fprintf(fp,"%5d ",dlist[i].resnr);
653 fprintf(fp,"%10.3f %10.3f ",S2Min,S2Max);
654 for (Dih=0; (Dih<NONCHI+maxchi); Dih++)
655 fprintf(fp,"%10.3f ",dlist[i].S2[Dih]);
656 fprintf(fp,"%12s\n",dlist[i].name);
657 }
658 ffclose(fp);
659
660 if (bPDB) {
661 int Xi;
662 real x0,y0,z0;
663
664 for(i=0; (i<top->atoms.nr); i++)
665 top->atoms.pdbinfo[i].bfac=bfac_init;
666
667 for(i=0; (i<nlist); i++) {
668 top->atoms.pdbinfo[dlist[i].atm.N].bfac=-dlist[i].S2[0];/* Phi */
669 top->atoms.pdbinfo[dlist[i].atm.H].bfac=-dlist[i].S2[0];/* Phi */
670 top->atoms.pdbinfo[dlist[i].atm.C].bfac=-dlist[i].S2[1];/* Psi */
671 top->atoms.pdbinfo[dlist[i].atm.O].bfac=-dlist[i].S2[1];/* Psi */
672 for (Xi=0; (Xi<maxchi); Xi++) { /* Chi's */
673 if (dlist[i].atm.Cn[Xi+3]!=-1) {
674 top->atoms.pdbinfo[dlist[i].atm.Cn[Xi+1]].bfac=-dlist[i].S2[NONCHI+Xi];
675 }
676 }
677 }
678
679 fp=ffopen(fn,"w");
680 fprintf(fp,"REMARK generated by g_chi\n");
681 fprintf(fp,"REMARK B-factor field contains negative of\n");
682 fprintf(fp,"REMARK dihedral order parameters\n");
683 write_pdbfile(fp,NULL,&(top->atoms),x,NULL,0,TRUE);
684 x0=y0=z0=1000.0;
685 for (i=0; (i<top->atoms.nr); i++) {
686 if (x[i][XX]<x0) x0=x[i][XX];
687 if (x[i][YY]<y0) y0=x[i][YY];
688 if (x[i][ZZ]<y0) z0=x[i][ZZ];
689 }
690 x0*=10.0;/* nm -> angstrom */
691 y0*=10.0;/* nm -> angstrom */
692 z0*=10.0;/* nm -> angstrom */
693 for (i=0; (i<10); i++)
694 fprintf(fp,"%6s%5d %-4.4s%3.3s %c%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
695 "ATOM ",10000+i,"CA","LEG",' ',1000,x0,y0,z0+(1.2*i),
696 0.0,-0.1*i);
697 ffclose(fp);
698 }
699
700 fprintf(log,"Dihedrals with S2 > 0.8\n");
701 fprintf(log,"Dihedral: ");
702 if (bPhi)
703 fprintf(log," Phi ");
704 if (bPsi)
705 fprintf(log," Psi ");
706 if (bChi)
707 for(Xi=0; (Xi<maxchi); Xi++)
708 fprintf(log,"Chi%d ",Xi+1);
709 fprintf(log,"\nNumber: ");
710 if (bPhi)
711 fprintf(log,"%4d ",nh[0]);
712 if (bPsi)
713 fprintf(log,"%4d ",nh[1]);
714 if (bChi)
715 for(Xi=0; (Xi<maxchi); Xi++)
716 fprintf(log,"%4d ",nh[NONCHI+Xi]);
717 fprintf(log,"\n");
718 }
719
720 int main(int argc,char *argv[])
721 {
722 static char *desc[] = {
723 "g_chi computes phi, psi, omega and chi dihedrals for all your ",
724 "amino acid backbone and sidechains.",
725 "It can compute dihedral angle as a function of time, and as",
726 "histogram distributions.",
727 "Output is in form of xvgr files, as well as a LaTeX table of the",
728 "number of transitions per nanosecond.[PAR]",
729 "Order parameters S2 for each of the dihedrals are calculated and",
730 "output as xvgr file and optionally as a pdb file with the S2",
731 "values as B-factor.[PAR]",
732 "If option [TT]-c[tt] is given, the program will",
733 "calculate dihedral autocorrelation functions. The function used",
734 "is C(t) = < cos(chi(tau)) cos(chi(tau+t)) >. The use of cosines",
735 "rather than angles themselves, resolves the problem of periodicity.",
736 "(Van der Spoel & Berendsen (1997), [BB]Biophys. J. 72[bb], 2032-2041).[PAR]",
737 "The option [TT]-r[tt] generates a contour plot of the average omega angle",
738 "as a function of the phi and psi angles, that is, in a Ramachandran plot",
739 "the average omega angle is plotted using color coding."
740 };
741
742 static char *bugs[] = {
743 "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."
744 };
745 static int r0=1,ndeg=1,maxchi=2;
746 static bool bAll=FALSE;
747 static bool bPhi=FALSE,bPsi=FALSE,bOmega=FALSE;
748 static real bfac_init=-1.0;
749 static char *maxchistr[] = { NULL, "0", "1", "2", "3", "4", "5", "6", NULL };
750 static bool bRama=FALSE,bShift=FALSE,bViol=FALSE,bRamOmega=FALSE;
751 t_pargs pa[] = {
752 { "-r0", FALSE, etINT, {&r0},
753 "starting residue" },
754 { "-phi", FALSE, etBOOL, {&bPhi},
755 "Output for Phi dihedral angles" },
756 { "-psi", FALSE, etBOOL, {&bPsi},
757 "Output for Psi dihedral angles" },
758 { "-omega",FALSE, etBOOL, {&bOmega},
759 "Output for Omega dihedrals (peptide bonds)" },
760 { "-rama", FALSE, etBOOL, {&bRama},
761 "Generate Phi/Psi and Chi1/Chi2 ramachandran plots" },
762 { "-viol", FALSE, etBOOL, {&bViol},
763 "Write a file that gives 0 or 1 for violated Ramachandran angles" },
764 { "-all", FALSE, etBOOL, {&bAll},
765 "Output separate files for every dihedral." },
766 { "-shift", FALSE, etBOOL, {&bShift},
767 "Compute chemical shifts from Phi/Psi angles" },
768 { "-run", FALSE, etINT, {&ndeg},
769 "perform running average over ndeg degrees for histograms" },
770 { "-maxchi", FALSE, etENUM, {maxchistr},
771 "calculate first ndih Chi dihedrals" },
772 { "-ramomega",FALSE,etBOOL, {&bRamOmega},
773 "compute average omega as a function of phi/psi and plot it in an xpm plot" },
774 { "-bfact", FALSE, etREAL, {&bfac_init},
775 "bfactor value for pdb file for atoms with no calculated dihedral order parameter"}
776 };
777
778 FILE *log;
779 int natoms,nlist,naa,idum;
780 t_topology top;
781 rvec *x;
782 t_dlist *dlist;
783 char **aa;
784 bool bChi,bCorr;
785 real dt=0;
786
787 atom_id isize,*index;
788 int ndih,nf;
789 real **dih,*trans_frac,*aver_angle,*time;
790
791 t_filenm fnm[] = {
792 { efTPX, NULL, NULL, ffREAD },
793 { efTRX, "-f", NULL, ffREAD },
794 { efXVG, "-o", "order", ffWRITE },
795 { efPDB, "-p", "order", ffOPTWR },
796 { efXVG, "-jc", "Jcoupling", ffWRITE },
797 { efXVG, "-c", "dihcorr",ffOPTWR },
798 { efLOG, "-g", "chi", ffWRITE }
799 };
800 #define NFILE asize(fnm)
801 int npargs;
802 t_pargs *ppa;
803
804 CopyRight(stderr,argv[0]);
805 npargs = asize(pa);
806 ppa = add_acf_pargs(&npargs,pa);
807 parse_common_args(&argc,argv,PCA_CAN_VIEW | PCA_CAN_TIME,TRUE,
808 NFILE,fnm,npargs,ppa,asize(desc),desc,asize(bugs),bugs);
809
810 /* Handle result from enumerated type */
811 sscanf(maxchistr[0],"%d",&maxchi);
812 bChi = (maxchi > 0);
813
814 log=ffopen(ftp2fn(efLOG,NFILE,fnm),"w");
815
816 if (bRamOmega) {
817 bOmega = TRUE;
818 bPhi = TRUE;
819 bPsi = TRUE;
820 }
821
822 bCorr=(opt2bSet("-c",NFILE,fnm));
823 if (bCorr)
824 fprintf(stderr,"Will calculate autocorrelation\n");
825
826 if (maxchi > MAXCHI) {
827 fprintf(stderr,
828 "Will only calculate first %d Chi dihedrals in stead of %d.\n",
829 MAXCHI, maxchi);
830 maxchi=MAXCHI;
831 }
832
833 /* Find the chi angles using topology and a list of amino acids */
834 {
835 t_tpxheader sh;
836 int dint;/* dummy */
837 real dreal;/* dummy */
838 char *fn;
839
840 fn=ftp2fn(efTPX,NFILE,fnm);
841 read_tpxheader(fn, &sh);
842 natoms = sh.natoms;
843 snew(x,natoms);
844 read_tpx(fn,&dint,&dreal,&dreal,NULL,NULL,
845 &natoms,x,NULL,NULL,&top);
846 snew(top.atoms.pdbinfo,top.atoms.nr);
847 fprintf(log,"topology name = %s\n",*top.name);
848 }
849
850 naa=get_strings("aminoacids.dat",&aa);
851 dlist=mk_dlist(log,&(top.atoms),&nlist,bPhi,bPsi,bChi,maxchi,r0,naa,aa);
852 fprintf(stderr,"%d residues with dihedrals found\n", nlist);
853
854 if (nlist == 0)
855 fatal_error(0,"No dihedrals in your topology!\n");
856
857 /* Make a linear index for reading all */
858 index=make_chi_ind(nlist,dlist,&ndih);
859 isize=4*ndih;
860 fprintf(stderr,"%d dihedrals found\n", ndih);
861
862 snew(dih,ndih);
863
864 /* COMPUTE ALL DIHEDRALS! */
865 read_ang_dih(opt2fn("-f",NFILE,fnm),ftp2fn(efTPX,NFILE,fnm),
866 FALSE,TRUE,FALSE,1,&idum,
867 &nf,&time,isize,index,&trans_frac,&aver_angle,dih);
868
869 if (bCorr) {
870 if (nf < 2)
871 fatal_error(0,"Need at least 2 frames for correlation");
872
873 dt=(time[nf-1]-time[0])/(nf-1);
874 }
875
876 /* put angles in -M_PI to M_PI ! and correct phase factor for phi and psi */
877 reset_em_all(nlist,dlist,nf,dih,maxchi,bPhi,bPsi,bChi);
878
879 if (bAll)
880 dump_em_all(nlist,dlist,nf,time,dih,maxchi,bPhi,bPsi,bChi,bOmega);
881
882 /* Histogramming & J coupling constants */
883 histogramming(log,naa,aa,nf,maxchi,dih,nlist,dlist,bPhi,bPsi,bOmega,bChi);
884
885 /* Order parameters */
886 order_params(log,opt2fn("-o",NFILE,fnm),maxchi,nlist,dlist,
887 opt2bSet("-p",NFILE,fnm),ftp2fn(efPDB,NFILE,fnm),bfac_init,
888 &top,x,bPhi,bPsi,bChi);
889
890 /* Print ramachandran maps! */
891 if (bRama)
892 do_rama(nf,nlist,dlist,dih,bViol,bRamOmega);
893
894 if (bShift)
895 do_pp2shifts(log,nf,nlist,dlist,dih);
896
897 pr_dlist(log,nlist,dlist,time[nf-1]-time[0]);
898 ffclose(log);
899
900 /* Correlation comes last because it fucks up the angles */
901 if (bCorr)
902 do_dihcorr(opt2fn("-c",NFILE,fnm),nf,ndih,dih,dt,nlist,dlist,time,
903 maxchi,bPhi,bPsi,bChi,bOmega);
904
905
906 xvgr_file(opt2fn("-o",NFILE,fnm),"-nxy");
907 xvgr_file(opt2fn("-jc",NFILE,fnm),"-nxy");
908 if (bCorr)
909 xvgr_file(opt2fn("-c",NFILE,fnm),"-nxy");
910
911 thanx(stdout);
912
913 return 0;
914 }
got merge a long time ago