3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
45 #include "gmx_random.h"
47 #include "mtop_util.h"
49 static void low_mspeed(real tempi,
50 gmx_mtop_t *mtop,rvec v[], gmx_rng_t rng)
54 real ekin,temp,mass,scal;
55 gmx_mtop_atomloop_all_t aloop;
61 aloop = gmx_mtop_atomloop_all_init(mtop);
62 while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
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];
73 temp=(2.0*ekin)/(nrdf*BOLTZ);
75 scal=sqrt(tempi/temp);
76 for(i=0; (i<mtop->natoms); i++)
77 for(m=0; (m<DIM); m++)
80 fprintf(stderr,"Velocities were taken from a Maxwell distribution at %g K\n",
84 "Velocities were taken from a Maxwell distribution\n"
85 "Initial generated temperature: %12.5e (scaled to: %12.5e)\n",
90 void maxwell_speed(real tempi,int seed,gmx_mtop_t *mtop, rvec v[])
98 fprintf(stderr,"Using random seed %d for generating velocities\n",seed);
101 rng = gmx_rng_init(seed);
103 low_mspeed(tempi,mtop,v,rng);
105 gmx_rng_destroy(rng);
108 real calc_cm(FILE *log,int natoms,real mass[],rvec x[],rvec v[],
109 rvec xcm,rvec vcm,rvec acm,matrix L)
119 for(i=0; (i<natoms); i++) {
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. */
130 for(m=0; (m<DIM); m++) {
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])
146 for(i=0; (i<natoms); 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;
166 void stop_cm(FILE *log,int natoms,real mass[],rvec x[],rvec v[])
173 fprintf(log,"stopping center of mass motion...\n");
175 (void)calc_cm(log,natoms,mass,x,v,xcm,vcm,acm,L);
177 /* Subtract center of mass velocity */
178 for(i=0; (i<natoms); i++)
179 for(m=0; (m<DIM); m++)
183 (void)calc_cm(log,natoms,mass,x,v,xcm,vcm,acm,L);