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,2015,2016,2018 by the GROMACS development team.
7 * Copyright (c) 2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
43 #include "gromacs/gmxana/gstat.h"
44 #include "gromacs/topology/residuetypes.h"
45 #include "gromacs/topology/topology.h"
46 #include "gromacs/utility/fatalerror.h"
47 #include "gromacs/utility/smalloc.h"
49 t_dlist* mk_dlist(FILE* log,
62 int nl = 0, nc[edMax];
66 snew(dl, atoms->nres + 1);
67 prev.C = prev.Cn[1] = -1; /* Keep the compiler quiet */
68 for (i = 0; (i < edMax); i++)
75 int ires = atoms->atom[i].resind;
77 /* Initiate all atom numbers to -1 */
78 atm.minC = atm.H = atm.N = atm.C = atm.O = atm.minCalpha = -1;
79 for (j = 0; (j < MAXCHI + 3); j++)
84 /* Look for atoms in this residue */
85 /* maybe should allow for chis to hydrogens? */
86 while ((i < atoms->nr) && (atoms->atom[i].resind == ires))
88 if ((std::strcmp(*(atoms->atomname[i]), "H") == 0)
89 || (std::strcmp(*(atoms->atomname[i]), "H1") == 0)
90 || (std::strcmp(*(atoms->atomname[i]), "HN") == 0))
94 else if (std::strcmp(*(atoms->atomname[i]), "N") == 0)
98 else if (std::strcmp(*(atoms->atomname[i]), "C") == 0)
102 else if ((std::strcmp(*(atoms->atomname[i]), "O") == 0)
103 || (std::strcmp(*(atoms->atomname[i]), "O1") == 0)
104 || (std::strcmp(*(atoms->atomname[i]), "OC1") == 0)
105 || (std::strcmp(*(atoms->atomname[i]), "OT1") == 0))
109 else if (std::strcmp(*(atoms->atomname[i]), "CA") == 0)
113 else if (std::strcmp(*(atoms->atomname[i]), "CB") == 0)
117 else if ((std::strcmp(*(atoms->atomname[i]), "CG") == 0)
118 || (std::strcmp(*(atoms->atomname[i]), "CG1") == 0)
119 || (std::strcmp(*(atoms->atomname[i]), "OG") == 0)
120 || (std::strcmp(*(atoms->atomname[i]), "OG1") == 0)
121 || (std::strcmp(*(atoms->atomname[i]), "SG") == 0))
125 else if ((std::strcmp(*(atoms->atomname[i]), "CD") == 0)
126 || (std::strcmp(*(atoms->atomname[i]), "CD1") == 0)
127 || (std::strcmp(*(atoms->atomname[i]), "SD") == 0)
128 || (std::strcmp(*(atoms->atomname[i]), "OD1") == 0)
129 || (std::strcmp(*(atoms->atomname[i]), "ND1") == 0))
130 { // NOLINT bugprone-branch-clone
133 /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
135 && ((std::strcmp(*(atoms->atomname[i]), "HG") == 0)
136 || (std::strcmp(*(atoms->atomname[i]), "HG1") == 0)))
140 else if ((std::strcmp(*(atoms->atomname[i]), "CE") == 0)
141 || (std::strcmp(*(atoms->atomname[i]), "CE1") == 0)
142 || (std::strcmp(*(atoms->atomname[i]), "OE1") == 0)
143 || (std::strcmp(*(atoms->atomname[i]), "NE") == 0))
147 else if ((std::strcmp(*(atoms->atomname[i]), "CZ") == 0)
148 || (std::strcmp(*(atoms->atomname[i]), "NZ") == 0))
152 /* HChi flag here too */
153 else if (bHChi && (std::strcmp(*(atoms->atomname[i]), "NH1") == 0))
160 thisres = *(atoms->resinfo[ires].name);
162 /* added by grs - special case for aromatics, whose chis above 2 are
163 not real and produce rubbish output - so set back to -1 */
164 if (std::strcmp(thisres, "PHE") == 0 || std::strcmp(thisres, "TYR") == 0
165 || std::strcmp(thisres, "PTR") == 0 || std::strcmp(thisres, "TRP") == 0
166 || std::strcmp(thisres, "HIS") == 0 || std::strcmp(thisres, "HISA") == 0
167 || std::strcmp(thisres, "HISB") == 0)
169 for (ii = 5; ii <= 7; ii++)
174 /* end fixing aromatics */
176 /* Special case for Pro, has no H */
177 if (std::strcmp(thisres, "PRO") == 0)
181 /* Carbon from previous residue */
186 /* Alpha-carbon from previous residue */
187 if (prev.Cn[1] != -1)
189 atm.minCalpha = prev.Cn[1];
193 /* Check how many dihedrals we have */
194 if ((atm.N != -1) && (atm.Cn[1] != -1) && (atm.C != -1) && (atm.O != -1)
195 && ((atm.H != -1) || (atm.minC != -1)))
197 dl[nl].resnr = ires + 1;
199 dl[nl].atm.Cn[0] = atm.N;
200 if ((atm.Cn[3] != -1) && (atm.Cn[2] != -1) && (atm.Cn[1] != -1))
224 if ((atm.minC != -1) && (atm.minCalpha != -1))
228 dl[nl].index = rt->indexFromResidueName(thisres);
230 /* Prevent seg fault from unknown residues. If one adds a custom residue to
231 * residuetypes.dat but somehow loses it, changes it, or does analysis on
232 * another machine, the residue type will be unknown. */
233 if (dl[nl].index == -1)
236 "Unknown residue %s when searching for residue type.\n"
237 "Maybe you need to add a custom residue in residuetypes.dat.",
241 sprintf(dl[nl].name, "%s%d", thisres, ires + r0);
247 "Could not find N atom but could find other atoms"
248 " in residue %s%d\n",
252 fprintf(stderr, "\n");
254 fprintf(log, "There are %d residues with dihedrals\n", nl);
266 for (i = 0; (i < maxchi); i++)
271 fprintf(log, "There are %d dihedrals\n", j);
272 fprintf(log, "Dihedral: ");
275 fprintf(log, " Phi ");
279 fprintf(log, " Psi ");
283 for (i = 0; (i < maxchi); i++)
285 fprintf(log, "Chi%d ", i + 1);
288 fprintf(log, "\nNumber: ");
291 fprintf(log, "%4d ", nl);
295 fprintf(log, "%4d ", nl);
299 for (i = 0; (i < maxchi); i++)
301 fprintf(log, "%4d ", nc[i]);
311 gmx_bool has_dihedral(int Dih, t_dlist* dl)
319 b = ((dl->atm.H != -1) && (dl->atm.N != -1) && (dl->atm.Cn[1] != -1) && (dl->atm.C != -1));
322 b = ((dl->atm.N != -1) && (dl->atm.Cn[1] != -1) && (dl->atm.C != -1) && (dl->atm.O != -1));
325 b = ((dl->atm.minCalpha != -1) && (dl->atm.minC != -1) && (dl->atm.N != -1)
326 && (dl->atm.Cn[1] != -1));
335 b = ((dl->atm.Cn[ddd] != -1) && (dl->atm.Cn[ddd + 1] != -1)
336 && (dl->atm.Cn[ddd + 2] != -1) && (dl->atm.Cn[ddd + 3] != -1));
339 pr_dlist(stdout, 1, dl, 1, 0, TRUE, TRUE, TRUE, TRUE, MAXCHI);
340 gmx_fatal(FARGS, "Non existent dihedral %d in file %s, line %d", Dih, __FILE__, __LINE__);
345 static void pr_one_ro(FILE* fp, t_dlist* dl, int nDih, real gmx_unused dt)
348 for (k = 0; k < NROT; k++)
350 fprintf(fp, " %6.2f", dl->rot_occ[nDih][k]);
355 static void pr_ntr_s2(FILE* fp, t_dlist* dl, int nDih, real dt)
357 fprintf(fp, " %6.2f %6.2f\n", (dt == 0) ? 0 : dl->ntr[nDih] / dt, dl->S2[nDih]);
360 void pr_dlist(FILE* fp,
373 void (*pr_props)(FILE*, t_dlist*, int, real);
375 /* Analysis of dihedral transitions etc */
377 if (printtype == edPrintST)
379 pr_props = pr_ntr_s2;
380 fprintf(stderr, "Now printing out transitions and OPs...\n");
384 pr_props = pr_one_ro;
385 fprintf(stderr, "Now printing out rotamer occupancies...\n");
386 fprintf(fp, "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
389 /* change atom numbers from 0 based to 1 based */
390 for (i = 0; (i < nl); i++)
392 fprintf(fp, "Residue %s\n", dl[i].name);
393 if (printtype == edPrintST)
396 " Angle [ AI, AJ, AK, AL] #tr/ns S^2D \n"
397 "--------------------------------------------\n");
402 " Angle [ AI, AJ, AK, AL] rotamers 0 g(-) t g(+)\n"
403 "--------------------------------------------\n");
407 fprintf(fp, " Phi [%5d,%5d,%5d,%5d]",
408 (dl[i].atm.H == -1) ? 1 + dl[i].atm.minC : 1 + dl[i].atm.H, 1 + dl[i].atm.N,
409 1 + dl[i].atm.Cn[1], 1 + dl[i].atm.C);
410 pr_props(fp, &dl[i], edPhi, dt);
414 fprintf(fp, " Psi [%5d,%5d,%5d,%5d]", 1 + dl[i].atm.N, 1 + dl[i].atm.Cn[1],
415 1 + dl[i].atm.C, 1 + dl[i].atm.O);
416 pr_props(fp, &dl[i], edPsi, dt);
418 if (bOmega && has_dihedral(edOmega, &(dl[i])))
420 fprintf(fp, " Omega [%5d,%5d,%5d,%5d]", 1 + dl[i].atm.minCalpha, 1 + dl[i].atm.minC,
421 1 + dl[i].atm.N, 1 + dl[i].atm.Cn[1]);
422 pr_props(fp, &dl[i], edOmega, dt);
424 for (Xi = 0; Xi < MAXCHI; Xi++)
426 if (bChi && (Xi < maxchi) && (dl[i].atm.Cn[Xi + 3] != -1))
428 fprintf(fp, " Chi%d[%5d,%5d,%5d,%5d]", Xi + 1, 1 + dl[i].atm.Cn[Xi],
429 1 + dl[i].atm.Cn[Xi + 1], 1 + dl[i].atm.Cn[Xi + 2], 1 + dl[i].atm.Cn[Xi + 3]);
430 pr_props(fp, &dl[i], Xi + edChi1, dt); /* Xi+2 was wrong here */
438 int pr_trans(FILE* fp, int nl, t_dlist dl[], real dt, int Xi)
440 /* never called at the moment */
445 fprintf(fp, "\\begin{table}[h]\n");
446 fprintf(fp, "\\caption{Number of dihedral transitions per nanosecond}\n");
447 fprintf(fp, "\\begin{tabular}{|l|l|}\n");
448 fprintf(fp, "\\hline\n");
449 fprintf(fp, "Residue\t&$\\chi_%d$\t\\\\\n", Xi + 1);
450 for (i = 0; (i < nl); i++)
452 nn = static_cast<int>(dl[i].ntr[Xi] / dt);
456 fprintf(fp, "%s\t&\\HL{%d}\t\\\\\n", dl[i].name, nn);
461 fprintf(fp, "%s\t&\\%d\t\\\\\n", dl[i].name, nn);
464 fprintf(fp, "\\hline\n");
465 fprintf(fp, "\\end{tabular}\n");
466 fprintf(fp, "\\end{table}\n\n");