mdtools/md_code/Thermo.cpp

#include <cmath>
#include <mpi++.h>

#include "Thermo.hpp"
#include "SRI.hpp"
#include "LRI.hpp"
#include "Integrator.hpp"


double Thermo::getKinetic(){

  const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
  double vx2, vy2, vz2;
  double kinetic, v_sqr;
  int kl;
  double jx2, jy2, jz2; // the square of the angular momentums

  DirectionalAtom *dAtom;

  int n_atoms;
  double kinetic_global;
  Atom** atoms;

  
  n_atoms = entry_plug->n_atoms;
  atoms = entry_plug->atoms;

  kinetic = 0.0;
  kinetic_global = 0.0;
  for( kl=0; kl < n_atoms; kl++ ){

    vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
    vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
    vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();

    v_sqr = vx2 + vy2 + vz2;
    kinetic += atoms[kl]->getMass() * v_sqr;

    if( atoms[kl]->isDirectional() ){
            
      dAtom = (DirectionalAtom *)atoms[kl];
      
      jx2 = dAtom->getJx() * dAtom->getJx();    
      jy2 = dAtom->getJy() * dAtom->getJy();
      jz2 = dAtom->getJz() * dAtom->getJz();
      
      kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy()) 
        + (jz2 / dAtom->getIzz());
    }
  }
#ifdef IS_MPI
  MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
  kinetic = kinetic_global;
#endif

  kinetic = kinetic * 0.5 / e_convert;

  return kinetic;
}

double Thermo::getPotential(){
  
  double potential;
  double potential_global;
  int el, nSRI;
  SRI** sris;

  sris = entry_plug->sr_interactions;
  nSRI = entry_plug->n_SRI;

  potential = 0.0;
  potential_global = 0.0;
  potential += entry_plug->longRange->get_potential();;

  // std::cerr << "long range potential: " << potential << "\n";
  for( el=0; el<nSRI; el++ ){
    
    potential += sris[el]->get_potential();
  }

  // Get total potential for entire system from MPI.
#ifdef IS_MPI
  MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
  potential = potential_global;
#endif

  return potential;
}

double Thermo::getTotalE(){

  double total;

  total = this->getKinetic() + this->getPotential();
  return total;
}

double Thermo::getTemperature(){

  const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
  double temperature;
  
  int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
    - entry_plug->n_constraints - 3;

  temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
  return temperature;
}

double Thermo::getPressure(){

//  const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
// const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
//  const double conv_A_m = 1.0E-10; //convert A -> m

  return 0.0;
}

