3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
44 #include "gmx_fatal.h"
47 t_dlist
*mk_dlist(FILE *log
,
48 t_atoms
*atoms
, int *nlist
,
49 gmx_bool bPhi
, gmx_bool bPsi
, gmx_bool bChi
, gmx_bool bHChi
,
50 int maxchi
, int r0
, gmx_residuetype_t rt
)
58 snew(dl
,atoms
->nres
+1);
60 for(i
=0; (i
<edMax
); i
++)
65 ires
= atoms
->atom
[i
].resind
;
67 /* Initiate all atom numbers to -1 */
68 atm
.minC
=atm
.H
=atm
.N
=atm
.C
=atm
.O
=atm
.minO
=-1;
69 for(j
=0; (j
<MAXCHI
+3); j
++)
72 /* Look for atoms in this residue */
73 /* maybe should allow for chis to hydrogens? */
74 while ((i
<atoms
->nr
) && (atoms
->atom
[i
].resind
== ires
)) {
75 if ((strcmp(*(atoms
->atomname
[i
]),"H") == 0) ||
76 (strcmp(*(atoms
->atomname
[i
]),"H1") == 0) )
78 else if (strcmp(*(atoms
->atomname
[i
]),"N") == 0)
80 else if (strcmp(*(atoms
->atomname
[i
]),"C") == 0)
82 else if ((strcmp(*(atoms
->atomname
[i
]),"O") == 0) ||
83 (strcmp(*(atoms
->atomname
[i
]),"O1") == 0))
85 else if (strcmp(*(atoms
->atomname
[i
]),"CA") == 0)
87 else if (strcmp(*(atoms
->atomname
[i
]),"CB") == 0)
89 else if ((strcmp(*(atoms
->atomname
[i
]),"CG") == 0) ||
90 (strcmp(*(atoms
->atomname
[i
]),"CG1") == 0) ||
91 (strcmp(*(atoms
->atomname
[i
]),"OG") == 0) ||
92 (strcmp(*(atoms
->atomname
[i
]),"OG1") == 0) ||
93 (strcmp(*(atoms
->atomname
[i
]),"SG") == 0))
95 else if ((strcmp(*(atoms
->atomname
[i
]),"CD") == 0) ||
96 (strcmp(*(atoms
->atomname
[i
]),"CD1") == 0) ||
97 (strcmp(*(atoms
->atomname
[i
]),"SD") == 0) ||
98 (strcmp(*(atoms
->atomname
[i
]),"OD1") == 0) ||
99 (strcmp(*(atoms
->atomname
[i
]),"ND1") == 0))
101 /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
102 else if (bHChi
&& ((strcmp(*(atoms
->atomname
[i
]),"HG") == 0) ||
103 (strcmp(*(atoms
->atomname
[i
]),"HG1") == 0)) )
105 else if ((strcmp(*(atoms
->atomname
[i
]),"CE") == 0) ||
106 (strcmp(*(atoms
->atomname
[i
]),"CE1") == 0) ||
107 (strcmp(*(atoms
->atomname
[i
]),"OE1") == 0) ||
108 (strcmp(*(atoms
->atomname
[i
]),"NE") == 0))
110 else if ((strcmp(*(atoms
->atomname
[i
]),"CZ") == 0) ||
111 (strcmp(*(atoms
->atomname
[i
]),"NZ") == 0))
113 /* HChi flag here too */
114 else if (bHChi
&& (strcmp(*(atoms
->atomname
[i
]),"NH1") == 0))
119 thisres
= *(atoms
->resinfo
[ires
].name
);
121 /* added by grs - special case for aromatics, whose chis above 2 are
122 not real and produce rubbish output - so set back to -1 */
123 if (strcmp(thisres
,"PHE") == 0 ||
124 strcmp(thisres
,"TYR") == 0 ||
125 strcmp(thisres
,"PTR") == 0 ||
126 strcmp(thisres
,"TRP") == 0 ||
127 strcmp(thisres
,"HIS") == 0 ||
128 strcmp(thisres
,"HISA") == 0 ||
129 strcmp(thisres
,"HISB") == 0 ) {
130 for (ii
=5 ; ii
<=7 ; ii
++)
133 /* end fixing aromatics */
135 /* Special case for Pro, has no H */
136 if (strcmp(thisres
,"PRO") == 0)
138 /* Carbon from previous residue */
145 /* Check how many dihedrals we have */
146 if ((atm
.N
!= -1) && (atm
.Cn
[1] != -1) && (atm
.C
!= -1) &&
147 (atm
.O
!= -1) && ((atm
.H
!= -1) || (atm
.minC
!= -1))) {
148 dl
[nl
].resnr
= ires
+1;
150 dl
[nl
].atm
.Cn
[0] = atm
.N
;
151 if ((atm
.Cn
[3] != -1) && (atm
.Cn
[2] != -1) && (atm
.Cn
[1] != -1)) {
153 if (atm
.Cn
[4] != -1) {
155 if (atm
.Cn
[5] != -1) {
157 if (atm
.Cn
[6] != -1) {
159 if (atm
.Cn
[7] != -1) {
161 if (atm
.Cn
[8] != -1) {
169 if ((atm
.minC
!= -1) && (atm
.minO
!= -1))
171 dl
[nl
].index
=gmx_residuetype_get_index(rt
,thisres
);
173 sprintf(dl
[nl
].name
,"%s%d",thisres
,ires
+r0
);
177 fprintf(debug
,"Could not find N atom but could find other atoms"
178 " in residue %s%d\n",thisres
,ires
+r0
);
180 fprintf(stderr
,"\n");
182 fprintf(log
,"There are %d residues with dihedrals\n",nl
);
187 for(i
=0; (i
<maxchi
); i
++)
189 fprintf(log
,"There are %d dihedrals\n",j
);
190 fprintf(log
,"Dihedral: ");
192 fprintf(log
," Phi ");
194 fprintf(log
," Psi ");
196 for(i
=0; (i
<maxchi
); i
++)
197 fprintf(log
,"Chi%d ",i
+1);
198 fprintf(log
,"\nNumber: ");
200 fprintf(log
,"%4d ",nl
);
202 fprintf(log
,"%4d ",nl
);
204 for(i
=0; (i
<maxchi
); i
++)
205 fprintf(log
,"%4d ",nc
[i
]);
213 gmx_bool
has_dihedral(int Dih
,t_dlist
*dl
)
220 b
= ((dl
->atm
.H
!=-1) && (dl
->atm
.N
!=-1) && (dl
->atm
.Cn
[1]!=-1) && (dl
->atm
.C
!=-1));
223 b
= ((dl
->atm
.N
!=-1) && (dl
->atm
.Cn
[1]!=-1) && (dl
->atm
.C
!=-1) && (dl
->atm
.O
!=-1));
226 b
= ((dl
->atm
.minO
!=-1) && (dl
->atm
.minC
!=-1) && (dl
->atm
.N
!=-1) && (dl
->atm
.Cn
[1]!=-1));
235 b
= ((dl
->atm
.Cn
[ddd
]!=-1) && (dl
->atm
.Cn
[ddd
+1]!=-1)&&
236 (dl
->atm
.Cn
[ddd
+2]!=-1) && (dl
->atm
.Cn
[ddd
+3]!=-1));
239 pr_dlist(stdout
,1,dl
,1,0,TRUE
,TRUE
,TRUE
,TRUE
,MAXCHI
);
240 gmx_fatal(FARGS
,"Non existant dihedral %d in file %s, line %d",
241 Dih
,__FILE__
,__LINE__
);
246 static void pr_one_ro(FILE *fp
,t_dlist
*dl
,int nDih
,real dt
)
250 fprintf(fp
," %6.2f",dl
->rot_occ
[nDih
][k
]);
254 static void pr_ntr_s2(FILE *fp
,t_dlist
*dl
,int nDih
,real dt
)
256 fprintf(fp
," %6.2f %6.2f\n",(dt
== 0) ? 0 : dl
->ntr
[nDih
]/dt
,dl
->S2
[nDih
]);
259 void pr_dlist(FILE *fp
,int nl
,t_dlist dl
[],real dt
, int printtype
,
260 gmx_bool bPhi
, gmx_bool bPsi
,gmx_bool bChi
,gmx_bool bOmega
, int maxchi
)
264 void (*pr_props
)(FILE *, t_dlist
*, int, real
);
266 /* Analysis of dihedral transitions etc */
268 if (printtype
== edPrintST
){
270 fprintf(stderr
,"Now printing out transitions and OPs...\n");
273 fprintf(stderr
,"Now printing out rotamer occupancies...\n");
274 fprintf(fp
,"\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
277 /* change atom numbers from 0 based to 1 based */
278 for(i
=0; (i
<nl
); i
++) {
279 fprintf(fp
,"Residue %s\n",dl
[i
].name
);
280 if (printtype
== edPrintST
){
281 fprintf(fp
," Angle [ AI, AJ, AK, AL] #tr/ns S^2D \n"
282 "--------------------------------------------\n");
284 fprintf(fp
," Angle [ AI, AJ, AK, AL] rotamers 0 g(-) t g(+)\n"
285 "--------------------------------------------\n");
288 fprintf(fp
," Phi [%5d,%5d,%5d,%5d]",
289 (dl
[i
].atm
.H
== -1) ? 1+dl
[i
].atm
.minC
: 1+dl
[i
].atm
.H
,
290 1+dl
[i
].atm
.N
, 1+dl
[i
].atm
.Cn
[1], 1+dl
[i
].atm
.C
);
291 pr_props(fp
,&dl
[i
],edPhi
,dt
);
294 fprintf(fp
," Psi [%5d,%5d,%5d,%5d]",1+dl
[i
].atm
.N
, 1+dl
[i
].atm
.Cn
[1],
295 1+dl
[i
].atm
.C
, 1+dl
[i
].atm
.O
);
296 pr_props(fp
,&dl
[i
],edPsi
,dt
);
298 if (bOmega
&& has_dihedral(edOmega
,&(dl
[i
]))) {
299 fprintf(fp
," Omega [%5d,%5d,%5d,%5d]",1+dl
[i
].atm
.minO
, 1+dl
[i
].atm
.minC
,
300 1+dl
[i
].atm
.N
, 1+dl
[i
].atm
.Cn
[1]);
301 pr_props(fp
,&dl
[i
],edOmega
,dt
);
303 for(Xi
=0; Xi
<MAXCHI
; Xi
++)
304 if (bChi
&& (Xi
< maxchi
) && (dl
[i
].atm
.Cn
[Xi
+3] != -1) ) {
305 fprintf(fp
," Chi%d[%5d,%5d,%5d,%5d]",Xi
+1, 1+dl
[i
].atm
.Cn
[Xi
],
306 1+dl
[i
].atm
.Cn
[Xi
+1], 1+dl
[i
].atm
.Cn
[Xi
+2],
307 1+dl
[i
].atm
.Cn
[Xi
+3]);
308 pr_props(fp
,&dl
[i
],Xi
+edChi1
,dt
); /* Xi+2 was wrong here */
316 int pr_trans(FILE *fp
,int nl
,t_dlist dl
[],real dt
,int Xi
)
318 /* never called at the moment */
323 fprintf(fp
,"\\begin{table}[h]\n");
324 fprintf(fp
,"\\caption{Number of dihedral transitions per nanosecond}\n");
325 fprintf(fp
,"\\begin{tabular}{|l|l|}\n");
326 fprintf(fp
,"\\hline\n");
327 fprintf(fp
,"Residue\t&$\\chi_%d$\t\\\\\n",Xi
+1);
328 for(i
=0; (i
<nl
); i
++) {
332 fprintf(fp
,"%s\t&\\HL{%d}\t\\\\\n",dl
[i
].name
,nn
);
336 fprintf(fp
,"%s\t&\\%d\t\\\\\n",dl
[i
].name
,nn
);
338 fprintf(fp
,"\\hline\n");
339 fprintf(fp
,"\\end{tabular}\n");
340 fprintf(fp
,"\\end{table}\n\n");