OOPSE/libmdtools/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 "Integrator.hpp"

#ifdef IS_MPI
#define __C
#include "mpiSimulation.hpp"
#endif // is_mpi


#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_local;
  double potential;
  int el, nSRI;
  SRI** sris;

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

  potential_local = 0.0;
  potential_local += entry_plug->lrPot;

  for( el=0; el<nSRI; el++ ){    
    potential_local += sris[el]->get_potential();
  }

  // Get total potential for entire system from MPI.
#ifdef IS_MPI
  MPI::COMM_WORLD.Allreduce(&potential_local,&potential,1,MPI_DOUBLE,MPI_SUM);
#else
  potential = potential_local; 
#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_local, ndf;
  
  ndf_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented
    - entry_plug->n_constraints;

#ifdef IS_MPI
  MPI::COMM_WORLD.Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM);
#else
  ndf = ndf_local;
#endif

  ndf = ndf - 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;
  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

    vx = vbar * gaussStream->getGaussian();
    vy = vbar * gaussStream->getGaussian();
    vz = vbar * gaussStream->getGaussian();

    atoms[vr]->set_vx( vx ); 
    atoms[vr]->set_vy( vy );
    atoms[vr]->set_vz( vz );
  }

  // Get the Center of Mass drift velocity.

  vdrift = getCOMVel();
  
  //  Corrects for the center of mass drift.
  // sums all the momentum and divides by total mass.

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

        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();
        
        dAtom->setJx( jx );
        dAtom->setJy( jy );
        dAtom->setJz( jz );
      }
    }   
  }
}

double* Thermo::getCOMVel(){

  double mtot, mtot_local;
  double* vdrift;
  double vdrift_local[3];
  int vd, n_atoms;
  Atom** atoms;

  vdrift = new double[3];
  // We are very careless here with the distinction between n_atoms and n_local
  // We should really fix this before someone pokes an eye out.

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

  mtot_local = 0.0;
  vdrift_local[0] = 0.0;
  vdrift_local[1] = 0.0;
  vdrift_local[2] = 0.0;
  
  for(vd = 0; vd < n_atoms; vd++){
    
    vdrift_local[0] += atoms[vd]->get_vx() * atoms[vd]->getMass();
    vdrift_local[1] += atoms[vd]->get_vy() * atoms[vd]->getMass();
    vdrift_local[2] += atoms[vd]->get_vz() * atoms[vd]->getMass();
    
    mtot_local += atoms[vd]->getMass();
  }

#ifdef IS_MPI
  MPI::COMM_WORLD.Allreduce(&mtot_local,&mtot,1,MPI_DOUBLE,MPI_SUM);
  MPI::COMM_WORLD.Allreduce(&vdrift_local,&vdrift,3,MPI_DOUBLE,MPI_SUM);
#else
  mtot = mtot_local;
  for(vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift_local[vd];
  }
#endif
    
  for (vd = 0; vd < 3; vd++) {
    vdrift[vd] = vdrift[vd] / mtot;
  }
  
  return vdrift;
}

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