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,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.
  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 <cmath>
  43 #include <cstring>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/math/units.h"
  48 #include "gromacs/math/utilities.h"
  49 #include "gromacs/math/vec.h"
  50 #include "gromacs/topology/ifunc.h"
  51 #include "gromacs/trajectory/energyframe.h"
  52 #include "gromacs/utility/basedefinitions.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/smalloc.h"
  56 #include "gromacs/utility/stringutil.h"
  57
  58 t_ebin* mk_ebin()
  59 {
  60     t_ebin* eb;
  61
  62     snew(eb, 1);
  63
  64     return eb;
  65 }
  66
  67 void done_ebin(t_ebin* eb)
  68 {
  69     for (int i = 0; i < eb->nener; i++)
  70     {
  71         sfree(eb->enm[i].name);
  72         sfree(eb->enm[i].unit);
  73     }
  74     sfree(eb->e);
  75     sfree(eb->e_sim);
  76     sfree(eb->enm);
  77     sfree(eb);
  78 }
  79
  80 int get_ebin_space(t_ebin* eb, int nener, const char* const enm[], const char* unit)
  81 {
  82     int         index;
  83     int         i, f;
  84     const char* u;
  85
  86     index = eb->nener;
  87     eb->nener += nener;
  88     srenew(eb->e, eb->nener);
  89     srenew(eb->e_sim, eb->nener);
  90     srenew(eb->enm, eb->nener);
  91     for (i = index; (i < eb->nener); i++)
  92     {
  93         eb->e[i].e        = 0;
  94         eb->e[i].eav      = 0;
  95         eb->e[i].esum     = 0;
  96         eb->e_sim[i].e    = 0;
  97         eb->e_sim[i].eav  = 0;
  98         eb->e_sim[i].esum = 0;
  99         eb->enm[i].name   = gmx_strdup(enm[i - index]);
 100         if (unit != nullptr)
 101         {
 102             eb->enm[i].unit = gmx_strdup(unit);
 103         }
 104         else
 105         {
 106             /* Determine the unit from the longname.
 107              * These units should have been defined in ifunc.c
 108              * But even better would be if all interactions functions
 109              * return energies and all non-interaction function
 110              * entries would be removed from the ifunc array.
 111              */
 112             u = unit_energy;
 113             for (f = 0; f < F_NRE; f++)
 114             {
 115                 if (strcmp(eb->enm[i].name, interaction_function[f].longname) == 0)
 116                 {
 117                     /* Only the terms in this list are not energies */
 118                     switch (f)
 119                     {
 120                         case F_DISRESVIOL: u = unit_length; break;
 121                         case F_ORIRESDEV: u = "obs"; break;
 122                         case F_TEMP: u = unit_temp_K; break;
 123                         case F_PDISPCORR:
 124                         case F_PRES: u = unit_pres_bar; break;
 125                     }
 126                 }
 127             }
 128             eb->enm[i].unit = gmx_strdup(u);
 129         }
 130     }
 131
 132     return index;
 133 }
 134
 135 // ICC 19 -O3 -msse2 generates wrong code. Lower optimization levels
 136 // and other SIMD levels seem fine, however.
 137 #if defined __ICC
 138 #    pragma intel optimization_level 2
 139 #endif
 140 void add_ebin(t_ebin* eb, int entryIndex, int nener, const real ener[], gmx_bool bSum)
 141 {
 142     int       i, m;
 143     double    e, invmm, diff;
 144     t_energy *eg, *egs;
 145
 146     if ((entryIndex + nener > eb->nener) || (entryIndex < 0))
 147     {
 148         gmx_fatal(FARGS, "%s-%d: Energies out of range: entryIndex=%d nener=%d maxener=%d",
 149                   __FILE__, __LINE__, entryIndex, nener, eb->nener);
 150     }
 151
 152     eg = &(eb->e[entryIndex]);
 153
 154     for (i = 0; (i < nener); i++)
 155     {
 156         eg[i].e = ener[i];
 157     }
 158
 159     if (bSum)
 160     {
 161         egs = &(eb->e_sim[entryIndex]);
 162
 163         m = eb->nsum;
 164
 165         if (m == 0)
 166         {
 167             for (i = 0; (i < nener); i++)
 168             {
 169                 eg[i].eav  = 0;
 170                 eg[i].esum = ener[i];
 171                 egs[i].esum += ener[i];
 172             }
 173         }
 174         else
 175         {
 176             invmm = (1.0 / m) / (m + 1.0);
 177
 178             for (i = 0; (i < nener); i++)
 179             {
 180                 /* Value for this component */
 181                 e = ener[i];
 182
 183                 /* first update sigma, then sum */
 184                 diff = eg[i].esum - m * e;
 185                 eg[i].eav += diff * diff * invmm;
 186                 eg[i].esum += e;
 187                 egs[i].esum += e;
 188             }
 189         }
 190     }
 191 }
 192
 193 // TODO It would be faster if this function was templated on both bSum
 194 // and whether eb->nsum was zero, to lift the branches out of the loop
 195 // over all possible energy terms, but that is true for a lot of the
 196 // ebin and mdebin functionality, so we should do it all or nothing.
 197 void add_ebin_indexed(t_ebin*                   eb,
 198                       int                       entryIndex,
 199                       gmx::ArrayRef<bool>       shouldUse,
 200                       gmx::ArrayRef<const real> ener,
 201                       gmx_bool                  bSum)
 202 {
 203
 204     GMX_ASSERT(shouldUse.size() == ener.size(), "View sizes must match");
 205     GMX_ASSERT(entryIndex + std::count(shouldUse.begin(), shouldUse.end(), true) <= eb->nener,
 206                gmx::formatString("Energies out of range: entryIndex=%d nener=%td maxener=%d", entryIndex,
 207                                  std::count(shouldUse.begin(), shouldUse.end(), true), eb->nener)
 208                        .c_str());
 209     GMX_ASSERT(entryIndex >= 0, "Must have non-negative entry");
 210
 211     const int    m              = eb->nsum;
 212     const double invmm          = (m == 0) ? 0 : (1.0 / m) / (m + 1.0);
 213     t_energy*    energyEntry    = &(eb->e[entryIndex]);
 214     t_energy*    simEnergyEntry = &(eb->e_sim[entryIndex]);
 215     auto         shouldUseIter  = shouldUse.begin();
 216     for (const auto& theEnergy : ener)
 217     {
 218         if (*shouldUseIter)
 219         {
 220             energyEntry->e = theEnergy;
 221             if (bSum)
 222             {
 223                 if (m == 0)
 224                 {
 225                     energyEntry->eav  = 0;
 226                     energyEntry->esum = theEnergy;
 227                     simEnergyEntry->esum += theEnergy;
 228                 }
 229                 else
 230                 {
 231                     /* first update sigma, then sum */
 232                     double diff = energyEntry->esum - m * theEnergy;
 233                     energyEntry->eav += diff * diff * invmm;
 234                     energyEntry->esum += theEnergy;
 235                     simEnergyEntry->esum += theEnergy;
 236                 }
 237                 ++simEnergyEntry;
 238             }
 239             ++energyEntry;
 240         }
 241         ++shouldUseIter;
 242     }
 243 }
 244
 245 void ebin_increase_count(int increment, t_ebin* eb, gmx_bool bSum)
 246 {
 247     eb->nsteps += increment;
 248     eb->nsteps_sim += increment;
 249
 250     if (bSum)
 251     {
 252         eb->nsum += increment;
 253         eb->nsum_sim += increment;
 254     }
 255 }
 256
 257 void reset_ebin_sums(t_ebin* eb)
 258 {
 259     eb->nsteps = 0;
 260     eb->nsum   = 0;
 261     /* The actual sums are cleared when the next frame is stored */
 262 }
 263
 264 void pr_ebin(FILE* fp, t_ebin* eb, int entryIndex, int nener, int nperline, int prmode, gmx_bool bPrHead)
 265 {
 266     int  i, j, i0;
 267     int  rc;
 268     char buf[30];
 269
 270     rc = 0;
 271
 272     if (entryIndex < 0 || entryIndex > eb->nener)
 273     {
 274         gmx_fatal(FARGS, "Invalid entryIndex in pr_ebin: %d", entryIndex);
 275     }
 276     int start = entryIndex;
 277     if (nener > eb->nener)
 278     {
 279         gmx_fatal(FARGS, "Invalid nener in pr_ebin: %d", nener);
 280     }
 281     int end = eb->nener;
 282     if (nener != -1)
 283     {
 284         end = entryIndex + nener;
 285     }
 286     for (i = start; (i < end) && rc >= 0;)
 287     {
 288         if (bPrHead)
 289         {
 290             i0 = i;
 291             for (j = 0; (j < nperline) && (i < end) && rc >= 0; j++, i++)
 292             {
 293                 if (strncmp(eb->enm[i].name, "Pres", 4) == 0)
 294                 {
 295                     /* Print the pressure unit to avoid confusion */
 296                     sprintf(buf, "%s (%s)", eb->enm[i].name, unit_pres_bar);
 297                     rc = fprintf(fp, "%15s", buf);
 298                 }
 299                 else
 300                 {
 301                     rc = fprintf(fp, "%15s", eb->enm[i].name);
 302                 }
 303             }
 304
 305             if (rc >= 0)
 306             {
 307                 rc = fprintf(fp, "\n");
 308             }
 309
 310             i = i0;
 311         }
 312         for (j = 0; (j < nperline) && (i < end) && rc >= 0; j++, i++)
 313         {
 314             switch (prmode)
 315             {
 316                 case eprNORMAL: rc = fprintf(fp, "   %12.5e", eb->e[i].e); break;
 317                 case eprAVER:
 318                     if (eb->nsum_sim > 0)
 319                     {
 320                         rc = fprintf(fp, "   %12.5e", eb->e_sim[i].esum / eb->nsum_sim);
 321                     }
 322                     else
 323                     {
 324                         rc = fprintf(fp, "    %-12s", "N/A");
 325                     }
 326                     break;
 327                 default: gmx_fatal(FARGS, "Invalid print mode %d in pr_ebin", prmode);
 328             }
 329         }
 330         if (rc >= 0)
 331         {
 332             rc = fprintf(fp, "\n");
 333         }
 334     }
 335     if (rc < 0)
 336     {
 337         gmx_fatal(FARGS, "Cannot write to logfile; maybe you are out of disk space?");
 338     }
 339 }