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,2017,2018,2019, 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/commandline/pargs.h"
43 #include "gromacs/fileio/enxio.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/trxio.h"
46 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/gmxana/gmx_ana.h"
48 #include "gromacs/gmxana/gstat.h"
49 #include "gromacs/math/functions.h"
50 #include "gromacs/math/units.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdlib/energyoutput.h"
53 #include "gromacs/mdtypes/enerdata.h"
54 #include "gromacs/mdtypes/md_enums.h"
55 #include "gromacs/trajectory/energyframe.h"
56 #include "gromacs/utility/arraysize.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
61 #include "gromacs/utility/strdb.h"
64 static int search_str2(int nstr, char** str, char* key)
67 int keylen = std::strlen(key);
70 while ((n < keylen) && ((key[n] < '0') || (key[n] > '9')))
74 for (i = 0; (i < nstr); i++)
76 if (gmx_strncasecmp(str[i], key, n) == 0)
85 int gmx_enemat(int argc, char* argv[])
87 const char* desc[] = {
88 "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
89 "With [TT]-groups[tt] a file must be supplied with on each",
90 "line a group of atoms to be used. For these groups matrix of",
91 "interaction energies will be extracted from the energy file",
92 "by looking for energy groups with names corresponding to pairs",
93 "of groups of atoms, e.g. if your [TT]-groups[tt] file contains::",
99 "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
100 "'LJ:Protein-SOL' are expected in the energy file (although",
101 "[THISMODULE] is most useful if many groups are analyzed",
102 "simultaneously). Matrices for different energy types are written",
103 "out separately, as controlled by the",
104 "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
105 "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
106 "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
107 "Finally, the total interaction energy energy per group can be ",
108 "calculated ([TT]-etot[tt]).[PAR]",
110 "An approximation of the free energy can be calculated using:",
111 "[MATH]E[SUB]free[sub] = E[SUB]0[sub] + kT ",
112 "[LOG][CHEVRON][EXP](E-E[SUB]0[sub])/kT[exp][chevron][log][math], where ",
113 "'[MATH][CHEVRON][chevron][math]'",
114 "stands for time-average. A file with reference free energies",
115 "can be supplied to calculate the free energy difference",
116 "with some reference state. Group names (e.g. residue names)",
117 "in the reference file should correspond to the group names",
118 "as used in the [TT]-groups[tt] file, but a appended number",
119 "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
122 static gmx_bool bSum = FALSE;
123 static gmx_bool bMeanEmtx = TRUE;
124 static int skip = 0, nlevels = 20;
125 static real cutmax = 1e20, cutmin = -1e20, reftemp = 300.0;
126 static gmx_bool bCoulSR = TRUE, bCoul14 = FALSE;
127 static gmx_bool bLJSR = TRUE, bLJ14 = FALSE, bBhamSR = FALSE, bFree = TRUE;
133 "Sum the energy terms selected rather than display them all" },
134 { "-skip", FALSE, etINT, { &skip }, "Skip number of frames between data points" },
139 "with [TT]-groups[tt] extracts matrix of mean energies instead of "
140 "matrix for each timestep" },
141 { "-nlevels", FALSE, etINT, { &nlevels }, "number of levels for matrix colors" },
142 { "-max", FALSE, etREAL, { &cutmax }, "max value for energies" },
143 { "-min", FALSE, etREAL, { &cutmin }, "min value for energies" },
144 { "-coulsr", FALSE, etBOOL, { &bCoulSR }, "extract Coulomb SR energies" },
145 { "-coul14", FALSE, etBOOL, { &bCoul14 }, "extract Coulomb 1-4 energies" },
146 { "-ljsr", FALSE, etBOOL, { &bLJSR }, "extract Lennard-Jones SR energies" },
147 { "-lj14", FALSE, etBOOL, { &bLJ14 }, "extract Lennard-Jones 1-4 energies" },
148 { "-bhamsr", FALSE, etBOOL, { &bBhamSR }, "extract Buckingham SR energies" },
149 { "-free", FALSE, etBOOL, { &bFree }, "calculate free energy" },
154 "reference temperature for free energy calculation" }
156 /* We will define egSP more energy-groups:
157 egTotal (total energy) */
160 gmx_bool egrp_use[egNR + egSP];
164 gmx_enxnm_t* enm = nullptr;
168 gmx_bool bCont, bRef;
169 gmx_bool bCutmax, bCutmin;
170 real ** eneset, *time = nullptr;
171 int * set, i, j, prevk, k, m = 0, n, nre, nset, nenergy;
172 char** groups = nullptr;
173 char groupname[255], fn[255];
175 t_rgb rlo, rhi, rmid;
176 real emax, emid, emin;
177 real *** emat, **etot, *groupnr;
178 double beta, expE, **e, *eaver, *efree = nullptr, edum;
180 char ** ereflines, **erefres = nullptr;
181 real * eref = nullptr, *edif = nullptr;
183 gmx_output_env_t* oenv;
185 t_filenm fnm[] = { { efEDR, "-f", nullptr, ffOPTRD },
186 { efDAT, "-groups", "groups", ffREAD },
187 { efDAT, "-eref", "eref", ffOPTRD },
188 { efXPM, "-emat", "emat", ffWRITE },
189 { efXVG, "-etot", "energy", ffWRITE } };
190 #define NFILE asize(fnm)
192 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME, NFILE, fnm, asize(pa), pa,
193 asize(desc), desc, 0, nullptr, &oenv))
198 for (i = 0; (i < egNR + egSP); i++)
202 egrp_use[egCOULSR] = bCoulSR;
203 egrp_use[egLJSR] = bLJSR;
204 egrp_use[egBHAMSR] = bBhamSR;
205 egrp_use[egCOUL14] = bCoul14;
206 egrp_use[egLJ14] = bLJ14;
207 egrp_use[egTotal] = TRUE;
209 bRef = opt2bSet("-eref", NFILE, fnm);
210 in = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
211 do_enxnms(in, &nre, &enm);
215 gmx_fatal(FARGS, "No energies!\n");
218 bCutmax = opt2parg_bSet("-max", asize(pa), pa);
219 bCutmin = opt2parg_bSet("-min", asize(pa), pa);
223 /* Read groupnames from input file and construct selection of
224 energy groups from it*/
226 fprintf(stderr, "Will read groupnames from inputfile\n");
227 ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
228 fprintf(stderr, "Read %d groups\n", ngroups);
229 snew(set, static_cast<size_t>(gmx::square(ngroups) * egNR / 2));
232 for (i = 0; (i < ngroups); i++)
234 for (j = i; (j < ngroups); j++)
236 for (m = 0; (m < egNR); m++)
240 sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
241 bool foundMatch = false;
242 for (k = prevk; (k < prevk + nre); k++)
244 if (std::strcmp(enm[k % nre].name, groupname) == 0)
254 "WARNING! could not find group %s (%d,%d) "
266 fprintf(stderr, "\n");
269 // Return an error, can't do what the user asked for
271 "None of the specified energy groups were found in this .edr file.\n"
272 "Perhaps you used the wrong groups, the wrong files, or didn't use a .tpr\n"
273 "that was made from an .mdp file that specified these energy groups.\n");
277 snew(eneset, nset + 1);
278 fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
280 /* Start reading energy frames */
286 bCont = do_enx(in, fr);
289 timecheck = check_times(fr->t);
291 } while (bCont && (timecheck < 0));
295 #define DONTSKIP(cnt) (skip) ? (((cnt) % skip) == 0) : TRUE
299 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
302 if ((nenergy % 1000) == 0)
304 srenew(time, nenergy + 1000);
305 for (i = 0; (i <= nset); i++)
307 srenew(eneset[i], nenergy + 1000);
310 time[nenergy] = fr->t;
312 for (i = 0; (i < nset); i++)
314 eneset[i][nenergy] = fr->ener[set[i]].e;
315 sum += fr->ener[set[i]].e;
319 eneset[nset][nenergy] = sum;
325 } while (bCont && (timecheck == 0));
327 fprintf(stderr, "\n");
330 "Will build energy half-matrix of %d groups, %d elements, "
332 ngroups, nset, nenergy);
334 snew(emat, egNR + egSP);
335 for (j = 0; (j < egNR + egSP); j++)
339 snew(emat[j], ngroups);
340 for (i = 0; (i < ngroups); i++)
342 snew(emat[j][i], ngroups);
346 snew(groupnr, ngroups);
347 for (i = 0; (i < ngroups); i++)
363 for (i = 0; (i < ngroups); i++)
368 for (i = 0; (i < ngroups); i++)
370 for (j = i; (j < ngroups); j++)
372 for (m = 0; (m < egNR); m++)
376 for (k = 0; (k < nenergy); k++)
378 emat[m][i][j] += eneset[n][k];
379 e[i][k] += eneset[n][k]; /* *0.5; */
380 e[j][k] += eneset[n][k]; /* *0.5; */
383 emat[egTotal][i][j] += emat[m][i][j];
384 emat[m][i][j] /= nenergy;
385 emat[m][j][i] = emat[m][i][j];
388 emat[egTotal][i][j] /= nenergy;
389 emat[egTotal][j][i] = emat[egTotal][i][j];
396 fprintf(stderr, "Will read reference energies from inputfile\n");
397 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
398 fprintf(stderr, "Read %d reference energies\n", neref);
400 snew(erefres, neref);
401 for (i = 0; (i < neref); i++)
404 sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
408 snew(eaver, ngroups);
409 for (i = 0; (i < ngroups); i++)
411 for (k = 0; (k < nenergy); k++)
417 beta = 1.0 / (BOLTZ * reftemp);
418 snew(efree, ngroups);
420 for (i = 0; (i < ngroups); i++)
423 for (k = 0; (k < nenergy); k++)
425 expE += std::exp(beta * (e[i][k] - eaver[i]));
427 efree[i] = std::log(expE / nenergy) / beta + eaver[i];
430 n = search_str2(neref, erefres, groups[i]);
433 edif[i] = efree[i] - eref[n];
439 "WARNING: group %s not found "
440 "in reference energies.\n",
451 emid = 0.0; /*(emin+emax)*0.5;*/
452 egrp_nm[egTotal] = "total";
453 for (m = 0; (m < egNR + egSP); m++)
459 for (i = 0; (i < ngroups); i++)
461 for (j = i; (j < ngroups); j++)
463 if (emat[m][i][j] > emax)
465 emax = emat[m][i][j];
467 else if (emat[m][i][j] < emin)
469 emin = emat[m][i][j];
476 "Matrix of %s energy is uniform at %f "
477 "(will not produce output).\n",
482 fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n", egrp_nm[m], emin, emax);
483 if ((bCutmax) || (emax > cutmax))
487 if ((bCutmin) || (emin < cutmin))
491 if ((emax == cutmax) || (emin == cutmin))
493 fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
496 sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
497 sprintf(label, "%s Interaction Energies", egrp_nm[m]);
498 out = gmx_ffopen(fn, "w");
501 write_xpm(out, 0, label, "Energy (kJ/mol)", "Residue Index",
502 "Residue Index", ngroups, ngroups, groupnr, groupnr, emat[m],
503 emid, emax, rmid, rhi, &nlevels);
505 else if (emax <= emid)
507 write_xpm(out, 0, label, "Energy (kJ/mol)", "Residue Index",
508 "Residue Index", ngroups, ngroups, groupnr, groupnr, emat[m],
509 emin, emid, rlo, rmid, &nlevels);
513 write_xpm3(out, 0, label, "Energy (kJ/mol)", "Residue Index",
514 "Residue Index", ngroups, ngroups, groupnr, groupnr, emat[m],
515 emin, emid, emax, rlo, rmid, rhi, &nlevels);
521 snew(etot, egNR + egSP);
522 for (m = 0; (m < egNR + egSP); m++)
524 snew(etot[m], ngroups);
525 for (i = 0; (i < ngroups); i++)
527 for (j = 0; (j < ngroups); j++)
529 etot[m][i] += emat[m][i][j];
534 out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol", oenv);
535 xvgr_legend(out, 0, nullptr, oenv);
537 if (output_env_get_print_xvgr_codes(oenv))
539 char str1[STRLEN], str2[STRLEN];
540 if (output_env_get_xvg_format(oenv) == exvgXMGR)
542 sprintf(str1, "@ legend string ");
547 sprintf(str1, "@ s");
548 sprintf(str2, " legend ");
551 for (m = 0; (m < egNR + egSP); m++)
555 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
560 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
564 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
566 fprintf(out, "@TYPE xy\n");
567 fprintf(out, "#%3s", "grp");
569 for (m = 0; (m < egNR + egSP); m++)
573 fprintf(out, " %9s", egrp_nm[m]);
578 fprintf(out, " %9s", "Free");
582 fprintf(out, " %9s", "Diff");
586 for (i = 0; (i < ngroups); i++)
588 fprintf(out, "%3.0f", groupnr[i]);
589 for (m = 0; (m < egNR + egSP); m++)
593 fprintf(out, " %9.5g", etot[m][i]);
598 fprintf(out, " %9.5g", efree[i]);
602 fprintf(out, " %9.5g", edif[i]);
611 "While typing at your keyboard, suddenly...\n"
612 "...nothing happens.\nWARNING: Not Implemented Yet\n");
614 out=ftp2FILE(efMAT,NFILE,fnm,"w");
617 for (k=0; (k<nenergy); k++) {
618 for (i=0; (i<ngroups); i++)
619 for (j=i+1; (j<ngroups); j++)
620 emat[i][j]=eneset[n][k];
621 sprintf(label,"t=%.0f ps",time[k]);
622 write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);