src/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  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 /* This file is completely threadsafe - keep it that way! */
  39 #ifdef HAVE_CONFIG_H
  40 #include <config.h>
  41 #endif
  42
  43 #include <math.h>
  44 #include <string.h>
  45 #include "sysstuff.h"
  46 #include "smalloc.h"
  47 #include "typedefs.h"
  48 #include "gmx_fatal.h"
  49 #include "string2.h"
  50 #include "ebin.h"
  51 #include "main.h"
  52 #include "maths.h"
  53 #include "vec.h"
  54 #include "physics.h"
  55
  56 t_ebin *mk_ebin(void)
  57 {
  58   t_ebin *eb;
  59
  60   snew(eb,1);
  61
  62   return eb;
  63 }
  64
  65 int get_ebin_space(t_ebin *eb,int nener,const char *enm[],const char *unit)
  66 {
  67     int  index;
  68     int  i,f;
  69     const char *u;
  70
  71     index = eb->nener;
  72     eb->nener += nener;
  73     srenew(eb->e,eb->nener);
  74     srenew(eb->e_sim,eb->nener);
  75     srenew(eb->enm,eb->nener);
  76     for(i=index; (i<eb->nener); i++)
  77     {
  78         eb->e[i].e        = 0;
  79         eb->e[i].eav      = 0;
  80         eb->e[i].esum     = 0;
  81         eb->e_sim[i].e    = 0;
  82         eb->e_sim[i].eav  = 0;
  83         eb->e_sim[i].esum = 0;
  84         eb->enm[i].name = strdup(enm[i-index]);
  85         if (unit != NULL)
  86         {
  87             eb->enm[i].unit = strdup(unit);
  88         }
  89         else
  90         {
  91             /* Determine the unit from the longname.
  92              * These units should have been defined in ifunc.c
  93              * But even better would be if all interactions functions
  94              * return energies and all non-interaction function
  95              * entries would be removed from the ifunc array.
  96              */
  97             u = unit_energy;
  98             for(f=0; f<F_NRE; f++)
  99             {
 100                 if (strcmp(eb->enm[i].name,
 101                            interaction_function[f].longname) == 0)
 102                 {
 103                     /* Only the terms in this list are not energies */
 104                     switch (f) {
 105                     case F_DISRESVIOL: u = unit_length;   break;
 106                     case F_ORIRESDEV:  u = "obs";         break;
 107                     case F_TEMP:       u = unit_temp_K;   break;
 108                     case F_VTEMP:       u = unit_temp_K;   break;
 109                     case F_PDISPCORR:
 110                     case F_PRES:       u = unit_pres_bar; break;
 111                     }
 112                 }
 113             }
 114             eb->enm[i].unit = strdup(u);
 115         }
 116     }
 117
 118     return index;
 119 }
 120
 121 void add_ebin(t_ebin *eb,int index,int nener,real ener[],gmx_bool bSum)
 122 {
 123     int      i,m;
 124     double   e,sum,sigma,invmm,diff;
 125     t_energy *eg,*egs;
 126
 127     if ((index+nener > eb->nener) || (index < 0))
 128     {
 129         gmx_fatal(FARGS,"%s-%d: Energies out of range: index=%d nener=%d maxener=%d",
 130                   __FILE__,__LINE__,index,nener,eb->nener);
 131     }
 132
 133     eg = &(eb->e[index]);
 134
 135     for(i=0; (i<nener); i++)
 136     {
 137         eg[i].e      = ener[i];
 138     }
 139
 140     if (bSum)
 141     {
 142         egs = &(eb->e_sim[index]);
 143
 144         m = eb->nsum;
 145
 146         if (m == 0)
 147         {
 148             for(i=0; (i<nener); i++)
 149             {
 150                 eg[i].eav    = 0;
 151                 eg[i].esum   = ener[i];
 152                 egs[i].esum += ener[i];
 153             }
 154         }
 155         else
 156         {
 157             invmm = (1.0/(double)m)/((double)m+1.0);
 158
 159             for(i=0; (i<nener); i++)
 160             {
 161                 /* Value for this component */
 162                 e = ener[i];
 163
 164                 /* first update sigma, then sum */
 165                 diff         = eg[i].esum - m*e;
 166                 eg[i].eav   += diff*diff*invmm;
 167                 eg[i].esum  += e;
 168                 egs[i].esum += e;
 169             }
 170         }
 171     }
 172 }
 173
 174 void ebin_increase_count(t_ebin *eb,gmx_bool bSum)
 175 {
 176     eb->nsteps++;
 177     eb->nsteps_sim++;
 178
 179     if (bSum)
 180     {
 181         eb->nsum++;
 182         eb->nsum_sim++;
 183     }
 184 }
 185
 186 void reset_ebin_sums(t_ebin *eb)
 187 {
 188     eb->nsteps = 0;
 189     eb->nsum   = 0;
 190     /* The actual sums are cleared when the next frame is stored */
 191 }
 192
 193 void pr_ebin(FILE *fp,t_ebin *eb,int index,int nener,int nperline,
 194              int prmode,gmx_bool bPrHead)
 195 {
 196     int  i,j,i0;
 197     real ee=0;
 198     int  rc;
 199     char buf[30];
 200
 201     rc = 0;
 202
 203     if (index < 0)
 204     {
 205         gmx_fatal(FARGS,"Invalid index in pr_ebin: %d",index);
 206     }
 207     if (nener == -1)
 208     {
 209         nener = eb->nener;
 210     }
 211     else
 212     {
 213         nener = index + nener;
 214     }
 215     for(i=index; (i<nener) && rc>=0; )
 216     {
 217         if (bPrHead)
 218         {
 219             i0=i;
 220             for(j=0; (j<nperline) && (i<nener) && rc>=0; j++,i++)
 221             {
 222                 if (strncmp(eb->enm[i].name,"Pres",4) == 0)
 223                 {
 224                     /* Print the pressure unit to avoid confusion */
 225                     sprintf(buf,"%s (%s)",eb->enm[i].name,unit_pres_bar);
 226                     rc = fprintf(fp,"%15s",buf);
 227                 }
 228                 else
 229                 {
 230                     rc = fprintf(fp,"%15s",eb->enm[i].name);
 231                 }
 232             }
 233
 234             if (rc >= 0)
 235             {
 236                 rc = fprintf(fp,"\n");
 237             }
 238
 239             i=i0;
 240         }
 241         for(j=0; (j<nperline) && (i<nener) && rc>=0; j++,i++)
 242         {
 243             switch (prmode) {
 244                 case eprNORMAL: ee = eb->e[i].e; break;
 245                 case eprAVER:   ee = eb->e_sim[i].esum/eb->nsum_sim; break;
 246                 default: gmx_fatal(FARGS,"Invalid print mode %d in pr_ebin",
 247                                    prmode);
 248             }
 249
 250             rc = fprintf(fp,"   %12.5e",ee);
 251         }
 252         if (rc >= 0)
 253         {
 254             rc = fprintf(fp,"\n");
 255         }
 256     }
 257     if (rc < 0)
 258     {
 259         gmx_fatal(FARGS,"Cannot write to logfile; maybe you are out of disk space?");
 260     }
 261 }
 262
 263 #ifdef DEBUGEBIN
 264 int main(int argc,char *argv[])
 265 {
 266 #define NE 12
 267 #define NT 7
 268 #define NS 5
 269
 270   t_ebin *eb;
 271   int    i;
 272   char   buf[25];
 273   char   *ce[NE],*ct[NT],*cs[NS];
 274   real   e[NE],t[NT],s[NS];
 275   int    ie,it,is;
 276
 277   eb=mk_ebin();
 278   for(i=0; (i<NE); i++) {
 279     e[i]=i;
 280     sprintf(buf,"e%d",i);
 281     ce[i]=strdup(buf);
 282   }
 283   ie=get_ebin_space(eb,NE,ce);
 284   add_ebin(eb,ie,NE,e,0);
 285   for(i=0; (i<NS); i++) {
 286     s[i]=i;
 287     sprintf(buf,"s%d",i);
 288     cs[i]=strdup(buf);
 289   }
 290   is=get_ebin_space(eb,NS,cs);
 291   add_ebin(eb,is,NS,s,0);
 292   for(i=0; (i<NT); i++) {
 293     t[i]=i;
 294     sprintf(buf,"t%d",i);
 295     ct[i]=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