void Thermo::velocitize() {
  
  double x,y;
  double vx, vy, vz;
  double jx, jy, jz;
  int i, vr, vd; // velocity randomizer loop counters
  double vdrift[3];
  double mtot = 0.0;
  double vbar;
  const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
  double av2;
  double kebar;
  int ndf; // number of degrees of freedom
  int ndfRaw; // the raw number of degrees of freedom
  int n_atoms;
  Atom** atoms;
  DirectionalAtom* dAtom;
  double temperature;
  int n_oriented;
  int n_constraints;

  atoms         = entry_plug->atoms;
  n_atoms       = entry_plug->n_atoms;
  temperature   = entry_plug->target_temp;
  n_oriented    = entry_plug->n_oriented;
  n_constraints = entry_plug->n_constraints;
  

  ndfRaw = 3 * n_atoms + 3 * n_oriented;
  ndf = ndfRaw - n_constraints - 3;
  kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
  
  for(vr = 0; vr < n_atoms; vr++){
    
    // uses equipartition theory to solve for vbar in angstrom/fs

    av2 = 2.0 * kebar / atoms[vr]->getMass();
    vbar = sqrt( av2 );

//     vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
    
    // picks random velocities from a gaussian distribution
    // centered on vbar
#ifndef USE_SPRNG
    /* If we are using mpi, we need to use the SPRNG random
       generator. The non drand48 generator will just repeat
       the same numbers for every node creating a non-gaussian
       distribution for the simulation. drand48 is fine for the
       single processor version of the code, but SPRNG should
       still be preferred for consistency.
    */
#ifdef IS_MPI
#error "SPRNG random number generator must be used for MPI"
#else
#warning "Using drand48 for random number generation"
#endif   
    x = drand48();
    y = drand48();
    vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
    
    x = drand48();
    y = drand48();
    vy = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
    
    x = drand48();
    y = drand48();
    vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
#endif

#ifdef USE_SPRNG
    vx = vbar * entry_plug->gaussStream->getGaussian();
    vy = vbar * entry_plug->gaussStream->getGaussian();
    vz = vbar * entry_plug->gaussStream->getGaussian();
#endif

    atoms[vr]->set_vx( vx ); 
    atoms[vr]->set_vy( vy );
    atoms[vr]->set_vz( vz );
  }
  
  //  Corrects for the center of mass drift.
  // sums all the momentum and divides by total mass.
  
  mtot = 0.0;
  vdrift[0] = 0.0;
  vdrift[1] = 0.0;
  vdrift[2] = 0.0;
  for(vd = 0; vd < n_atoms; vd++){
    
    vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
    vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
    vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
    
    mtot = mtot + atoms[vd]->getMass();
  }
  
  for (vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift[vd] / mtot;
  }
  
  for(vd = 0; vd < n_atoms; vd++){
    
    vx = atoms[vd]->get_vx();
    vy = atoms[vd]->get_vy();
    vz = atoms[vd]->get_vz();
    
    
    vx -= vdrift[0];
    vy -= vdrift[1];
    vz -= vdrift[2];
    
    atoms[vd]->set_vx(vx);
    atoms[vd]->set_vy(vy);
    atoms[vd]->set_vz(vz);
  }
  if( n_oriented ){
  
    for( i=0; i<n_atoms; i++ ){
      
      if( atoms[i]->isDirectional() ){
        
        dAtom = (DirectionalAtom *)atoms[i];
#ifdef IS_MPI
#error "SPRNG random number generator must be used for MPI"
#else
#warning "Using drand48 for random number generation"
#endif   
        
        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
        x = drand48();
        y = drand48();
        jx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
        
        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
        x = drand48();
        y = drand48();
        jy = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
        
        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
        x = drand48();
        y = drand48();
        jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
#endif
#ifdef USE_SPRNG
        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
        jx = vbar * entry_plug->gaussStream->getGaussian();

        vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
        jy = vbar * entry_plug->gaussStream->getGaussian();

        vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
        jz = vbar * entry_plug->gaussStream->getGaussian();
#endif
        
        dAtom->setJx( jx );
        dAtom->setJy( jy );
        dAtom->setJz( jz );
      }
    }   
  }
}
Revision:	221
Committed:	Thu Jan 2 20:14:08 2003 UTC (21 years, 6 months ago) by chuckv
File size:	6918 byte(s)
Log Message:	Thermo now can use SPRNG or rand48 (if not MPI). Finished up work on randomSPRNG. Still need to fix access to MPISimulation.
#	Content
1	#include <cmath>
2	#include <mpi++.h>
3
4	#include "Thermo.hpp"
5	#include "SRI.hpp"
6	#include "LRI.hpp"
7	#include "Integrator.hpp"
8
9
10	double Thermo::getKinetic(){
11
12	const double e_convert = 4.184E-4; // convert kcal/mol -> (amu A^2)/fs^2
13	double vx2, vy2, vz2;
14	double kinetic, v_sqr;
15	int kl;
16	double jx2, jy2, jz2; // the square of the angular momentums
17
18	DirectionalAtom *dAtom;
19
20	int n_atoms;
21	double kinetic_global;
22	Atom** atoms;
23
24
25	n_atoms = entry_plug->n_atoms;
26	atoms = entry_plug->atoms;
27
28	kinetic = 0.0;
29	kinetic_global = 0.0;
30	for( kl=0; kl < n_atoms; kl++ ){
31
32	vx2 = atoms[kl]->get_vx() * atoms[kl]->get_vx();
33	vy2 = atoms[kl]->get_vy() * atoms[kl]->get_vy();
34	vz2 = atoms[kl]->get_vz() * atoms[kl]->get_vz();
35
36	v_sqr = vx2 + vy2 + vz2;
37	kinetic += atoms[kl]->getMass() * v_sqr;
38
39	if( atoms[kl]->isDirectional() ){
40
41	dAtom = (DirectionalAtom *)atoms[kl];
42
43	jx2 = dAtom->getJx() * dAtom->getJx();
44	jy2 = dAtom->getJy() * dAtom->getJy();
45	jz2 = dAtom->getJz() * dAtom->getJz();
46
47	kinetic += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
48	+ (jz2 / dAtom->getIzz());
49	}
50	}
51	#ifdef IS_MPI
52	MPI_COMM_WORLD.Allreduce(&kinetic,&kinetic_global,1,MPI_DOUBLE,MPI_SUM);
53	kinetic = kinetic_global;
54	#endif
55
56	kinetic = kinetic * 0.5 / e_convert;
57
58	return kinetic;
59	}
60
61	double Thermo::getPotential(){
62
63	double potential;
64	double potential_global;
65	int el, nSRI;
66	SRI** sris;
67
68	sris = entry_plug->sr_interactions;
69	nSRI = entry_plug->n_SRI;
70
71	potential = 0.0;
72	potential_global = 0.0;
73	potential += entry_plug->longRange->get_potential();;
74
75	// std::cerr << "long range potential: " << potential << "\n";
76	for( el=0; el<nSRI; el++ ){
77
78	potential += sris[el]->get_potential();
79	}
80
81	// Get total potential for entire system from MPI.
82	#ifdef IS_MPI
83	MPI_COMM_WORLD.Allreduce(&potential,&potential_global,1,MPI_DOUBLE,MPI_SUM);
84	potential = potential_global;
85	#endif
86
87	return potential;
88	}
89
90	double Thermo::getTotalE(){
91
92	double total;
93
94	total = this->getKinetic() + this->getPotential();
95	return total;
96	}
97
98	double Thermo::getTemperature(){
99
100	const double kb = 1.9872179E-3; // boltzman's constant in kcal/(mol K)
101	double temperature;
102
103	int ndf = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
104	- entry_plug->n_constraints - 3;
105
106	temperature = ( 2.0 * this->getKinetic() ) / ( ndf * kb );
107	return temperature;
108	}
109
110	double Thermo::getPressure(){
111
112	// const double conv_Pa_atm = 9.901E-6; // convert Pa -> atm
113	// const double conv_internal_Pa = 1.661E-7; //convert amu/(fs^2 A) -> Pa
114	// const double conv_A_m = 1.0E-10; //convert A -> m
115
116	return 0.0;
117	}
118
119	void Thermo::velocitize() {
120
121	double x,y;
122	double vx, vy, vz;
123	double jx, jy, jz;
124	int i, vr, vd; // velocity randomizer loop counters
125	double vdrift[3];
126	double mtot = 0.0;
127	double vbar;
128	const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
129	double av2;
130	double kebar;
131	int ndf; // number of degrees of freedom
132	int ndfRaw; // the raw number of degrees of freedom
133	int n_atoms;
134	Atom** atoms;
135	DirectionalAtom* dAtom;
136	double temperature;
137	int n_oriented;
138	int n_constraints;
139
140	atoms = entry_plug->atoms;
141	n_atoms = entry_plug->n_atoms;
142	temperature = entry_plug->target_temp;
143	n_oriented = entry_plug->n_oriented;
144	n_constraints = entry_plug->n_constraints;
145
146
147	ndfRaw = 3 * n_atoms + 3 * n_oriented;
148	ndf = ndfRaw - n_constraints - 3;
149	kebar = kb * temperature * (double)ndf / ( 2.0 * (double)ndfRaw );
150
151	for(vr = 0; vr < n_atoms; vr++){
152
153	// uses equipartition theory to solve for vbar in angstrom/fs
154
155	av2 = 2.0 * kebar / atoms[vr]->getMass();
156	vbar = sqrt( av2 );
157
158	// vbar = sqrt( 8.31451e-7 * temperature / atoms[vr]->getMass() );
159
160	// picks random velocities from a gaussian distribution
161	// centered on vbar
162	#ifndef USE_SPRNG
163	/* If we are using mpi, we need to use the SPRNG random
164	generator. The non drand48 generator will just repeat
165	the same numbers for every node creating a non-gaussian
166	distribution for the simulation. drand48 is fine for the
167	single processor version of the code, but SPRNG should
168	still be preferred for consistency.
169	*/
170	#ifdef IS_MPI
171	#error "SPRNG random number generator must be used for MPI"
172	#else
173	#warning "Using drand48 for random number generation"
174	#endif
175	x = drand48();
176	y = drand48();
177	vx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
178
179	x = drand48();
180	y = drand48();
181	vy = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
182
183	x = drand48();
184	y = drand48();
185	vz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
186	#endif
187
188	#ifdef USE_SPRNG
189	vx = vbar * entry_plug->gaussStream->getGaussian();
190	vy = vbar * entry_plug->gaussStream->getGaussian();
191	vz = vbar * entry_plug->gaussStream->getGaussian();
192	#endif
193
194	atoms[vr]->set_vx( vx );
195	atoms[vr]->set_vy( vy );
196	atoms[vr]->set_vz( vz );
197	}
198
199	// Corrects for the center of mass drift.
200	// sums all the momentum and divides by total mass.
201
202	mtot = 0.0;
203	vdrift[0] = 0.0;
204	vdrift[1] = 0.0;
205	vdrift[2] = 0.0;
206	for(vd = 0; vd < n_atoms; vd++){
207
208	vdrift[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
209	vdrift[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
210	vdrift[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
211
212	mtot = mtot + atoms[vd]->getMass();
213	}
214
215	for (vd = 0; vd < 3; vd++) {
216	vdrift[vd] = vdrift[vd] / mtot;
217	}
218
219	for(vd = 0; vd < n_atoms; vd++){
220
221	vx = atoms[vd]->get_vx();
222	vy = atoms[vd]->get_vy();
223	vz = atoms[vd]->get_vz();
224
225
226	vx -= vdrift[0];
227	vy -= vdrift[1];
228	vz -= vdrift[2];
229
230	atoms[vd]->set_vx(vx);
231	atoms[vd]->set_vy(vy);
232	atoms[vd]->set_vz(vz);
233	}
234	if( n_oriented ){
235
236	for( i=0; i<n_atoms; i++ ){
237
238	if( atoms[i]->isDirectional() ){
239
240	dAtom = (DirectionalAtom *)atoms[i];
241	#ifdef IS_MPI
242	#error "SPRNG random number generator must be used for MPI"
243	#else
244	#warning "Using drand48 for random number generation"
245	#endif
246
247	vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
248	x = drand48();
249	y = drand48();
250	jx = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
251
252	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
253	x = drand48();
254	y = drand48();
255	jy = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
256
257	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
258	x = drand48();
259	y = drand48();
260	jz = vbar * sqrt( -2.0 * log(x)) * cos(2 * M_PI * y);
261	#endif
262	#ifdef USE_SPRNG
263	vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
264	jx = vbar * entry_plug->gaussStream->getGaussian();
265
266	vbar = sqrt( 2.0 * kebar * dAtom->getIyy() );
267	jy = vbar * entry_plug->gaussStream->getGaussian();
268
269	vbar = sqrt( 2.0 * kebar * dAtom->getIzz() );
270	jz = vbar * entry_plug->gaussStream->getGaussian();
271	#endif
272
273	dAtom->setJx( jx );
274	dAtom->setJy( jy );
275	dAtom->setJz( jz );
276	}
277	}
278	}
279	}