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/utility/cstringutil.h"
43 #include "gromacs/legacyheaders/typedefs.h"
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/utility/smalloc.h"
47 #include "gromacs/fileio/enxio.h"
48 #include "gromacs/commandline/pargs.h"
49 #include "gromacs/legacyheaders/names.h"
50 #include "gromacs/legacyheaders/macros.h"
51 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/fileio/matio.h"
55 #include "gromacs/fileio/strdb.h"
57 #include "gromacs/fileio/trxio.h"
60 static int search_str2(int nstr, char **str, char *key)
63 int keylen = strlen(key);
66 while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
70 for (i = 0; (i < nstr); i++)
72 if (gmx_strncasecmp(str[i], key, n) == 0)
81 int gmx_enemat(int argc, char *argv[])
83 const char *desc[] = {
84 "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
85 "With [TT]-groups[tt] a file must be supplied with on each",
86 "line a group of atoms to be used. For these groups matrix of",
87 "interaction energies will be extracted from the energy file",
88 "by looking for energy groups with names corresponding to pairs",
89 "of groups of atoms, e.g. if your [TT]-groups[tt] file contains:[BR]",
91 "[TT]Protein[tt][BR]",
93 "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
94 "'LJ:Protein-SOL' are expected in the energy file (although",
95 "[THISMODULE] is most useful if many groups are analyzed",
96 "simultaneously). Matrices for different energy types are written",
97 "out separately, as controlled by the",
98 "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
99 "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
100 "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
101 "Finally, the total interaction energy energy per group can be ",
102 "calculated ([TT]-etot[tt]).[PAR]",
104 "An approximation of the free energy can be calculated using:",
105 "[MATH]E[SUB]free[sub] = E[SUB]0[sub] + kT [LOG][CHEVRON][EXP](E-E[SUB]0[sub])/kT[exp][chevron][log][math], where '[MATH][CHEVRON][chevron][math]'",
106 "stands for time-average. A file with reference free energies",
107 "can be supplied to calculate the free energy difference",
108 "with some reference state. Group names (e.g. residue names)",
109 "in the reference file should correspond to the group names",
110 "as used in the [TT]-groups[tt] file, but a appended number",
111 "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
114 static gmx_bool bSum = FALSE;
115 static gmx_bool bMeanEmtx = TRUE;
116 static int skip = 0, nlevels = 20;
117 static real cutmax = 1e20, cutmin = -1e20, reftemp = 300.0;
118 static gmx_bool bCoulSR = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
119 static gmx_bool bLJSR = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
122 { "-sum", FALSE, etBOOL, {&bSum},
123 "Sum the energy terms selected rather than display them all" },
124 { "-skip", FALSE, etINT, {&skip},
125 "Skip number of frames between data points" },
126 { "-mean", FALSE, etBOOL, {&bMeanEmtx},
127 "with [TT]-groups[tt] extracts matrix of mean energies instead of "
128 "matrix for each timestep" },
129 { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
130 { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
131 { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
132 { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
133 { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
134 { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
135 { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
136 { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
137 { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
138 { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
139 { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
140 { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
141 { "-temp", FALSE, etREAL, {&reftemp},
142 "reference temperature for free energy calculation"}
144 /* We will define egSP more energy-groups:
145 egTotal (total energy) */
148 gmx_bool egrp_use[egNR+egSP];
152 gmx_enxnm_t *enm = NULL;
156 gmx_bool bCont, bRef;
157 gmx_bool bCutmax, bCutmin;
158 real **eneset, *time = NULL;
159 int *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
160 char **groups = NULL;
161 char groupname[255], fn[255];
163 t_rgb rlo, rhi, rmid;
164 real emax, emid, emin;
165 real ***emat, **etot, *groupnr;
166 double beta, expE, **e, *eaver, *efree = NULL, edum;
168 char **ereflines, **erefres = NULL;
169 real *eref = NULL, *edif = NULL;
174 { efEDR, "-f", NULL, ffOPTRD },
175 { efDAT, "-groups", "groups", ffREAD },
176 { efDAT, "-eref", "eref", ffOPTRD },
177 { efXPM, "-emat", "emat", ffWRITE },
178 { efXVG, "-etot", "energy", ffWRITE }
180 #define NFILE asize(fnm)
182 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
183 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
188 for (i = 0; (i < egNR+egSP); i++)
192 egrp_use[egCOULSR] = bCoulSR;
193 egrp_use[egLJSR] = bLJSR;
194 egrp_use[egBHAMSR] = bBhamSR;
195 egrp_use[egCOULLR] = bCoulLR;
196 egrp_use[egLJLR] = bLJLR;
197 egrp_use[egBHAMLR] = bBhamLR;
198 egrp_use[egCOUL14] = bCoul14;
199 egrp_use[egLJ14] = bLJ14;
200 egrp_use[egTotal] = TRUE;
202 bRef = opt2bSet("-eref", NFILE, fnm);
203 in = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
204 do_enxnms(in, &nre, &enm);
208 gmx_fatal(FARGS, "No energies!\n");
211 bCutmax = opt2parg_bSet("-max", asize(pa), pa);
212 bCutmin = opt2parg_bSet("-min", asize(pa), pa);
216 /* Read groupnames from input file and construct selection of
217 energy groups from it*/
219 fprintf(stderr, "Will read groupnames from inputfile\n");
220 ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
221 fprintf(stderr, "Read %d groups\n", ngroups);
222 snew(set, sqr(ngroups)*egNR/2);
225 for (i = 0; (i < ngroups); i++)
227 fprintf(stderr, "\rgroup %d", i);
228 for (j = i; (j < ngroups); j++)
230 for (m = 0; (m < egNR); m++)
234 sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
236 fprintf(stderr, "\r%-15s %5d", groupname, n);
238 for (k = prevk; (k < prevk+nre); k++)
240 if (strcmp(enm[k%nre].name, groupname) == 0)
248 fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
249 "in energy file\n", groupname, i, j);
259 fprintf(stderr, "\n");
261 snew(eneset, nset+1);
262 fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
264 /* Start reading energy frames */
270 bCont = do_enx(in, fr);
273 timecheck = check_times(fr->t);
276 while (bCont && (timecheck < 0));
280 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
284 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
286 if ((nenergy % 1000) == 0)
288 srenew(time, nenergy+1000);
289 for (i = 0; (i <= nset); i++)
291 srenew(eneset[i], nenergy+1000);
294 time[nenergy] = fr->t;
296 for (i = 0; (i < nset); i++)
298 eneset[i][nenergy] = fr->ener[set[i]].e;
299 sum += fr->ener[set[i]].e;
303 eneset[nset][nenergy] = sum;
310 while (bCont && (timecheck == 0));
312 fprintf(stderr, "\n");
314 fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
315 "over %d frames\n", ngroups, nset, nenergy);
317 snew(emat, egNR+egSP);
318 for (j = 0; (j < egNR+egSP); j++)
322 snew(emat[j], ngroups);
323 for (i = 0; (i < ngroups); i++)
325 snew(emat[j][i], ngroups);
329 snew(groupnr, ngroups);
330 for (i = 0; (i < ngroups); i++)
334 rlo.r = 1.0, rlo.g = 0.0, rlo.b = 0.0;
335 rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
336 rhi.r = 0.0, rhi.g = 0.0, rhi.b = 1.0;
340 for (i = 0; (i < ngroups); i++)
345 for (i = 0; (i < ngroups); i++)
347 for (j = i; (j < ngroups); j++)
349 for (m = 0; (m < egNR); m++)
353 for (k = 0; (k < nenergy); k++)
355 emat[m][i][j] += eneset[n][k];
356 e[i][k] += eneset[n][k]; /* *0.5; */
357 e[j][k] += eneset[n][k]; /* *0.5; */
360 emat[egTotal][i][j] += emat[m][i][j];
361 emat[m][i][j] /= nenergy;
362 emat[m][j][i] = emat[m][i][j];
365 emat[egTotal][i][j] /= nenergy;
366 emat[egTotal][j][i] = emat[egTotal][i][j];
373 fprintf(stderr, "Will read reference energies from inputfile\n");
374 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
375 fprintf(stderr, "Read %d reference energies\n", neref);
377 snew(erefres, neref);
378 for (i = 0; (i < neref); i++)
381 sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
385 snew(eaver, ngroups);
386 for (i = 0; (i < ngroups); i++)
388 for (k = 0; (k < nenergy); k++)
394 beta = 1.0/(BOLTZ*reftemp);
395 snew(efree, ngroups);
397 for (i = 0; (i < ngroups); i++)
400 for (k = 0; (k < nenergy); k++)
402 expE += exp(beta*(e[i][k]-eaver[i]));
404 efree[i] = log(expE/nenergy)/beta + eaver[i];
407 n = search_str2(neref, erefres, groups[i]);
410 edif[i] = efree[i]-eref[n];
415 fprintf(stderr, "WARNING: group %s not found "
416 "in reference energies.\n", groups[i]);
426 emid = 0.0; /*(emin+emax)*0.5;*/
427 egrp_nm[egTotal] = "total";
428 for (m = 0; (m < egNR+egSP); m++)
434 for (i = 0; (i < ngroups); i++)
436 for (j = i; (j < ngroups); j++)
438 if (emat[m][i][j] > emax)
440 emax = emat[m][i][j];
442 else if (emat[m][i][j] < emin)
444 emin = emat[m][i][j];
450 fprintf(stderr, "Matrix of %s energy is uniform at %f "
451 "(will not produce output).\n", egrp_nm[m], emax);
455 fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
456 egrp_nm[m], emin, emax);
457 if ((bCutmax) || (emax > cutmax))
461 if ((bCutmin) || (emin < cutmin))
465 if ((emax == cutmax) || (emin == cutmin))
467 fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
470 sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
471 sprintf(label, "%s Interaction Energies", egrp_nm[m]);
472 out = gmx_ffopen(fn, "w");
475 write_xpm(out, 0, label, "Energy (kJ/mol)",
476 "Residue Index", "Residue Index",
477 ngroups, ngroups, groupnr, groupnr, emat[m],
478 emid, emax, rmid, rhi, &nlevels);
480 else if (emax <= emid)
482 write_xpm(out, 0, label, "Energy (kJ/mol)",
483 "Residue Index", "Residue Index",
484 ngroups, ngroups, groupnr, groupnr, emat[m],
485 emin, emid, rlo, rmid, &nlevels);
489 write_xpm3(out, 0, label, "Energy (kJ/mol)",
490 "Residue Index", "Residue Index",
491 ngroups, ngroups, groupnr, groupnr, emat[m],
492 emin, emid, emax, rlo, rmid, rhi, &nlevels);
498 snew(etot, egNR+egSP);
499 for (m = 0; (m < egNR+egSP); m++)
501 snew(etot[m], ngroups);
502 for (i = 0; (i < ngroups); i++)
504 for (j = 0; (j < ngroups); j++)
506 etot[m][i] += emat[m][i][j];
511 out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
513 xvgr_legend(out, 0, NULL, oenv);
515 if (output_env_get_print_xvgr_codes(oenv))
517 char str1[STRLEN], str2[STRLEN];
518 if (output_env_get_xvg_format(oenv) == exvgXMGR)
520 sprintf(str1, "@ legend string ");
525 sprintf(str1, "@ s");
526 sprintf(str2, " legend ");
529 for (m = 0; (m < egNR+egSP); m++)
533 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
538 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
542 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
544 fprintf(out, "@TYPE xy\n");
545 fprintf(out, "#%3s", "grp");
547 for (m = 0; (m < egNR+egSP); m++)
551 fprintf(out, " %9s", egrp_nm[m]);
556 fprintf(out, " %9s", "Free");
560 fprintf(out, " %9s", "Diff");
564 for (i = 0; (i < ngroups); i++)
566 fprintf(out, "%3.0f", groupnr[i]);
567 for (m = 0; (m < egNR+egSP); m++)
571 fprintf(out, " %9.5g", etot[m][i]);
576 fprintf(out, " %9.5g", efree[i]);
580 fprintf(out, " %9.5g", edif[i]);
588 fprintf(stderr, "While typing at your keyboard, suddenly...\n"
589 "...nothing happens.\nWARNING: Not Implemented Yet\n");
591 out=ftp2FILE(efMAT,NFILE,fnm,"w");
594 for (k=0; (k<nenergy); k++) {
595 for (i=0; (i<ngroups); i++)
596 for (j=i+1; (j<ngroups); j++)
597 emat[i][j]=eneset[n][k];
598 sprintf(label,"t=%.0f ps",time[k]);
599 write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);