src/mdlib/ebin.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_ebin_c = "$Id$";
  30
  31 #include <math.h>
  32 #include <string.h>
  33 #include "sysstuff.h"
  34 #include "smalloc.h"
  35 #include "typedefs.h"
  36 #include "fatal.h"
  37 #include "string2.h"
  38 #include "ebin.h"
  39 #include "main.h"
  40 #include "maths.h"
  41 #include "vec.h"
  42
  43 static real rms_ener(t_energy *e,int nsteps)
  44 {
  45   double eav,rms;
  46
  47   eav=e->esum/nsteps;
  48   rms=sqrt(e->eav/nsteps);
  49
  50   if (eav != 0.0)
  51     if (fabs(rms/eav) < 1e-6)
  52       rms=0.0;
  53
  54   return rms;
  55 }
  56
  57 t_ebin *mk_ebin(void)
  58 {
  59   t_ebin *eb;
  60
  61   snew(eb,1);
  62
  63   return eb;
  64 }
  65
  66 int get_ebin_space(t_ebin *eb,int nener,char *enm[])
  67 {
  68   int index;
  69   int i;
  70
  71   index=eb->nener;
  72   eb->nener+=nener;
  73   srenew(eb->e,eb->nener);
  74   srenew(eb->enm,eb->nener);
  75   for(i=index; (i<eb->nener); i++) {
  76     eb->e[i].e=0;
  77     eb->e[i].eav=0;
  78     eb->e[i].esum=0;
  79     eb->e[i].e2sum=0;
  80     eb->enm[i]=strdup(enm[i-index]);
  81   }
  82   return index;
  83 }
  84
  85 void add_ebin(t_ebin *eb,int index,int nener,real ener[],int step)
  86 {
  87   int      i,m;
  88   double   e,sum,sigma,invmm;
  89   t_energy *eg;
  90
  91   if ((index+nener > eb->nener) || (index < 0))
  92     fatal_error(0,"%s-%d: Energies out of range: index=%d nener=%d maxener=%d",
  93                 __FILE__,__LINE__,index,nener,eb->nener);
  94
  95   m      = step;
  96   if (m > 0)
  97     invmm = (1.0/(double)m)/((double)m+1.0);
  98   else
  99     invmm = 0.0;
 100
 101   eg=&(eb->e[index]);
 102
 103   for(i=0; (i<nener); i++) {
 104     /* Value for this component */
 105     e      = ener[i];
 106
 107     /* first update sigma, then sum */
 108     eg[i].e    = e;
 109     eg[i].eav  += sqr(eg[i].esum - m*e)*invmm;;
 110     eg[i].esum += e;
 111   }
 112 }
 113
 114 void pr_ebin(FILE *fp,t_ebin *eb,int index,int nener,int nperline,
 115              int prmode,int tsteps,bool bPrHead)
 116 {
 117   int  i,j,i0;
 118   real ee=0;
 119
 120   if (index < 0)
 121     fatal_error(0,"Invalid index in pr_ebin: %d",index);
 122   if (nener == -1)
 123     nener=eb->nener;
 124   else
 125     nener=index+nener;
 126   for(i=index; (i<nener); ) {
 127     if (bPrHead) {
 128       i0=i;
 129       for(j=0; (j<nperline) && (i<nener); j++,i++)
 130         fprintf(fp,"%15s",eb->enm[i]);
 131       fprintf(fp,"\n");
 132       i=i0;
 133     }
 134     for(j=0; (j<nperline) && (i<nener); j++,i++) {
 135       if (prmode == eprNORMAL)
 136         ee=eb->e[i].e;
 137       else if (prmode == eprRMS)
 138         ee=rms_ener(&(eb->e[i]),tsteps);
 139       else if (prmode == eprAVER)
 140         ee=eb->e[i].esum/tsteps;
 141       else
 142         fatal_error(0,"Invalid print mode %d in pr_ebin",prmode);
 143
 144       fprintf(fp,"   %12.5e",ee);
 145     }
 146     fprintf(fp,"\n");
 147   }
 148 }
 149
 150 #ifdef DEBUGEBIN
 151 int main(int argc,char *argv[])
 152 {
 153 #define NE 12
 154 #define NT 7
 155 #define NS 5
 156
 157   t_ebin *eb;
 158   int    i;
 159   char   buf[25];
 160   char   *ce[NE],*ct[NT],*cs[NS];
 161   real   e[NE],t[NT],s[NS];
 162   int    ie,it,is;
 163
 164   eb=mk_ebin();
 165   for(i=0; (i<NE); i++) {
 166     e[i]=i;
 167     sprintf(buf,"e%d",i);
 168     ce[i]=strdup(buf);
 169   }
 170   ie=get_ebin_space(eb,NE,ce);
 171   add_ebin(eb,ie,NE,e,0);
 172   for(i=0; (i<NS); i++) {
 173     s[i]=i;
 174     sprintf(buf,"s%d",i);
 175     cs[i]=strdup(buf);
 176   }
 177   is=get_ebin_space(eb,NS,cs);
 178   add_ebin(eb,is,NS,s,0);
 179   for(i=0; (i<NT); i++) {
 180     t[i]=i;
 181     sprintf(buf,"t%d",i);
 182     ct[i]=strdup(buf);
 183   }
 184   it=get_ebin_space(eb,NT,ct);
 185   add_ebin(eb,it,NT,t,0);
 186
 187   printf("Normal:\n");
 188   pr_ebin(stdout,eb,0,-1,5,eprNORMAL,1);
 189
 190   printf("Average:\n");
 191   pr_ebin(stdout,eb,ie,NE,5,eprAVER,1);
 192   pr_ebin(stdout,eb,is,NS,3,eprAVER,1);
 193   pr_ebin(stdout,eb,it,NT,4,eprAVER,1);
 194
 195   printf("RMS:\n");
 196   pr_ebin(stdout,eb,0,-1,5,eprRMS,1);
 197 }
 198 #endif