mdtools/md_code/Thermo.cpp

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

#ifdef IS_MPI
#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;
  
  cerr << "creating thermo stream\n";
  gaussStream = new gaussianSPRNG( baseSeed );
  cerr << "created thermo stream\n";
}

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