src/gromacs/mdlib/ebin.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) 2012,2014,2015,2017,2018, 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 /* This file is completely threadsafe - keep it that way! */
  38 #include "gmxpre.h"
  39
  40 #include "ebin.h"
  41
  42 #include <string.h>
  43
  44 #include <cmath>
  45
  46 #include "gromacs/math/units.h"
  47 #include "gromacs/math/utilities.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/topology/ifunc.h"
  50 #include "gromacs/trajectory/energyframe.h"
  51 #include "gromacs/utility/cstringutil.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 t_ebin *mk_ebin(void)
  56 {
  57     t_ebin *eb;
  58
  59     snew(eb, 1);
  60
  61     return eb;
  62 }
  63
  64 void done_ebin(t_ebin *eb)
  65 {
  66     for (int i = 0; i < eb->nener; i++)
  67     {
  68         sfree(eb->enm[i].name);
  69         sfree(eb->enm[i].unit);
  70     }
  71     sfree(eb->e);
  72     sfree(eb->e_sim);
  73     sfree(eb->enm);
  74 }
  75
  76 int get_ebin_space(t_ebin *eb, int nener, const char *enm[], const char *unit)
  77 {
  78     int         index;
  79     int         i, f;
  80     const char *u;
  81
  82     index      = eb->nener;
  83     eb->nener += nener;
  84     srenew(eb->e, eb->nener);
  85     srenew(eb->e_sim, eb->nener);
  86     srenew(eb->enm, eb->nener);
  87     for (i = index; (i < eb->nener); i++)
  88     {
  89         eb->e[i].e        = 0;
  90         eb->e[i].eav      = 0;
  91         eb->e[i].esum     = 0;
  92         eb->e_sim[i].e    = 0;
  93         eb->e_sim[i].eav  = 0;
  94         eb->e_sim[i].esum = 0;
  95         eb->enm[i].name   = gmx_strdup(enm[i-index]);
  96         if (unit != nullptr)
  97         {
  98             eb->enm[i].unit = gmx_strdup(unit);
  99         }
 100         else
 101         {
 102             /* Determine the unit from the longname.
 103              * These units should have been defined in ifunc.c
 104              * But even better would be if all interactions functions
 105              * return energies and all non-interaction function
 106              * entries would be removed from the ifunc array.
 107              */
 108             u = unit_energy;
 109             for (f = 0; f < F_NRE; f++)
 110             {
 111                 if (strcmp(eb->enm[i].name,
 112                            interaction_function[f].longname) == 0)
 113                 {
 114                     /* Only the terms in this list are not energies */
 115                     switch (f)
 116                     {
 117                         case F_DISRESVIOL: u = unit_length;   break;
 118                         case F_ORIRESDEV:  u = "obs";         break;
 119                         case F_TEMP:       u = unit_temp_K;   break;
 120                         case F_PDISPCORR:
 121                         case F_PRES:       u = unit_pres_bar; break;
 122                     }
 123                 }
 124             }
 125             eb->enm[i].unit = gmx_strdup(u);
 126         }
 127     }
 128
 129     return index;
 130 }
 131
 132 void add_ebin(t_ebin *eb, int index, int nener, const real ener[], gmx_bool bSum)
 133 {
 134     int       i, m;
 135     double    e, invmm, diff;
 136     t_energy *eg, *egs;
 137
 138     if ((index+nener > eb->nener) || (index < 0))
 139     {
 140         gmx_fatal(FARGS, "%s-%d: Energies out of range: index=%d nener=%d maxener=%d",
 141                   __FILE__, __LINE__, index, nener, eb->nener);
 142     }
 143
 144     eg = &(eb->e[index]);
 145
 146     for (i = 0; (i < nener); i++)
 147     {
 148         eg[i].e      = ener[i];
 149     }
 150
 151     if (bSum)
 152     {
 153         egs = &(eb->e_sim[index]);
 154
 155         m = eb->nsum;
 156
 157         if (m == 0)
 158         {
 159             for (i = 0; (i < nener); i++)
 160             {
 161                 eg[i].eav    = 0;
 162                 eg[i].esum   = ener[i];
 163                 egs[i].esum += ener[i];
 164             }
 165         }
 166         else
 167         {
 168             invmm = (1.0/(double)m)/((double)m+1.0);
 169
 170             for (i = 0; (i < nener); i++)
 171             {
 172                 /* Value for this component */
 173                 e = ener[i];
 174
 175                 /* first update sigma, then sum */
 176                 diff         = eg[i].esum - m*e;
 177                 eg[i].eav   += diff*diff*invmm;
 178                 eg[i].esum  += e;
 179                 egs[i].esum += e;
 180             }
 181         }
 182     }
 183 }
 184
 185 void ebin_increase_count(t_ebin *eb, gmx_bool bSum)
 186 {
 187     eb->nsteps++;
 188     eb->nsteps_sim++;
 189
 190     if (bSum)
 191     {
 192         eb->nsum++;
 193         eb->nsum_sim++;
 194     }
 195 }
 196
 197 void reset_ebin_sums(t_ebin *eb)
 198 {
 199     eb->nsteps = 0;
 200     eb->nsum   = 0;
 201     /* The actual sums are cleared when the next frame is stored */
 202 }
 203
 204 void pr_ebin(FILE *fp, t_ebin *eb, int index, int nener, int nperline,
 205              int prmode, gmx_bool bPrHead)
 206 {
 207     int  i, j, i0;
 208     int  rc;
 209     char buf[30];
 210
 211     rc = 0;
 212
 213     if (index < 0 || index > eb->nener)
 214     {
 215         gmx_fatal(FARGS, "Invalid index in pr_ebin: %d", index);
 216     }
 217     int start = index;
 218     if (nener > eb->nener)
 219     {
 220         gmx_fatal(FARGS, "Invalid nener in pr_ebin: %d", nener);
 221     }
 222     int end = eb->nener;
 223     if (nener != -1)
 224     {
 225         end = index + nener;
 226     }
 227     for (i = start; (i < end) && rc >= 0; )
 228     {
 229         if (bPrHead)
 230         {
 231             i0 = i;
 232             for (j = 0; (j < nperline) && (i < end) && rc >= 0; j++, i++)
 233             {
 234                 if (strncmp(eb->enm[i].name, "Pres", 4) == 0)
 235                 {
 236                     /* Print the pressure unit to avoid confusion */
 237                     sprintf(buf, "%s (%s)", eb->enm[i].name, unit_pres_bar);
 238                     rc = fprintf(fp, "%15s", buf);
 239                 }
 240                 else
 241                 {
 242                     rc = fprintf(fp, "%15s", eb->enm[i].name);
 243                 }
 244             }
 245
 246             if (rc >= 0)
 247             {
 248                 rc = fprintf(fp, "\n");
 249             }
 250
 251             i = i0;
 252         }
 253         for (j = 0; (j < nperline) && (i < end) && rc >= 0; j++, i++)
 254         {
 255             switch (prmode)
 256             {
 257                 case eprNORMAL:
 258                     rc = fprintf(fp, "   %12.5e", eb->e[i].e);
 259                     break;
 260                 case eprAVER:
 261                     if (eb->nsum_sim > 0)
 262                     {
 263                         rc = fprintf(fp, "   %12.5e", eb->e_sim[i].esum/eb->nsum_sim);
 264                     }
 265                     else
 266                     {
 267                         rc = fprintf(fp, "    %-12s", "N/A");
 268                     }
 269                     break;
 270                 default: gmx_fatal(FARGS, "Invalid print mode %d in pr_ebin",
 271                                    prmode);
 272             }
 273
 274         }
 275         if (rc >= 0)
 276         {
 277             rc = fprintf(fp, "\n");
 278         }
 279     }
 280     if (rc < 0)
 281     {
 282         gmx_fatal(FARGS, "Cannot write to logfile; maybe you are out of disk space?");
 283     }
 284 }