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"
#include "simError.h"

#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;
  Molecule* molecules;

  molecules = entry_plug->molecules;
  nSRI = entry_plug->n_SRI;

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

  for( el=0; el<entry_plug->n_mol; el++ ){    
    potential_local += molecules[el].getPotential();
  }

#ifdef IS_MPI
  /*
  std::cerr << "node " << worldRank << ": before LONG RANGE pot = " << entry_plug->lrPot 
            << "; pot_local = " << potential_local 
            << "; pot = " << potential << "\n";
  */
#endif

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

#ifdef IS_MPI
  /*
  std::cerr << "node " << worldRank << ": after pot = " << potential << "\n";
  */
#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_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(){
  // returns pressure in units amu*fs^-2*Ang^-1
  // routine derived via viral theorem description in:
  // Paci, E. and Marchi, M. J.Phys.Chem. 1996, 100, 4314-4322

  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 vbar;
  const double kb = 8.31451e-7; // kb in amu, angstroms, fs, etc.
  double av2;
  double kebar;
  int ndf, ndf_local; // number of degrees of freedom
  int ndfRaw, ndfRaw_local; // 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;
  
  // Raw degrees of freedom that we have to set
  ndfRaw_local = 3 * entry_plug->n_atoms + 3 * entry_plug->n_oriented;

  // Degrees of freedom that can contain kinetic energy
  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);
  MPI::COMM_WORLD.Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM);
#else
  ndfRaw = ndfRaw_local;
  ndf = ndf_local;
#endif
  ndf = ndf - 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.

  getCOMVel(vdrift);
  
  //  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 );
      }
    }   
  }
}

void Thermo::getCOMVel(double vdrift[3]){

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

  // 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;
  }
  
}

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