mdtools/md_code/Thermo.cpp

#include <cmath>

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


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;
  Atom** atoms;
  
  n_atoms = entry_plug->n_atoms;
  atoms = entry_plug->atoms;

  kinetic = 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());
    }
  }
  
  kinetic = kinetic * 0.5 / e_convert;

  return kinetic;
}

double Thermo::getPotential(){
  
  double potential;
  int el, nSRI;
  SRI** sris;

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

  potential = 0.0;

  potential += entry_plug->longRange->get_potential();;

  // std::cerr << "long range potential: " << potential << "\n";

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

  return potential;
}

double Thermo::getTotalE(){

  double total;

  total = this->getKinetic() + this->getPotential();
  return total;
}

double Thermo::getTemperature(){

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