2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
44 #include "gromacs/topology/residuetypes.h"
45 #include "gromacs/utility/fatalerror.h"
46 #include "gromacs/utility/smalloc.h"
48 t_dlist
*mk_dlist(FILE *log
,
49 t_atoms
*atoms
, int *nlist
,
50 gmx_bool bPhi
, gmx_bool bPsi
, gmx_bool bChi
, gmx_bool bHChi
,
51 int maxchi
, int r0
, gmx_residuetype_t
*rt
)
53 int ires
, i
, j
, k
, ii
;
55 int nl
= 0, nc
[edMax
];
59 snew(dl
, atoms
->nres
+1);
60 prev
.C
= prev
.Cn
[1] = -1; /* Keep the compiler quiet */
61 for (i
= 0; (i
< edMax
); i
++)
69 ires
= atoms
->atom
[i
].resind
;
71 /* Initiate all atom numbers to -1 */
72 atm
.minC
= atm
.H
= atm
.N
= atm
.C
= atm
.O
= atm
.minCalpha
= -1;
73 for (j
= 0; (j
< MAXCHI
+3); j
++)
78 /* Look for atoms in this residue */
79 /* maybe should allow for chis to hydrogens? */
80 while ((i
< atoms
->nr
) && (atoms
->atom
[i
].resind
== ires
))
82 if ((strcmp(*(atoms
->atomname
[i
]), "H") == 0) ||
83 (strcmp(*(atoms
->atomname
[i
]), "H1") == 0) ||
84 (strcmp(*(atoms
->atomname
[i
]), "HN") == 0) )
88 else if (strcmp(*(atoms
->atomname
[i
]), "N") == 0)
92 else if (strcmp(*(atoms
->atomname
[i
]), "C") == 0)
96 else if ((strcmp(*(atoms
->atomname
[i
]), "O") == 0) ||
97 (strcmp(*(atoms
->atomname
[i
]), "O1") == 0) ||
98 (strcmp(*(atoms
->atomname
[i
]), "OC1") == 0) ||
99 (strcmp(*(atoms
->atomname
[i
]), "OT1") == 0))
103 else if (strcmp(*(atoms
->atomname
[i
]), "CA") == 0)
107 else if (strcmp(*(atoms
->atomname
[i
]), "CB") == 0)
111 else if ((strcmp(*(atoms
->atomname
[i
]), "CG") == 0) ||
112 (strcmp(*(atoms
->atomname
[i
]), "CG1") == 0) ||
113 (strcmp(*(atoms
->atomname
[i
]), "OG") == 0) ||
114 (strcmp(*(atoms
->atomname
[i
]), "OG1") == 0) ||
115 (strcmp(*(atoms
->atomname
[i
]), "SG") == 0))
119 else if ((strcmp(*(atoms
->atomname
[i
]), "CD") == 0) ||
120 (strcmp(*(atoms
->atomname
[i
]), "CD1") == 0) ||
121 (strcmp(*(atoms
->atomname
[i
]), "SD") == 0) ||
122 (strcmp(*(atoms
->atomname
[i
]), "OD1") == 0) ||
123 (strcmp(*(atoms
->atomname
[i
]), "ND1") == 0))
127 /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
128 else if (bHChi
&& ((strcmp(*(atoms
->atomname
[i
]), "HG") == 0) ||
129 (strcmp(*(atoms
->atomname
[i
]), "HG1") == 0)) )
133 else if ((strcmp(*(atoms
->atomname
[i
]), "CE") == 0) ||
134 (strcmp(*(atoms
->atomname
[i
]), "CE1") == 0) ||
135 (strcmp(*(atoms
->atomname
[i
]), "OE1") == 0) ||
136 (strcmp(*(atoms
->atomname
[i
]), "NE") == 0))
140 else if ((strcmp(*(atoms
->atomname
[i
]), "CZ") == 0) ||
141 (strcmp(*(atoms
->atomname
[i
]), "NZ") == 0))
145 /* HChi flag here too */
146 else if (bHChi
&& (strcmp(*(atoms
->atomname
[i
]), "NH1") == 0))
153 thisres
= *(atoms
->resinfo
[ires
].name
);
155 /* added by grs - special case for aromatics, whose chis above 2 are
156 not real and produce rubbish output - so set back to -1 */
157 if (strcmp(thisres
, "PHE") == 0 ||
158 strcmp(thisres
, "TYR") == 0 ||
159 strcmp(thisres
, "PTR") == 0 ||
160 strcmp(thisres
, "TRP") == 0 ||
161 strcmp(thisres
, "HIS") == 0 ||
162 strcmp(thisres
, "HISA") == 0 ||
163 strcmp(thisres
, "HISB") == 0)
165 for (ii
= 5; ii
<= 7; ii
++)
170 /* end fixing aromatics */
172 /* Special case for Pro, has no H */
173 if (strcmp(thisres
, "PRO") == 0)
177 /* Carbon from previous residue */
182 /* Alpha-carbon from previous residue */
183 if (prev
.Cn
[1] != -1)
185 atm
.minCalpha
= prev
.Cn
[1];
189 /* Check how many dihedrals we have */
190 if ((atm
.N
!= -1) && (atm
.Cn
[1] != -1) && (atm
.C
!= -1) &&
191 (atm
.O
!= -1) && ((atm
.H
!= -1) || (atm
.minC
!= -1)))
193 dl
[nl
].resnr
= ires
+1;
195 dl
[nl
].atm
.Cn
[0] = atm
.N
;
196 if ((atm
.Cn
[3] != -1) && (atm
.Cn
[2] != -1) && (atm
.Cn
[1] != -1))
220 if ((atm
.minC
!= -1) && (atm
.minCalpha
!= -1))
224 dl
[nl
].index
= gmx_residuetype_get_index(rt
, thisres
);
226 sprintf(dl
[nl
].name
, "%s%d", thisres
, ires
+r0
);
231 fprintf(debug
, "Could not find N atom but could find other atoms"
232 " in residue %s%d\n", thisres
, ires
+r0
);
235 fprintf(stderr
, "\n");
237 fprintf(log
, "There are %d residues with dihedrals\n", nl
);
249 for (i
= 0; (i
< maxchi
); i
++)
254 fprintf(log
, "There are %d dihedrals\n", j
);
255 fprintf(log
, "Dihedral: ");
258 fprintf(log
, " Phi ");
262 fprintf(log
, " Psi ");
266 for (i
= 0; (i
< maxchi
); i
++)
268 fprintf(log
, "Chi%d ", i
+1);
271 fprintf(log
, "\nNumber: ");
274 fprintf(log
, "%4d ", nl
);
278 fprintf(log
, "%4d ", nl
);
282 for (i
= 0; (i
< maxchi
); i
++)
284 fprintf(log
, "%4d ", nc
[i
]);
294 gmx_bool
has_dihedral(int Dih
, t_dlist
*dl
)
302 b
= ((dl
->atm
.H
!= -1) && (dl
->atm
.N
!= -1) && (dl
->atm
.Cn
[1] != -1) && (dl
->atm
.C
!= -1));
305 b
= ((dl
->atm
.N
!= -1) && (dl
->atm
.Cn
[1] != -1) && (dl
->atm
.C
!= -1) && (dl
->atm
.O
!= -1));
308 b
= ((dl
->atm
.minCalpha
!= -1) && (dl
->atm
.minC
!= -1) && (dl
->atm
.N
!= -1) && (dl
->atm
.Cn
[1] != -1));
317 b
= ((dl
->atm
.Cn
[ddd
] != -1) && (dl
->atm
.Cn
[ddd
+1] != -1) &&
318 (dl
->atm
.Cn
[ddd
+2] != -1) && (dl
->atm
.Cn
[ddd
+3] != -1));
321 pr_dlist(stdout
, 1, dl
, 1, 0, TRUE
, TRUE
, TRUE
, TRUE
, MAXCHI
);
322 gmx_fatal(FARGS
, "Non existant dihedral %d in file %s, line %d",
323 Dih
, __FILE__
, __LINE__
);
328 static void pr_one_ro(FILE *fp
, t_dlist
*dl
, int nDih
, real gmx_unused dt
)
331 for (k
= 0; k
< NROT
; k
++)
333 fprintf(fp
, " %6.2f", dl
->rot_occ
[nDih
][k
]);
338 static void pr_ntr_s2(FILE *fp
, t_dlist
*dl
, int nDih
, real dt
)
340 fprintf(fp
, " %6.2f %6.2f\n", (dt
== 0) ? 0 : dl
->ntr
[nDih
]/dt
, dl
->S2
[nDih
]);
343 void pr_dlist(FILE *fp
, int nl
, t_dlist dl
[], real dt
, int printtype
,
344 gmx_bool bPhi
, gmx_bool bPsi
, gmx_bool bChi
, gmx_bool bOmega
, int maxchi
)
348 void (*pr_props
)(FILE *, t_dlist
*, int, real
);
350 /* Analysis of dihedral transitions etc */
352 if (printtype
== edPrintST
)
354 pr_props
= pr_ntr_s2
;
355 fprintf(stderr
, "Now printing out transitions and OPs...\n");
359 pr_props
= pr_one_ro
;
360 fprintf(stderr
, "Now printing out rotamer occupancies...\n");
361 fprintf(fp
, "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
364 /* change atom numbers from 0 based to 1 based */
365 for (i
= 0; (i
< nl
); i
++)
367 fprintf(fp
, "Residue %s\n", dl
[i
].name
);
368 if (printtype
== edPrintST
)
370 fprintf(fp
, " Angle [ AI, AJ, AK, AL] #tr/ns S^2D \n"
371 "--------------------------------------------\n");
375 fprintf(fp
, " Angle [ AI, AJ, AK, AL] rotamers 0 g(-) t g(+)\n"
376 "--------------------------------------------\n");
380 fprintf(fp
, " Phi [%5d,%5d,%5d,%5d]",
381 (dl
[i
].atm
.H
== -1) ? 1+dl
[i
].atm
.minC
: 1+dl
[i
].atm
.H
,
382 1+dl
[i
].atm
.N
, 1+dl
[i
].atm
.Cn
[1], 1+dl
[i
].atm
.C
);
383 pr_props(fp
, &dl
[i
], edPhi
, dt
);
387 fprintf(fp
, " Psi [%5d,%5d,%5d,%5d]", 1+dl
[i
].atm
.N
, 1+dl
[i
].atm
.Cn
[1],
388 1+dl
[i
].atm
.C
, 1+dl
[i
].atm
.O
);
389 pr_props(fp
, &dl
[i
], edPsi
, dt
);
391 if (bOmega
&& has_dihedral(edOmega
, &(dl
[i
])))
393 fprintf(fp
, " Omega [%5d,%5d,%5d,%5d]", 1+dl
[i
].atm
.minCalpha
, 1+dl
[i
].atm
.minC
,
394 1+dl
[i
].atm
.N
, 1+dl
[i
].atm
.Cn
[1]);
395 pr_props(fp
, &dl
[i
], edOmega
, dt
);
397 for (Xi
= 0; Xi
< MAXCHI
; Xi
++)
399 if (bChi
&& (Xi
< maxchi
) && (dl
[i
].atm
.Cn
[Xi
+3] != -1) )
401 fprintf(fp
, " Chi%d[%5d,%5d,%5d,%5d]", Xi
+1, 1+dl
[i
].atm
.Cn
[Xi
],
402 1+dl
[i
].atm
.Cn
[Xi
+1], 1+dl
[i
].atm
.Cn
[Xi
+2],
403 1+dl
[i
].atm
.Cn
[Xi
+3]);
404 pr_props(fp
, &dl
[i
], Xi
+edChi1
, dt
); /* Xi+2 was wrong here */
413 int pr_trans(FILE *fp
, int nl
, t_dlist dl
[], real dt
, int Xi
)
415 /* never called at the moment */
420 fprintf(fp
, "\\begin{table}[h]\n");
421 fprintf(fp
, "\\caption{Number of dihedral transitions per nanosecond}\n");
422 fprintf(fp
, "\\begin{tabular}{|l|l|}\n");
423 fprintf(fp
, "\\hline\n");
424 fprintf(fp
, "Residue\t&$\\chi_%d$\t\\\\\n", Xi
+1);
425 for (i
= 0; (i
< nl
); i
++)
427 nn
= dl
[i
].ntr
[Xi
]/dt
;
431 fprintf(fp
, "%s\t&\\HL{%d}\t\\\\\n", dl
[i
].name
, nn
);
436 fprintf(fp
, "%s\t&\\%d\t\\\\\n", dl
[i
].name
, nn
);
439 fprintf(fp
, "\\hline\n");
440 fprintf(fp
, "\\end{tabular}\n");
441 fprintf(fp
, "\\end{table}\n\n");