OOPSE/libmdtools/NVT.cpp

#include "Atom.hpp"
#include "SRI.hpp"
#include "AbstractClasses.hpp"
#include "SimInfo.hpp"
#include "ForceFields.hpp"
#include "Thermo.hpp"
#include "ReadWrite.hpp"
#include "Integrator.hpp"
#include "simError.h" 


// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697

NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff):
  Integrator( theInfo, the_ff )
{
  zeta = 0.0;
  have_tau_thermostat = 0;
  have_target_temp = 0;
  have_qmass = 0;
}

void NVT::moveA() {
  
  int i,j,k;
  int atomIndex, aMatIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double ke;
  double angle;


  ke = tStats->getKinetic() * eConvert;
  zeta += dt2 * ( (2.0 * ke  -  NkBT) / qmass );

  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    aMatIndex = i * 9;
    
    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*zeta);

    // position whole step    
    for( j=atomIndex; j<(atomIndex+3); j++ )
      pos[j] += dt * vel[j];

   
    if( atoms[i]->isDirectional() ){

      dAtom = (DirectionalAtom *)atoms[i];
          
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and propagate a half step

      ji[0] = dAtom->getJx();
      ji[1] = dAtom->getJy();
      ji[2] = dAtom->getJz();
      
      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*zeta);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*zeta);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*zeta);
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step
      
      // rotate about the x-axis      
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] ); 
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the z-axis
      angle = dt * ji[2] / dAtom->getIzz();
      this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / dAtom->getIyy();
      this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / dAtom->getIxx();
      this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
    
  }
}

void NVT::moveB( void ){
  int i,j,k;
  int atomIndex;
  DirectionalAtom* dAtom;
  double Tb[3];
  double ji[3];
  double ke;
  

  ke = tStats->getKinetic() * eConvert;
  zeta += dt2 * ( (2.0 * ke  -  NkBT) / qmass );
  
  for( i=0; i<nAtoms; i++ ){
    atomIndex = i * 3;
    
    // velocity half step
    for( j=atomIndex; j<(atomIndex+3); j++ )
      vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert - vel[j]*zeta);
    
    if( atoms[i]->isDirectional() ){
      
      dAtom = (DirectionalAtom *)atoms[i];
      
      // get and convert the torque to body frame
      
      Tb[0] = dAtom->getTx();
      Tb[1] = dAtom->getTy();
      Tb[2] = dAtom->getTz();
      
      dAtom->lab2Body( Tb );
      
      // get the angular momentum, and complete the angular momentum
      // half step
      
      ji[0] = dAtom->getJx();
      ji[1] = dAtom->getJy();
      ji[2] = dAtom->getJz();
      
      ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*zeta);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*zeta);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*zeta);
      
      dAtom->setJx( ji[0] );
      dAtom->setJy( ji[1] );
      dAtom->setJz( ji[2] );
    }
  }
}

