Moved first batch of analysis tool source to C++
[alexxy/gromacs.git] / src / gromacs / gmxana / gmx_enemat.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
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, 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.
10  *
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.
15  *
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.
20  *
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.
25  *
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.
33  *
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.
36  */
37 #include "gmxpre.h"
38
39 #include <cmath>
40 #include <cstring>
41
42 #include "gromacs/commandline/pargs.h"
43 #include "gromacs/fileio/enxio.h"
44 #include "gromacs/fileio/matio.h"
45 #include "gromacs/fileio/strdb.h"
46 #include "gromacs/fileio/trxio.h"
47 #include "gromacs/fileio/xvgr.h"
48 #include "gromacs/gmxana/gmx_ana.h"
49 #include "gromacs/gmxana/gstat.h"
50 #include "gromacs/legacyheaders/macros.h"
51 #include "gromacs/legacyheaders/names.h"
52 #include "gromacs/legacyheaders/typedefs.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/futil.h"
58 #include "gromacs/utility/smalloc.h"
59
60
61 static int search_str2(int nstr, char **str, char *key)
62 {
63     int  i, n;
64     int  keylen = std::strlen(key);
65     /* Linear search */
66     n = 0;
67     while ( (n < keylen) && ((key[n] < '0') || (key[n] > '9')) )
68     {
69         n++;
70     }
71     for (i = 0; (i < nstr); i++)
72     {
73         if (gmx_strncasecmp(str[i], key, n) == 0)
74         {
75             return i;
76         }
77     }
78
79     return -1;
80 }
81
82 int gmx_enemat(int argc, char *argv[])
83 {
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::",
91         "",
92         "    2",
93         "    Protein",
94         "    SOL",
95         "",
96         "then energy groups with names like 'Coul-SR:Protein-SOL' and ",
97         "'LJ:Protein-SOL' are expected in the energy file (although",
98         "[THISMODULE] is most useful if many groups are analyzed",
99         "simultaneously). Matrices for different energy types are written",
100         "out separately, as controlled by the",
101         "[TT]-[no]coul[tt], [TT]-[no]coulr[tt], [TT]-[no]coul14[tt], ",
102         "[TT]-[no]lj[tt], [TT]-[no]lj14[tt], ",
103         "[TT]-[no]bham[tt] and [TT]-[no]free[tt] options.",
104         "Finally, the total interaction energy energy per group can be ",
105         "calculated ([TT]-etot[tt]).[PAR]",
106
107         "An approximation of the free energy can be calculated using:",
108         "[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]'",
109         "stands for time-average. A file with reference free energies",
110         "can be supplied to calculate the free energy difference",
111         "with some reference state. Group names (e.g. residue names)",
112         "in the reference file should correspond to the group names",
113         "as used in the [TT]-groups[tt] file, but a appended number",
114         "(e.g. residue number) in the [TT]-groups[tt] will be ignored",
115         "in the comparison."
116     };
117     static gmx_bool bSum      = FALSE;
118     static gmx_bool bMeanEmtx = TRUE;
119     static int      skip      = 0, nlevels = 20;
120     static real     cutmax    = 1e20, cutmin = -1e20, reftemp = 300.0;
121     static gmx_bool bCoulSR   = TRUE, bCoulLR = FALSE, bCoul14 = FALSE;
122     static gmx_bool bLJSR     = TRUE, bLJLR = FALSE, bLJ14 = FALSE, bBhamSR = FALSE, bBhamLR = FALSE,
123                     bFree     = TRUE;
124     t_pargs         pa[]      = {
125         { "-sum",  FALSE, etBOOL, {&bSum},
126           "Sum the energy terms selected rather than display them all" },
127         { "-skip", FALSE, etINT,  {&skip},
128           "Skip number of frames between data points" },
129         { "-mean", FALSE, etBOOL, {&bMeanEmtx},
130           "with [TT]-groups[tt] extracts matrix of mean energies instead of "
131           "matrix for each timestep" },
132         { "-nlevels", FALSE, etINT, {&nlevels}, "number of levels for matrix colors"},
133         { "-max", FALSE, etREAL, {&cutmax}, "max value for energies"},
134         { "-min", FALSE, etREAL, {&cutmin}, "min value for energies"},
135         { "-coulsr", FALSE, etBOOL, {&bCoulSR}, "extract Coulomb SR energies"},
136         { "-coullr", FALSE, etBOOL, {&bCoulLR}, "extract Coulomb LR energies"},
137         { "-coul14", FALSE, etBOOL, {&bCoul14}, "extract Coulomb 1-4 energies"},
138         { "-ljsr", FALSE, etBOOL, {&bLJSR}, "extract Lennard-Jones SR energies"},
139         { "-ljlr", FALSE, etBOOL, {&bLJLR}, "extract Lennard-Jones LR energies"},
140         { "-lj14", FALSE, etBOOL, {&bLJ14}, "extract Lennard-Jones 1-4 energies"},
141         { "-bhamsr", FALSE, etBOOL, {&bBhamSR}, "extract Buckingham SR energies"},
142         { "-bhamlr", FALSE, etBOOL, {&bBhamLR}, "extract Buckingham LR energies"},
143         { "-free", FALSE, etBOOL, {&bFree}, "calculate free energy"},
144         { "-temp", FALSE, etREAL, {&reftemp},
145           "reference temperature for free energy calculation"}
146     };
147     /* We will define egSP more energy-groups:
148        egTotal (total energy) */
149 #define egTotal egNR
150 #define egSP 1
151     gmx_bool       egrp_use[egNR+egSP];
152     ener_file_t    in;
153     FILE          *out;
154     int            timecheck = 0;
155     gmx_enxnm_t   *enm       = NULL;
156     t_enxframe    *fr;
157     int            teller = 0;
158     real           sum;
159     gmx_bool       bCont, bRef;
160     gmx_bool       bCutmax, bCutmin;
161     real         **eneset, *time = NULL;
162     int           *set, i, j, k, prevk, m = 0, n, nre, nset, nenergy;
163     char         **groups = NULL;
164     char           groupname[255], fn[255];
165     int            ngroups;
166     t_rgb          rlo, rhi, rmid;
167     real           emax, emid, emin;
168     real        ***emat, **etot, *groupnr;
169     double         beta, expE, **e, *eaver, *efree = NULL, edum;
170     char           label[234];
171     char         **ereflines, **erefres = NULL;
172     real          *eref  = NULL, *edif = NULL;
173     int            neref = 0;
174     output_env_t   oenv;
175
176     t_filenm       fnm[] = {
177         { efEDR, "-f", NULL, ffOPTRD },
178         { efDAT, "-groups", "groups", ffREAD },
179         { efDAT, "-eref",   "eref",   ffOPTRD },
180         { efXPM, "-emat",   "emat",   ffWRITE },
181         { efXVG, "-etot",   "energy", ffWRITE }
182     };
183 #define NFILE asize(fnm)
184
185     if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
186                            NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
187     {
188         return 0;
189     }
190
191     for (i = 0; (i < egNR+egSP); i++)
192     {
193         egrp_use[i] = FALSE;
194     }
195     egrp_use[egCOULSR] = bCoulSR;
196     egrp_use[egLJSR]   = bLJSR;
197     egrp_use[egBHAMSR] = bBhamSR;
198     egrp_use[egCOULLR] = bCoulLR;
199     egrp_use[egLJLR]   = bLJLR;
200     egrp_use[egBHAMLR] = bBhamLR;
201     egrp_use[egCOUL14] = bCoul14;
202     egrp_use[egLJ14]   = bLJ14;
203     egrp_use[egTotal]  = TRUE;
204
205     bRef = opt2bSet("-eref", NFILE, fnm);
206     in   = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
207     do_enxnms(in, &nre, &enm);
208
209     if (nre == 0)
210     {
211         gmx_fatal(FARGS, "No energies!\n");
212     }
213
214     bCutmax = opt2parg_bSet("-max", asize(pa), pa);
215     bCutmin = opt2parg_bSet("-min", asize(pa), pa);
216
217     nenergy = 0;
218
219     /* Read groupnames from input file and construct selection of
220        energy groups from it*/
221
222     fprintf(stderr, "Will read groupnames from inputfile\n");
223     ngroups = get_lines(opt2fn("-groups", NFILE, fnm), &groups);
224     fprintf(stderr, "Read %d groups\n", ngroups);
225     snew(set, static_cast<size_t>(sqr(ngroups)*egNR/2));
226     n     = 0;
227     prevk = 0;
228     for (i = 0; (i < ngroups); i++)
229     {
230         fprintf(stderr, "\rgroup %d", i);
231         for (j = i; (j < ngroups); j++)
232         {
233             for (m = 0; (m < egNR); m++)
234             {
235                 if (egrp_use[m])
236                 {
237                     sprintf(groupname, "%s:%s-%s", egrp_nm[m], groups[i], groups[j]);
238 #ifdef DEBUG
239                     fprintf(stderr, "\r%-15s %5d", groupname, n);
240 #endif
241                     for (k = prevk; (k < prevk+nre); k++)
242                     {
243                         if (std::strcmp(enm[k%nre].name, groupname) == 0)
244                         {
245                             set[n++] = k;
246                             break;
247                         }
248                     }
249                     if (k == prevk+nre)
250                     {
251                         fprintf(stderr, "WARNING! could not find group %s (%d,%d)"
252                                 "in energy file\n", groupname, i, j);
253                     }
254                     else
255                     {
256                         prevk = k;
257                     }
258                 }
259             }
260         }
261     }
262     fprintf(stderr, "\n");
263     nset = n;
264     snew(eneset, nset+1);
265     fprintf(stderr, "Will select half-matrix of energies with %d elements\n", n);
266
267     /* Start reading energy frames */
268     snew(fr, 1);
269     do
270     {
271         do
272         {
273             bCont = do_enx(in, fr);
274             if (bCont)
275             {
276                 timecheck = check_times(fr->t);
277             }
278         }
279         while (bCont && (timecheck < 0));
280
281         if (timecheck == 0)
282         {
283 #define DONTSKIP(cnt) (skip) ? ((cnt % skip) == 0) : TRUE
284
285             if (bCont)
286             {
287                 fprintf(stderr, "\rRead frame: %d, Time: %.3f", teller, fr->t);
288
289                 if ((nenergy % 1000) == 0)
290                 {
291                     srenew(time, nenergy+1000);
292                     for (i = 0; (i <= nset); i++)
293                     {
294                         srenew(eneset[i], nenergy+1000);
295                     }
296                 }
297                 time[nenergy] = fr->t;
298                 sum           = 0;
299                 for (i = 0; (i < nset); i++)
300                 {
301                     eneset[i][nenergy] = fr->ener[set[i]].e;
302                     sum               += fr->ener[set[i]].e;
303                 }
304                 if (bSum)
305                 {
306                     eneset[nset][nenergy] = sum;
307                 }
308                 nenergy++;
309             }
310             teller++;
311         }
312     }
313     while (bCont && (timecheck == 0));
314
315     fprintf(stderr, "\n");
316
317     fprintf(stderr, "Will build energy half-matrix of %d groups, %d elements, "
318             "over %d frames\n", ngroups, nset, nenergy);
319
320     snew(emat, egNR+egSP);
321     for (j = 0; (j < egNR+egSP); j++)
322     {
323         if (egrp_use[m])
324         {
325             snew(emat[j], ngroups);
326             for (i = 0; (i < ngroups); i++)
327             {
328                 snew(emat[j][i], ngroups);
329             }
330         }
331     }
332     snew(groupnr, ngroups);
333     for (i = 0; (i < ngroups); i++)
334     {
335         groupnr[i] = i+1;
336     }
337     rlo.r  = 1.0, rlo.g  = 0.0, rlo.b  = 0.0;
338     rmid.r = 1.0, rmid.g = 1.0, rmid.b = 1.0;
339     rhi.r  = 0.0, rhi.g  = 0.0, rhi.b  = 1.0;
340     if (bMeanEmtx)
341     {
342         snew(e, ngroups);
343         for (i = 0; (i < ngroups); i++)
344         {
345             snew(e[i], nenergy);
346         }
347         n = 0;
348         for (i = 0; (i < ngroups); i++)
349         {
350             for (j = i; (j < ngroups); j++)
351             {
352                 for (m = 0; (m < egNR); m++)
353                 {
354                     if (egrp_use[m])
355                     {
356                         for (k = 0; (k < nenergy); k++)
357                         {
358                             emat[m][i][j] += eneset[n][k];
359                             e[i][k]       += eneset[n][k]; /* *0.5; */
360                             e[j][k]       += eneset[n][k]; /* *0.5; */
361                         }
362                         n++;
363                         emat[egTotal][i][j] += emat[m][i][j];
364                         emat[m][i][j]       /= nenergy;
365                         emat[m][j][i]        = emat[m][i][j];
366                     }
367                 }
368                 emat[egTotal][i][j] /= nenergy;
369                 emat[egTotal][j][i]  = emat[egTotal][i][j];
370             }
371         }
372         if (bFree)
373         {
374             if (bRef)
375             {
376                 fprintf(stderr, "Will read reference energies from inputfile\n");
377                 neref = get_lines(opt2fn("-eref", NFILE, fnm), &ereflines);
378                 fprintf(stderr, "Read %d reference energies\n", neref);
379                 snew(eref, neref);
380                 snew(erefres, neref);
381                 for (i = 0; (i < neref); i++)
382                 {
383                     snew(erefres[i], 5);
384                     sscanf(ereflines[i], "%s %lf", erefres[i], &edum);
385                     eref[i] = edum;
386                 }
387             }
388             snew(eaver, ngroups);
389             for (i = 0; (i < ngroups); i++)
390             {
391                 for (k = 0; (k < nenergy); k++)
392                 {
393                     eaver[i] += e[i][k];
394                 }
395                 eaver[i] /= nenergy;
396             }
397             beta = 1.0/(BOLTZ*reftemp);
398             snew(efree, ngroups);
399             snew(edif, ngroups);
400             for (i = 0; (i < ngroups); i++)
401             {
402                 expE = 0;
403                 for (k = 0; (k < nenergy); k++)
404                 {
405                     expE += std::exp(beta*(e[i][k]-eaver[i]));
406                 }
407                 efree[i] = std::log(expE/nenergy)/beta + eaver[i];
408                 if (bRef)
409                 {
410                     n = search_str2(neref, erefres, groups[i]);
411                     if (n != -1)
412                     {
413                         edif[i] = efree[i]-eref[n];
414                     }
415                     else
416                     {
417                         edif[i] = efree[i];
418                         fprintf(stderr, "WARNING: group %s not found "
419                                 "in reference energies.\n", groups[i]);
420                     }
421                 }
422                 else
423                 {
424                     edif[i] = 0;
425                 }
426             }
427         }
428
429         emid             = 0.0; /*(emin+emax)*0.5;*/
430         egrp_nm[egTotal] = "total";
431         for (m = 0; (m < egNR+egSP); m++)
432         {
433             if (egrp_use[m])
434             {
435                 emin = 1e10;
436                 emax = -1e10;
437                 for (i = 0; (i < ngroups); i++)
438                 {
439                     for (j = i; (j < ngroups); j++)
440                     {
441                         if (emat[m][i][j] > emax)
442                         {
443                             emax = emat[m][i][j];
444                         }
445                         else if (emat[m][i][j] < emin)
446                         {
447                             emin = emat[m][i][j];
448                         }
449                     }
450                 }
451                 if (emax == emin)
452                 {
453                     fprintf(stderr, "Matrix of %s energy is uniform at %f "
454                             "(will not produce output).\n", egrp_nm[m], emax);
455                 }
456                 else
457                 {
458                     fprintf(stderr, "Matrix of %s energy ranges from %f to %f\n",
459                             egrp_nm[m], emin, emax);
460                     if ((bCutmax) || (emax > cutmax))
461                     {
462                         emax = cutmax;
463                     }
464                     if ((bCutmin) || (emin < cutmin))
465                     {
466                         emin = cutmin;
467                     }
468                     if ((emax == cutmax) || (emin == cutmin))
469                     {
470                         fprintf(stderr, "Energy range adjusted: %f to %f\n", emin, emax);
471                     }
472
473                     sprintf(fn, "%s%s", egrp_nm[m], ftp2fn(efXPM, NFILE, fnm));
474                     sprintf(label, "%s Interaction Energies", egrp_nm[m]);
475                     out = gmx_ffopen(fn, "w");
476                     if (emin >= emid)
477                     {
478                         write_xpm(out, 0, label, "Energy (kJ/mol)",
479                                   "Residue Index", "Residue Index",
480                                   ngroups, ngroups, groupnr, groupnr, emat[m],
481                                   emid, emax, rmid, rhi, &nlevels);
482                     }
483                     else if (emax <= emid)
484                     {
485                         write_xpm(out, 0, label, "Energy (kJ/mol)",
486                                   "Residue Index", "Residue Index",
487                                   ngroups, ngroups, groupnr, groupnr, emat[m],
488                                   emin, emid, rlo, rmid, &nlevels);
489                     }
490                     else
491                     {
492                         write_xpm3(out, 0, label, "Energy (kJ/mol)",
493                                    "Residue Index", "Residue Index",
494                                    ngroups, ngroups, groupnr, groupnr, emat[m],
495                                    emin, emid, emax, rlo, rmid, rhi, &nlevels);
496                     }
497                     gmx_ffclose(out);
498                 }
499             }
500         }
501         snew(etot, egNR+egSP);
502         for (m = 0; (m < egNR+egSP); m++)
503         {
504             snew(etot[m], ngroups);
505             for (i = 0; (i < ngroups); i++)
506             {
507                 for (j = 0; (j < ngroups); j++)
508                 {
509                     etot[m][i] += emat[m][i][j];
510                 }
511             }
512         }
513
514         out = xvgropen(ftp2fn(efXVG, NFILE, fnm), "Mean Energy", "Residue", "kJ/mol",
515                        oenv);
516         xvgr_legend(out, 0, NULL, oenv);
517         j = 0;
518         if (output_env_get_print_xvgr_codes(oenv))
519         {
520             char str1[STRLEN], str2[STRLEN];
521             if (output_env_get_xvg_format(oenv) == exvgXMGR)
522             {
523                 sprintf(str1, "@ legend string ");
524                 sprintf(str2, " ");
525             }
526             else
527             {
528                 sprintf(str1, "@ s");
529                 sprintf(str2, " legend ");
530             }
531
532             for (m = 0; (m < egNR+egSP); m++)
533             {
534                 if (egrp_use[m])
535                 {
536                     fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, egrp_nm[m]);
537                 }
538             }
539             if (bFree)
540             {
541                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Free");
542             }
543             if (bFree)
544             {
545                 fprintf(out, "%s%d%s \"%s\"\n", str1, j++, str2, "Diff");
546             }
547             fprintf(out, "@TYPE xy\n");
548             fprintf(out, "#%3s", "grp");
549
550             for (m = 0; (m < egNR+egSP); m++)
551             {
552                 if (egrp_use[m])
553                 {
554                     fprintf(out, " %9s", egrp_nm[m]);
555                 }
556             }
557             if (bFree)
558             {
559                 fprintf(out, " %9s", "Free");
560             }
561             if (bFree)
562             {
563                 fprintf(out, " %9s", "Diff");
564             }
565             fprintf(out, "\n");
566         }
567         for (i = 0; (i < ngroups); i++)
568         {
569             fprintf(out, "%3.0f", groupnr[i]);
570             for (m = 0; (m < egNR+egSP); m++)
571             {
572                 if (egrp_use[m])
573                 {
574                     fprintf(out, " %9.5g", etot[m][i]);
575                 }
576             }
577             if (bFree)
578             {
579                 fprintf(out, " %9.5g", efree[i]);
580             }
581             if (bRef)
582             {
583                 fprintf(out, " %9.5g", edif[i]);
584             }
585             fprintf(out, "\n");
586         }
587         xvgrclose(out);
588     }
589     else
590     {
591         fprintf(stderr, "While typing at your keyboard, suddenly...\n"
592                 "...nothing happens.\nWARNING: Not Implemented Yet\n");
593 /*
594     out=ftp2FILE(efMAT,NFILE,fnm,"w");
595     n=0;
596     emin=emax=0.0;
597     for (k=0; (k<nenergy); k++) {
598       for (i=0; (i<ngroups); i++)
599     for (j=i+1; (j<ngroups); j++)
600       emat[i][j]=eneset[n][k];
601       sprintf(label,"t=%.0f ps",time[k]);
602       write_matrix(out,ngroups,1,ngroups,groupnr,emat,label,emin,emax,nlevels);
603       n++;
604     }
605     gmx_ffclose(out);
606  */
607     }
608     close_enx(in);
609
610     return 0;
611 }