mdtools/md_code/Thermo.cpp

#include <cmath>
#include <iostream>
using namespace std;

#ifdef IS_MPI
#include <mpi.h>
#include <mpi++.h>
#endif //is_mpi

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

#define BASE_SEED 123456789

Thermo::Thermo( SimInfo* the_entry_plug ) { 
  entry_plug = the_entry_plug;
  int baseSeed = BASE_SEED;
  
  gaussStream = new gaussianSPRNG( baseSeed );
}

Thermo::~Thermo(){
  delete gaussStream;
}

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 //is_mpi

  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->lrPot;

  // 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 // is_mpi

  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  // is_mpi 

    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 // use_spring

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

    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];

#ifndef USE_SPRNG

#ifdef IS_MPI
#error "SPRNG random number generator must be used for MPI"
#else  // is_mpi
        //warning "Using drand48 for random number generation"
#endif   // is_MPI
        
        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);

#else //use_sprng

        vbar = sqrt( 2.0 * kebar * dAtom->getIxx() );
        jx = vbar * gaussStream->getGaussian();

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

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