src/mdlib/ebin.c

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