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"

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

  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

    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 );
  }
  
  //  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 += 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];

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