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.
43 #include "gromacs/utility/cstringutil.h"
45 #include "gmx_fatal.h"
47 #include "gromacs/utility/smalloc.h"
48 #include "gromacs/fileio/enxio.h"
49 #include "gromacs/commandline/pargs.h"
55 #include "gromacs/fileio/matio.h"
56 #include "gromacs/fileio/strdb.h"
58 #include "gromacs/fileio/trxio.h"
61 static int search_str2(int nstr, char **str, char *key)
64 int keylen = strlen(key);
67 while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
71 for (i = 0; (i < nstr); i++)
73 if (gmx_strncasecmp(str[i], key, n) == 0)
82 int gmx_enemat(int argc, char *argv[])
84 const char *desc[] = {
85 "[THISMODULE] extracts an energy matrix from the energy file ([TT]-f[tt]).",
86 "With [TT]-groups[tt] a file must be supplied with on each",
87 "line a group of atoms to be used. For these groups matrix of",
88 "interaction energies will be extracted from the energy file",
89 "by looking for energy groups with names corresponding to pairs",
90 "of groups of atoms, e.g. if your [TT]-groups[tt] file contains:[BR]",
92 "[TT]Protein[tt][BR]",
94 "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
95 "'LJ:Protein-SOL' are expected in the energy file (although",
96 "[THISMODULE] is most useful if many groups are analyzed",
97 "simultaneously). Matrices for different energy types are written",
98 "out separately, as controlled by the",
99 "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
100 "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
101 "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
102 "Finally, the total interaction energy energy per group can be ",
103 "calculated ([TT]-etot[tt]).[PAR]",
105 "An approximation of the free energy can be calculated using:",
106 "[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]'",
107 "stands for time-average. A file with reference free energies",
108 "can be supplied to calculate the free energy difference",
109 "with some reference state. Group names (e.g. residue names)",
110 "in the reference file should correspond to the group names",
111 "as used in the [TT]-groups[tt] file, but a appended number",
112 "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
115 static gmx_bool bSum = FALSE;
116 static gmx_bool bMeanEmtx = TRUE;
117 static int skip = 0, nlevels = 20;
118 static real cutmax = 1e20, cutmin = -1e20, reftemp = 300.0;
119 static gmx_bool bCoulSR = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
120 static gmx_bool bLJSR = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
123 { "-sum", FALSE, etBOOL, {&bSum},
124 "Sum the energy terms selected rather than display them all" },
125 { "-skip", FALSE, etINT, {&skip},
126 "Skip number of frames between data points" },
127 { "-mean", FALSE, etBOOL, {&bMeanEmtx},
128 "with [TT]-groups[tt] extracts matrix of mean energies instead of "
129 "matrix for each timestep" },
130 { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
131 { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
132 { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
133 { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
134 { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
135 { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
136 { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
137 { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
138 { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
139 { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
140 { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
141 { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
142 { "-temp", FALSE, etREAL, {&reftemp},
143 "reference temperature for free energy calculation"}
145 /* We will define egSP more energy-groups:
146 egTotal (total energy) */
149 gmx_bool egrp_use[egNR+egSP];
153 gmx_enxnm_t *enm = NULL;
157 gmx_bool bCont, bRef;
158 gmx_bool bCutmax, bCutmin;
159 real **eneset, *time = NULL;
160 int *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
161 char **groups = NULL;
162 char groupname[255], fn[255];
164 t_rgb rlo, rhi, rmid;
165 real emax, emid, emin;
166 real ***emat, **etot, *groupnr;
167 double beta, expE, **e, *eaver, *efree = NULL, edum;
169 char **ereflines, **erefres = NULL;
170 real *eref = NULL, *edif = NULL;
175 { efEDR, "-f", NULL, ffOPTRD },
176 { efDAT, "-groups", "groups.dat", ffREAD },
177 { efDAT, "-eref", "eref.dat", ffOPTRD },
178 { efXPM, "-emat", "emat", ffWRITE },
179 { efXVG, "-etot", "energy", ffWRITE }
181 #define NFILE asize(fnm)
183 if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
184 NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
189 for (i = 0; (i < egNR+egSP); i++)
193 egrp_use[egCOULSR] = bCoulSR;
194 egrp_use[egLJSR] = bLJSR;
195 egrp_use[egBHAMSR] = bBhamSR;
196 egrp_use[egCOULLR] = bCoulLR;
197 egrp_use[egLJLR] = bLJLR;
198 egrp_use[egBHAMLR] = bBhamLR;
199 egrp_use[egCOUL14] = bCoul14;
200 egrp_use[egLJ14] = bLJ14;
201 egrp_use[egTotal] = TRUE;
203 bRef = opt2bSet("-eref", NFILE, fnm);
204 in = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
205 do_enxnms(in, &nre, &enm);
209 gmx_fatal(FARGS, "No energies!\n");
212 bCutmax = opt2parg_bSet("-max", asize(pa), pa);
213 bCutmin = opt2parg_bSet("-min", asize(pa), pa);
217 /* Read groupnames from input file and construct selection of
218 energy groups from it*/
220 fprintf(stderr, "Will read groupnames from inputfile\n");
221 ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
222 fprintf(stderr, "Read %d groups\n", ngroups);
223 snew(set, sqr(ngroups)*egNR/2);
226 for (i = 0; (i < ngroups); i++)
228 fprintf(stderr, "\rgroup %d", i);
229 for (j = i; (j < ngroups); j++)
231 for (m = 0; (m < egNR); m++)
235 sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
237 fprintf(stderr, "\r%-15s %5d", groupname, n);
239 for (k = prevk; (k < prevk+nre); k++)
241 if (strcmp(enm[k%nre].name, groupname) == 0)
249 fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
250 "in energy file\n", groupname, i, j);
260 fprintf(stderr, "\n");
262 snew(eneset, nset+1);
263 fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
265 /* Start reading energy frames */
271 bCont = do_enx(in, fr);
274 timecheck = check_times(fr->t);
277 while (bCont && (timecheck < 0));
281 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
285 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
287 if ((nenergy % 1000) == 0)
289 srenew(time, nenergy+1000);
290 for (i = 0; (i <= nset); i++)
292 srenew(eneset[i], nenergy+1000);
295 time[nenergy] = fr->t;
297 for (i = 0; (i < nset); i++)
299 eneset[i][nenergy] = fr->ener[set[i]].e;
300 sum += fr->ener[set[i]].e;
304 eneset[nset][nenergy] = sum;
311 while (bCont && (timecheck == 0));
313 fprintf(stderr, "\n");
315 fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
316 "over %d frames\n", ngroups, nset, nenergy);
318 snew(emat, egNR+egSP);
319 for (j = 0; (j < egNR+egSP); j++)
323 snew(emat[j], ngroups);
324 for (i = 0; (i < ngroups); i++)
326 snew(emat[j][i], ngroups);
330 snew(groupnr, ngroups);
331 for (i = 0; (i < ngroups); i++)
335 rlo.r = 1.0, rlo.g = 0.0, rlo.b = 0.0;
336 rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
337 rhi.r = 0.0, rhi.g = 0.0, rhi.b = 1.0;
341 for (i = 0; (i < ngroups); i++)
346 for (i = 0; (i < ngroups); i++)
348 for (j = i; (j < ngroups); j++)
350 for (m = 0; (m < egNR); m++)
354 for (k = 0; (k < nenergy); k++)
356 emat[m][i][j] += eneset[n][k];
357 e[i][k] += eneset[n][k]; /* *0.5; */
358 e[j][k] += eneset[n][k]; /* *0.5; */
361 emat[egTotal][i][j] += emat[m][i][j];
362 emat[m][i][j] /= nenergy;
363 emat[m][j][i] = emat[m][i][j];
366 emat[egTotal][i][j] /= nenergy;
367 emat[egTotal][j][i] = emat[egTotal][i][j];
374 fprintf(stderr, "Will read reference energies from inputfile\n");
375 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
376 fprintf(stderr, "Read %d reference energies\n", neref);
378 snew(erefres, neref);
379 for (i = 0; (i < neref); i++)
382 sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
386 snew(eaver, ngroups);
387 for (i = 0; (i < ngroups); i++)
389 for (k = 0; (k < nenergy); k++)
395 beta = 1.0/(BOLTZ*reftemp);
396 snew(efree, ngroups);
398 for (i = 0; (i < ngroups); i++)
401 for (k = 0; (k < nenergy); k++)
403 expE += exp(beta*(e[i][k]-eaver[i]));
405 efree[i] = log(expE/nenergy)/beta + eaver[i];
408 n = search_str2(neref, erefres, groups[i]);
411 edif[i] = efree[i]-eref[n];
416 fprintf(stderr, "WARNING: group %s not found "
417 "in reference energies.\n", groups[i]);
427 emid = 0.0; /*(emin+emax)*0.5;*/
428 egrp_nm[egTotal] = "total";
429 for (m = 0; (m < egNR+egSP); m++)
435 for (i = 0; (i < ngroups); i++)
437 for (j = i; (j < ngroups); j++)
439 if (emat[m][i][j] > emax)
441 emax = emat[m][i][j];
443 else if (emat[m][i][j] < emin)
445 emin = emat[m][i][j];
451 fprintf(stderr, "Matrix of %s energy is uniform at %f "
452 "(will not produce output).\n", egrp_nm[m], emax);
456 fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
457 egrp_nm[m], emin, emax);
458 if ((bCutmax) || (emax > cutmax))
462 if ((bCutmin) || (emin < cutmin))
466 if ((emax == cutmax) || (emin == cutmin))
468 fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
471 sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
472 sprintf(label, "%s Interaction Energies", egrp_nm[m]);
473 out = gmx_ffopen(fn, "w");
476 write_xpm(out, 0, label, "Energy (kJ/mol)",
477 "Residue Index", "Residue Index",
478 ngroups, ngroups, groupnr, groupnr, emat[m],
479 emid, emax, rmid, rhi, &nlevels);
481 else if (emax <= emid)
483 write_xpm(out, 0, label, "Energy (kJ/mol)",
484 "Residue Index", "Residue Index",
485 ngroups, ngroups, groupnr, groupnr, emat[m],
486 emin, emid, rlo, rmid, &nlevels);
490 write_xpm3(out, 0, label, "Energy (kJ/mol)",
491 "Residue Index", "Residue Index",
492 ngroups, ngroups, groupnr, groupnr, emat[m],
493 emin, emid, emax, rlo, rmid, rhi, &nlevels);
499 snew(etot, egNR+egSP);
500 for (m = 0; (m < egNR+egSP); m++)
502 snew(etot[m], ngroups);
503 for (i = 0; (i < ngroups); i++)
505 for (j = 0; (j < ngroups); j++)
507 etot[m][i] += emat[m][i][j];
512 out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
514 xvgr_legend(out, 0, NULL, oenv);
516 if (output_env_get_print_xvgr_codes(oenv))
518 char str1[STRLEN], str2[STRLEN];
519 if (output_env_get_xvg_format(oenv) == exvgXMGR)
521 sprintf(str1, "@ legend string ");
526 sprintf(str1, "@ s");
527 sprintf(str2, " legend ");
530 for (m = 0; (m < egNR+egSP); m++)
534 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
539 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
543 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
545 fprintf(out, "@TYPE xy\n");
546 fprintf(out, "#%3s", "grp");
548 for (m = 0; (m < egNR+egSP); m++)
552 fprintf(out, " %9s", egrp_nm[m]);
557 fprintf(out, " %9s", "Free");
561 fprintf(out, " %9s", "Diff");
565 for (i = 0; (i < ngroups); i++)
567 fprintf(out, "%3.0f", groupnr[i]);
568 for (m = 0; (m < egNR+egSP); m++)
572 fprintf(out, " %9.5g", etot[m][i]);
577 fprintf(out, " %9.5g", efree[i]);
581 fprintf(out, " %9.5g", edif[i]);
589 fprintf(stderr, "While typing at your keyboard, suddenly...\n"
590 "...nothing happens.\nWARNING: Not Implemented Yet\n");
592 out=ftp2FILE(efMAT,NFILE,fnm,"w");
595 for (k=0; (k<nenergy); k++) {
596 for (i=0; (i<ngroups); i++)
597 for (j=i+1; (j<ngroups); j++)
598 emat[i][j]=eneset[n][k];
599 sprintf(label,"t=%.0f ps",time[k]);
600 write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);