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.
42 #include "gromacs/gmxana/gstat.h"
43 #include "gromacs/topology/residuetypes.h"
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/utility/smalloc.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)
52 int ires, i, j, k, ii;
54 int nl = 0, nc[edMax];
58 snew(dl, atoms->nres+1);
59 prev.C = prev.Cn[1] = -1; /* Keep the compiler quiet */
60 for (i = 0; (i < edMax); i++)
68 ires = atoms->atom[i].resind;
70 /* Initiate all atom numbers to -1 */
71 atm.minC = atm.H = atm.N = atm.C = atm.O = atm.minCalpha = -1;
72 for (j = 0; (j < MAXCHI+3); j++)
77 /* Look for atoms in this residue */
78 /* maybe should allow for chis to hydrogens? */
79 while ((i < atoms->nr) && (atoms->atom[i].resind == ires))
81 if ((strcmp(*(atoms->atomname[i]), "H") == 0) ||
82 (strcmp(*(atoms->atomname[i]), "H1") == 0) ||
83 (strcmp(*(atoms->atomname[i]), "HN") == 0) )
87 else if (strcmp(*(atoms->atomname[i]), "N") == 0)
91 else if (strcmp(*(atoms->atomname[i]), "C") == 0)
95 else if ((strcmp(*(atoms->atomname[i]), "O") == 0) ||
96 (strcmp(*(atoms->atomname[i]), "O1") == 0) ||
97 (strcmp(*(atoms->atomname[i]), "OC1") == 0) ||
98 (strcmp(*(atoms->atomname[i]), "OT1") == 0))
102 else if (strcmp(*(atoms->atomname[i]), "CA") == 0)
106 else if (strcmp(*(atoms->atomname[i]), "CB") == 0)
110 else if ((strcmp(*(atoms->atomname[i]), "CG") == 0) ||
111 (strcmp(*(atoms->atomname[i]), "CG1") == 0) ||
112 (strcmp(*(atoms->atomname[i]), "OG") == 0) ||
113 (strcmp(*(atoms->atomname[i]), "OG1") == 0) ||
114 (strcmp(*(atoms->atomname[i]), "SG") == 0))
118 else if ((strcmp(*(atoms->atomname[i]), "CD") == 0) ||
119 (strcmp(*(atoms->atomname[i]), "CD1") == 0) ||
120 (strcmp(*(atoms->atomname[i]), "SD") == 0) ||
121 (strcmp(*(atoms->atomname[i]), "OD1") == 0) ||
122 (strcmp(*(atoms->atomname[i]), "ND1") == 0))
126 /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
127 else if (bHChi && ((strcmp(*(atoms->atomname[i]), "HG") == 0) ||
128 (strcmp(*(atoms->atomname[i]), "HG1") == 0)) )
132 else if ((strcmp(*(atoms->atomname[i]), "CE") == 0) ||
133 (strcmp(*(atoms->atomname[i]), "CE1") == 0) ||
134 (strcmp(*(atoms->atomname[i]), "OE1") == 0) ||
135 (strcmp(*(atoms->atomname[i]), "NE") == 0))
139 else if ((strcmp(*(atoms->atomname[i]), "CZ") == 0) ||
140 (strcmp(*(atoms->atomname[i]), "NZ") == 0))
144 /* HChi flag here too */
145 else if (bHChi && (strcmp(*(atoms->atomname[i]), "NH1") == 0))
152 thisres = *(atoms->resinfo[ires].name);
154 /* added by grs - special case for aromatics, whose chis above 2 are
155 not real and produce rubbish output - so set back to -1 */
156 if (strcmp(thisres, "PHE") == 0 ||
157 strcmp(thisres, "TYR") == 0 ||
158 strcmp(thisres, "PTR") == 0 ||
159 strcmp(thisres, "TRP") == 0 ||
160 strcmp(thisres, "HIS") == 0 ||
161 strcmp(thisres, "HISA") == 0 ||
162 strcmp(thisres, "HISB") == 0)
164 for (ii = 5; ii <= 7; ii++)
169 /* end fixing aromatics */
171 /* Special case for Pro, has no H */
172 if (strcmp(thisres, "PRO") == 0)
176 /* Carbon from previous residue */
181 /* Alpha-carbon from previous residue */
182 if (prev.Cn[1] != -1)
184 atm.minCalpha = prev.Cn[1];
188 /* Check how many dihedrals we have */
189 if ((atm.N != -1) && (atm.Cn[1] != -1) && (atm.C != -1) &&
190 (atm.O != -1) && ((atm.H != -1) || (atm.minC != -1)))
192 dl[nl].resnr = ires+1;
194 dl[nl].atm.Cn[0] = atm.N;
195 if ((atm.Cn[3] != -1) && (atm.Cn[2] != -1) && (atm.Cn[1] != -1))
219 if ((atm.minC != -1) && (atm.minCalpha != -1))
223 dl[nl].index = gmx_residuetype_get_index(rt, thisres);
225 sprintf(dl[nl].name, "%s%d", thisres, ires+r0);
230 fprintf(debug, "Could not find N atom but could find other atoms"
231 " in residue %s%d\n", thisres, ires+r0);
234 fprintf(stderr, "\n");
236 fprintf(log, "There are %d residues with dihedrals\n", nl);
248 for (i = 0; (i < maxchi); i++)
253 fprintf(log, "There are %d dihedrals\n", j);
254 fprintf(log, "Dihedral: ");
257 fprintf(log, " Phi ");
261 fprintf(log, " Psi ");
265 for (i = 0; (i < maxchi); i++)
267 fprintf(log, "Chi%d ", i+1);
270 fprintf(log, "\nNumber: ");
273 fprintf(log, "%4d ", nl);
277 fprintf(log, "%4d ", nl);
281 for (i = 0; (i < maxchi); i++)
283 fprintf(log, "%4d ", nc[i]);
293 gmx_bool has_dihedral(int Dih, t_dlist *dl)
301 b = ((dl->atm.H != -1) && (dl->atm.N != -1) && (dl->atm.Cn[1] != -1) && (dl->atm.C != -1));
304 b = ((dl->atm.N != -1) && (dl->atm.Cn[1] != -1) && (dl->atm.C != -1) && (dl->atm.O != -1));
307 b = ((dl->atm.minCalpha != -1) && (dl->atm.minC != -1) && (dl->atm.N != -1) && (dl->atm.Cn[1] != -1));
316 b = ((dl->atm.Cn[ddd] != -1) && (dl->atm.Cn[ddd+1] != -1) &&
317 (dl->atm.Cn[ddd+2] != -1) && (dl->atm.Cn[ddd+3] != -1));
320 pr_dlist(stdout, 1, dl, 1, 0, TRUE, TRUE, TRUE, TRUE, MAXCHI);
321 gmx_fatal(FARGS, "Non existant dihedral %d in file %s, line %d",
322 Dih, __FILE__, __LINE__);
327 static void pr_one_ro(FILE *fp, t_dlist *dl, int nDih, real gmx_unused dt)
330 for (k = 0; k < NROT; k++)
332 fprintf(fp, " %6.2f", dl->rot_occ[nDih][k]);
337 static void pr_ntr_s2(FILE *fp, t_dlist *dl, int nDih, real dt)
339 fprintf(fp, " %6.2f %6.2f\n", (dt == 0) ? 0 : dl->ntr[nDih]/dt, dl->S2[nDih]);
342 void pr_dlist(FILE *fp, int nl, t_dlist dl[], real dt, int printtype,
343 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bOmega, int maxchi)
347 void (*pr_props)(FILE *, t_dlist *, int, real);
349 /* Analysis of dihedral transitions etc */
351 if (printtype == edPrintST)
353 pr_props = pr_ntr_s2;
354 fprintf(stderr, "Now printing out transitions and OPs...\n");
358 pr_props = pr_one_ro;
359 fprintf(stderr, "Now printing out rotamer occupancies...\n");
360 fprintf(fp, "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
363 /* change atom numbers from 0 based to 1 based */
364 for (i = 0; (i < nl); i++)
366 fprintf(fp, "Residue %s\n", dl[i].name);
367 if (printtype == edPrintST)
369 fprintf(fp, " Angle [ AI, AJ, AK, AL] #tr/ns S^2D \n"
370 "--------------------------------------------\n");
374 fprintf(fp, " Angle [ AI, AJ, AK, AL] rotamers 0 g(-) t g(+)\n"
375 "--------------------------------------------\n");
379 fprintf(fp, " Phi [%5d,%5d,%5d,%5d]",
380 (dl[i].atm.H == -1) ? 1+dl[i].atm.minC : 1+dl[i].atm.H,
381 1+dl[i].atm.N, 1+dl[i].atm.Cn[1], 1+dl[i].atm.C);
382 pr_props(fp, &dl[i], edPhi, dt);
386 fprintf(fp, " Psi [%5d,%5d,%5d,%5d]", 1+dl[i].atm.N, 1+dl[i].atm.Cn[1],
387 1+dl[i].atm.C, 1+dl[i].atm.O);
388 pr_props(fp, &dl[i], edPsi, dt);
390 if (bOmega && has_dihedral(edOmega, &(dl[i])))
392 fprintf(fp, " Omega [%5d,%5d,%5d,%5d]", 1+dl[i].atm.minCalpha, 1+dl[i].atm.minC,
393 1+dl[i].atm.N, 1+dl[i].atm.Cn[1]);
394 pr_props(fp, &dl[i], edOmega, dt);
396 for (Xi = 0; Xi < MAXCHI; Xi++)
398 if (bChi && (Xi < maxchi) && (dl[i].atm.Cn[Xi+3] != -1) )
400 fprintf(fp, " Chi%d[%5d,%5d,%5d,%5d]", Xi+1, 1+dl[i].atm.Cn[Xi],
401 1+dl[i].atm.Cn[Xi+1], 1+dl[i].atm.Cn[Xi+2],
402 1+dl[i].atm.Cn[Xi+3]);
403 pr_props(fp, &dl[i], Xi+edChi1, dt); /* Xi+2 was wrong here */
412 int pr_trans(FILE *fp, int nl, t_dlist dl[], real dt, int Xi)
414 /* never called at the moment */
419 fprintf(fp, "\\begin{table}[h]\n");
420 fprintf(fp, "\\caption{Number of dihedral transitions per nanosecond}\n");
421 fprintf(fp, "\\begin{tabular}{|l|l|}\n");
422 fprintf(fp, "\\hline\n");
423 fprintf(fp, "Residue\t&$\\chi_%d$\t\\\\\n", Xi+1);
424 for (i = 0; (i < nl); i++)
426 nn = dl[i].ntr[Xi]/dt;
430 fprintf(fp, "%s\t&\\HL{%d}\t\\\\\n", dl[i].name, nn);
435 fprintf(fp, "%s\t&\\%d\t\\\\\n", dl[i].name, nn);
438 fprintf(fp, "\\hline\n");
439 fprintf(fp, "\\end{tabular}\n");
440 fprintf(fp, "\\end{table}\n\n");