src/gromacs/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     {
  64         mass = atom->m;
  65         if (mass > 0)
  66         {
  67             sd = sqrt(boltz/mass);
  68             for (m = 0; (m < DIM); m++)
  69             {
  70                 v[i][m] = sd*gmx_rng_gaussian_real(rng);
  71                 ekin   += 0.5*mass*v[i][m]*v[i][m];
  72             }
  73             nrdf += DIM;
  74         }
  75     }
  76     temp = (2.0*ekin)/(nrdf*BOLTZ);
  77     if (temp > 0)
  78     {
  79         scal = sqrt(tempi/temp);
  80         for (i = 0; (i < mtop->natoms); i++)
  81         {
  82             for (m = 0; (m < DIM); m++)
  83             {
  84                 v[i][m] *= scal;
  85             }
  86         }
  87     }
  88     fprintf(stderr, "Velocities were taken from a Maxwell distribution at %g K\n",
  89             tempi);
  90     if (debug)
  91     {
  92         fprintf(debug,
  93                 "Velocities were taken from a Maxwell distribution\n"
  94                 "Initial generated temperature: %12.5e (scaled to: %12.5e)\n",
  95                 temp, tempi);
  96     }
  97 }
  98
  99 void maxwell_speed(real tempi, int seed, gmx_mtop_t *mtop, rvec v[])
 100 {
 101     atom_id  *dummy;
 102     int       i;
 103     gmx_rng_t rng;
 104
 105     if (seed == -1)
 106     {
 107         seed = make_seed();
 108         fprintf(stderr, "Using random seed %d for generating velocities\n", seed);
 109     }
 110
 111     rng = gmx_rng_init(seed);
 112
 113     low_mspeed(tempi, mtop, v, rng);
 114
 115     gmx_rng_destroy(rng);
 116 }
 117
 118 real calc_cm(FILE *log, int natoms, real mass[], rvec x[], rvec v[],
 119              rvec xcm, rvec vcm, rvec acm, matrix L)
 120 {
 121     rvec dx, a0;
 122     real tm, m0;
 123     int  i, m;
 124
 125     clear_rvec(xcm);
 126     clear_rvec(vcm);
 127     clear_rvec(acm);
 128     tm = 0.0;
 129     for (i = 0; (i < natoms); i++)
 130     {
 131         m0  = mass[i];
 132         tm += m0;
 133         cprod(x[i], v[i], a0);
 134         for (m = 0; (m < DIM); m++)
 135         {
 136             xcm[m] += m0*x[i][m]; /* c.o.m. position */
 137             vcm[m] += m0*v[i][m]; /* c.o.m. velocity */
 138             acm[m] += m0*a0[m];   /* rotational velocity around c.o.m. */
 139         }
 140     }
 141     cprod(xcm, vcm, a0);
 142     for (m = 0; (m < DIM); m++)
 143     {
 144         xcm[m] /= tm;
 145         vcm[m] /= tm;
 146         acm[m] -= a0[m]/tm;
 147     }
 148
 149 #define PVEC(str, v) fprintf(log, \
 150                              "%s[X]: %10.5e  %s[Y]: %10.5e  %s[Z]: %10.5e\n", \
 151                              str, v[0], str, v[1], str, v[2])
 152 #ifdef DEBUG
 153     PVEC("xcm", xcm);
 154     PVEC("acm", acm);
 155     PVEC("vcm", vcm);
 156 #endif
 157
 158     clear_mat(L);
 159     for (i = 0; (i < natoms); i++)
 160     {
 161         m0 = mass[i];
 162         for (m = 0; (m < DIM); m++)
 163         {
 164             dx[m] = x[i][m]-xcm[m];
 165         }
 166         L[XX][XX] += dx[XX]*dx[XX]*m0;
 167         L[XX][YY] += dx[XX]*dx[YY]*m0;
 168         L[XX][ZZ] += dx[XX]*dx[ZZ]*m0;
 169         L[YY][YY] += dx[YY]*dx[YY]*m0;
 170         L[YY][ZZ] += dx[YY]*dx[ZZ]*m0;
 171         L[ZZ][ZZ] += dx[ZZ]*dx[ZZ]*m0;
 172     }
 173 #ifdef DEBUG
 174     PVEC("L-x", L[XX]);
 175     PVEC("L-y", L[YY]);
 176     PVEC("L-z", L[ZZ]);
 177 #endif
 178
 179     return tm;
 180 }
 181
 182 void stop_cm(FILE *log, int natoms, real mass[], rvec x[], rvec v[])
 183 {
 184     rvec   xcm, vcm, acm;
 185     tensor L;
 186     int    i, m;
 187
 188 #ifdef DEBUG
 189     fprintf(log, "stopping center of mass motion...\n");
 190 #endif
 191     (void)calc_cm(log, natoms, mass, x, v, xcm, vcm, acm, L);
 192
 193     /* Subtract center of mass velocity */
 194     for (i = 0; (i < natoms); i++)
 195     {
 196         for (m = 0; (m < DIM); m++)
 197         {
 198             v[i][m] -= vcm[m];
 199         }
 200     }
 201
 202 #ifdef DEBUG
 203     (void)calc_cm(log, natoms, mass, x, v, xcm, vcm, acm, L);
 204 #endif
 205 }