int NVT::readyCheck() {
 
  // First check to see if we have a target temperature. 
  // Not having one is fatal. 
  
  if (!have_target_temp) {
    sprintf( painCave.errMsg,
             "NVT error: You can't use the NVT integrator without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }
     
  // Next check to see that we have a reasonable number of degrees of freedom
  // and then set NkBT if we do have it.   Unreasonable numbers of DOFs
  // are also fatal.

  if (info->ndf > 0) {
    NkBT = (double)info->ndf * kB * targetTemp;
  } else {
    sprintf( painCave.errMsg,
             "NVT error: We got a silly number of degrees of freedom!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }
    
  // We have our choice on setting qmass or tauThermostat.  One of them
  // must be set.

  if (!have_qmass) {
    if (have_tau_thermostat) {
      sprintf( painCave.errMsg,
               "NVT info: Setting qMass = %lf\n", tauThermostat * NkBT);
      this->setQmass(tauThermostat * NkBT);      
      painCave.isFatal = 0;
      simError();
    } else {
      sprintf( painCave.errMsg,
               "NVT error: If you use the constant temperature\n"
               "   integrator, you must set either tauThermostat or qMass.\n");
      painCave.isFatal = 1;
      simError();
      return -1;
    }
  }
  
  return 1;
}

Revision:	561
Committed:	Fri Jun 20 20:29:36 2003 UTC (21 years ago) by mmeineke
File size:	5157 byte(s)
Log Message:	Most of the integrator and NVT seem to be working now.
#	User	Rev	Content
1	gezelter	560	#include "Atom.hpp"
2			#include "SRI.hpp"
3			#include "AbstractClasses.hpp"
4			#include "SimInfo.hpp"
5			#include "ForceFields.hpp"
6			#include "Thermo.hpp"
7			#include "ReadWrite.hpp"
8			#include "Integrator.hpp"
9	mmeineke	561	#include "simError.h"
10
11
12	gezelter	560	// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
13
14	mmeineke	561	NVT::NVT ( SimInfo theInfo, ForceFields the_ff):
15			Integrator( theInfo, the_ff )
16			{
17	gezelter	560	zeta = 0.0;
18			have_tau_thermostat = 0;
19			have_target_temp = 0;
20			have_qmass = 0;
21			}
22
23			void NVT::moveA() {
24
25			int i,j,k;
26			int atomIndex, aMatIndex;
27			DirectionalAtom* dAtom;
28			double Tb[3];
29			double ji[3];
30	mmeineke	561	double ke;
31			double angle;
32	gezelter	560
33	mmeineke	561
34	gezelter	560	ke = tStats->getKinetic() * eConvert;
35			zeta += dt2 * ( (2.0 * ke - NkBT) / qmass );
36
37			for( i=0; i<nAtoms; i++ ){
38			atomIndex = i * 3;
39			aMatIndex = i * 9;
40
41			// velocity half step
42			for( j=atomIndex; j<(atomIndex+3); j++ )
43			vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())eConvert - vel[j]zeta);
44
45			// position whole step
46			for( j=atomIndex; j<(atomIndex+3); j++ )
47			pos[j] += dt * vel[j];
48
49
50			if( atoms[i]->isDirectional() ){
51
52			dAtom = (DirectionalAtom *)atoms[i];
53
54			// get and convert the torque to body frame
55
56			Tb[0] = dAtom->getTx();
57			Tb[1] = dAtom->getTy();
58			Tb[2] = dAtom->getTz();
59
60			dAtom->lab2Body( Tb );
61
62			// get the angular momentum, and propagate a half step
63
64			ji[0] = dAtom->getJx();
65			ji[1] = dAtom->getJy();
66			ji[2] = dAtom->getJz();
67
68			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*zeta);
69			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*zeta);
70			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*zeta);
71
72			// use the angular velocities to propagate the rotation matrix a
73			// full time step
74
75			// rotate about the x-axis
76			angle = dt2 * ji[0] / dAtom->getIxx();
77	mmeineke	561	this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
78	gezelter	560
79			// rotate about the y-axis
80			angle = dt2 * ji[1] / dAtom->getIyy();
81	mmeineke	561	this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
82	gezelter	560
83			// rotate about the z-axis
84			angle = dt * ji[2] / dAtom->getIzz();
85	mmeineke	561	this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
86	gezelter	560
87			// rotate about the y-axis
88			angle = dt2 * ji[1] / dAtom->getIyy();
89	mmeineke	561	this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
90	gezelter	560
91			// rotate about the x-axis
92			angle = dt2 * ji[0] / dAtom->getIxx();
93	mmeineke	561	this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
94	gezelter	560
95			dAtom->setJx( ji[0] );
96			dAtom->setJy( ji[1] );
97			dAtom->setJz( ji[2] );
98			}
99
100			}
101			}
102
103	mmeineke	561	void NVT::moveB( void ){
104	gezelter	560	int i,j,k;
105			int atomIndex;
106			DirectionalAtom* dAtom;
107			double Tb[3];
108			double ji[3];
109	mmeineke	561	double ke;
110
111	gezelter	560
112			ke = tStats->getKinetic() * eConvert;
113			zeta += dt2 * ( (2.0 * ke - NkBT) / qmass );
114
115			for( i=0; i<nAtoms; i++ ){
116			atomIndex = i * 3;
117
118			// velocity half step
119			for( j=atomIndex; j<(atomIndex+3); j++ )
120			vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())eConvert - vel[j]zeta);
121
122			if( atoms[i]->isDirectional() ){
123
124			dAtom = (DirectionalAtom *)atoms[i];
125
126			// get and convert the torque to body frame
127
128			Tb[0] = dAtom->getTx();
129			Tb[1] = dAtom->getTy();
130			Tb[2] = dAtom->getTz();
131
132			dAtom->lab2Body( Tb );
133
134			// get the angular momentum, and complete the angular momentum
135			// half step
136
137			ji[0] = dAtom->getJx();
138			ji[1] = dAtom->getJy();
139			ji[2] = dAtom->getJz();
140
141			ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*zeta);
142			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*zeta);
143			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*zeta);
144
145			dAtom->setJx( ji[0] );
146			dAtom->setJy( ji[1] );
147			dAtom->setJz( ji[2] );
148			}
149			}
150			}
151
152			int NVT::readyCheck() {
153	mmeineke	561
154	gezelter	560	// First check to see if we have a target temperature.
155			// Not having one is fatal.
156
157			if (!have_target_temp) {
158			sprintf( painCave.errMsg,
159			"NVT error: You can't use the NVT integrator without a targetTemp!\n"
160			);
161			painCave.isFatal = 1;
162			simError();
163			return -1;
164			}
165
166			// Next check to see that we have a reasonable number of degrees of freedom
167			// and then set NkBT if we do have it. Unreasonable numbers of DOFs
168			// are also fatal.
169
170	mmeineke	561	if (info->ndf > 0) {
171			NkBT = (double)info->ndf * kB * targetTemp;
172	gezelter	560	} else {
173			sprintf( painCave.errMsg,
174			"NVT error: We got a silly number of degrees of freedom!\n"
175			);
176			painCave.isFatal = 1;
177			simError();
178			return -1;
179			}
180
181			// We have our choice on setting qmass or tauThermostat. One of them
182			// must be set.
183
184			if (!have_qmass) {
185			if (have_tau_thermostat) {
186			sprintf( painCave.errMsg,
187	mmeineke	561	"NVT info: Setting qMass = %lf\n", tauThermostat * NkBT);
188	gezelter	560	this->setQmass(tauThermostat * NkBT);
189			painCave.isFatal = 0;
190			simError();
191			} else {
192			sprintf( painCave.errMsg,
193			"NVT error: If you use the constant temperature\n"
194			" integrator, you must set either tauThermostat or qMass.\n");
195			painCave.isFatal = 1;
196			simError();
197			return -1;
198			}
199			}
200
201			return 1;
202			}
203