src/gmxlib/random.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "sysstuff.h"
  41 #include "smalloc.h"
  42 #include "physics.h"
  43 #include "typedefs.h"
  44 #include "vec.h"
  45 #include "gmx_random.h"
  46 #include "random.h"
  47 #include "mtop_util.h"
  48
  49 static void low_mspeed(real tempi,
  50                        gmx_mtop_t *mtop,rvec v[], gmx_rng_t rng)
  51 {
  52   int  i,m,nrdf;
  53   real boltz,sd;
  54   real ekin,temp,mass,scal;
  55   gmx_mtop_atomloop_all_t aloop;
  56   t_atom *atom;
  57
  58   boltz=BOLTZ*tempi;
  59   ekin=0.0;
  60   nrdf=0;
  61   aloop = gmx_mtop_atomloop_all_init(mtop);
  62   while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
  63     mass = atom->m;
  64     if (mass > 0) {
  65       sd=sqrt(boltz/mass);
  66       for (m=0; (m<DIM); m++) {
  67         v[i][m]=sd*gmx_rng_gaussian_real(rng);
  68         ekin += 0.5*mass*v[i][m]*v[i][m];
  69       }
  70       nrdf += DIM;
  71     }
  72   }
  73   temp=(2.0*ekin)/(nrdf*BOLTZ);
  74   if (temp > 0) {
  75     scal=sqrt(tempi/temp);
  76     for(i=0; (i<mtop->natoms); i++)
  77       for(m=0; (m<DIM); m++)
  78         v[i][m]*=scal;
  79   }
  80   fprintf(stderr,"Velocities were taken from a Maxwell distribution at %g K\n",
  81           tempi);
  82   if (debug) {
  83     fprintf(debug,
  84             "Velocities were taken from a Maxwell distribution\n"
  85             "Initial generated temperature: %12.5e (scaled to: %12.5e)\n",
  86             temp,tempi);
  87   }
  88 }
  89
  90 void maxwell_speed(real tempi,int seed,gmx_mtop_t *mtop, rvec v[])
  91 {
  92   atom_id *dummy;
  93   int     i;
  94   gmx_rng_t rng;
  95
  96   if (seed == -1) {
  97     seed = make_seed();
  98     fprintf(stderr,"Using random seed %d for generating velocities\n",seed);
  99   }
 100
 101   rng = gmx_rng_init(seed);
 102
 103   low_mspeed(tempi,mtop,v,rng);
 104
 105   gmx_rng_destroy(rng);
 106 }
 107
 108 real calc_cm(FILE *log,int natoms,real mass[],rvec x[],rvec v[],
 109              rvec xcm,rvec vcm,rvec acm,matrix L)
 110 {
 111   rvec dx,a0;
 112   real tm,m0;
 113   int  i,m;
 114
 115   clear_rvec(xcm);
 116   clear_rvec(vcm);
 117   clear_rvec(acm);
 118   tm=0.0;
 119   for(i=0; (i<natoms); i++) {
 120     m0=mass[i];
 121     tm+=m0;
 122     cprod(x[i],v[i],a0);
 123     for(m=0; (m<DIM); m++) {
 124       xcm[m]+=m0*x[i][m]; /* c.o.m. position */
 125       vcm[m]+=m0*v[i][m]; /* c.o.m. velocity */
 126       acm[m]+=m0*a0[m];   /* rotational velocity around c.o.m. */
 127     }
 128   }
 129   cprod(xcm,vcm,a0);
 130   for(m=0; (m<DIM); m++) {
 131     xcm[m]/=tm;
 132     vcm[m]/=tm;
 133     acm[m]-=a0[m]/tm;
 134   }
 135
 136 #define PVEC(str,v) fprintf(log,\
 137                             "%s[X]: %10.5e  %s[Y]: %10.5e  %s[Z]: %10.5e\n", \
 138                            str,v[0],str,v[1],str,v[2])
 139 #ifdef DEBUG
 140   PVEC("xcm",xcm);
 141   PVEC("acm",acm);
 142   PVEC("vcm",vcm);
 143 #endif
 144
 145   clear_mat(L);
 146   for(i=0; (i<natoms); i++) {
 147     m0=mass[i];
 148     for(m=0; (m<DIM); m++)
 149       dx[m]=x[i][m]-xcm[m];
 150     L[XX][XX]+=dx[XX]*dx[XX]*m0;
 151     L[XX][YY]+=dx[XX]*dx[YY]*m0;
 152     L[XX][ZZ]+=dx[XX]*dx[ZZ]*m0;
 153     L[YY][YY]+=dx[YY]*dx[YY]*m0;
 154     L[YY][ZZ]+=dx[YY]*dx[ZZ]*m0;
 155     L[ZZ][ZZ]+=dx[ZZ]*dx[ZZ]*m0;
 156   }
 157 #ifdef DEBUG
 158   PVEC("L-x",L[XX]);
 159   PVEC("L-y",L[YY]);
 160   PVEC("L-z",L[ZZ]);
 161 #endif
 162
 163   return tm;
 164 }
 165
 166 void stop_cm(FILE *log,int natoms,real mass[],rvec x[],rvec v[])
 167 {
 168   rvec   xcm,vcm,acm;
 169   tensor L;
 170   int    i,m;
 171
 172 #ifdef DEBUG
 173   fprintf(log,"stopping center of mass motion...\n");
 174 #endif
 175   (void)calc_cm(log,natoms,mass,x,v,xcm,vcm,acm,L);
 176
 177   /* Subtract center of mass velocity */
 178   for(i=0; (i<natoms); i++)
 179     for(m=0; (m<DIM); m++)
 180       v[i][m]-=vcm[m];
 181
 182 #ifdef DEBUG
 183   (void)calc_cm(log,natoms,mass,x,v,xcm,vcm,acm,L);
 184 #endif
 185 }