src/gromacs/mdlib/ebin.c

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