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 )
{
  chi = 0.0;
  have_tau_thermostat = 0;
  have_target_temp = 0;
}

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

  instTemp = tStats->getTemperature();

  // first evolve chi a half step
  
  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);

  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]*chi);

    // 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]*chi);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
      
      // 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 instTemp;
  
  instTemp = tStats->getTemperature();
  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
  
  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]*chi);
    
    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]*chi);
      ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
      ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
      
      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;
  }
  
  // We must set tauThermostat.
   
  if (!have_tau_thermostat) {
    sprintf( painCave.errMsg,
             "NVT error: If you use the constant temperature\n"
             "   integrator, you must set tauThermostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    
  return 1;
}

Revision:	565
Committed:	Tue Jun 24 22:51:57 2003 UTC (22 years, 4 months ago) by gezelter
File size:	4470 byte(s)
Log Message:	Fixes to NVT. Check them!
#	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	565	chi = 0.0;
18	gezelter	560	have_tau_thermostat = 0;
19			have_target_temp = 0;
20			}
21
22			void NVT::moveA() {
23
24			int i,j,k;
25			int atomIndex, aMatIndex;
26			DirectionalAtom* dAtom;
27			double Tb[3];
28			double ji[3];
29	gezelter	565	double instTemp;
30	mmeineke	561	double angle;
31	gezelter	560
32	gezelter	565	instTemp = tStats->getTemperature();
33	mmeineke	561
34	gezelter	565	// first evolve chi a half step
35
36			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
37	gezelter	560
38			for( i=0; i<nAtoms; i++ ){
39			atomIndex = i * 3;
40			aMatIndex = i * 9;
41
42			// velocity half step
43			for( j=atomIndex; j<(atomIndex+3); j++ )
44	gezelter	565	vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())eConvert - vel[j]chi);
45	gezelter	560
46			// position whole step
47			for( j=atomIndex; j<(atomIndex+3); j++ )
48			pos[j] += dt * vel[j];
49
50
51			if( atoms[i]->isDirectional() ){
52
53			dAtom = (DirectionalAtom *)atoms[i];
54
55			// get and convert the torque to body frame
56
57			Tb[0] = dAtom->getTx();
58			Tb[1] = dAtom->getTy();
59			Tb[2] = dAtom->getTz();
60
61			dAtom->lab2Body( Tb );
62
63			// get the angular momentum, and propagate a half step
64
65			ji[0] = dAtom->getJx();
66			ji[1] = dAtom->getJy();
67			ji[2] = dAtom->getJz();
68
69	gezelter	565	ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
70			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
71			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
72	gezelter	560
73			// use the angular velocities to propagate the rotation matrix a
74			// full time step
75
76			// rotate about the x-axis
77			angle = dt2 * ji[0] / dAtom->getIxx();
78	mmeineke	561	this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
79	gezelter	560
80			// rotate about the y-axis
81			angle = dt2 * ji[1] / dAtom->getIyy();
82	mmeineke	561	this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
83	gezelter	560
84			// rotate about the z-axis
85			angle = dt * ji[2] / dAtom->getIzz();
86	mmeineke	561	this->rotate( 0, 1, angle, ji, &Amat[aMatIndex] );
87	gezelter	560
88			// rotate about the y-axis
89			angle = dt2 * ji[1] / dAtom->getIyy();
90	mmeineke	561	this->rotate( 2, 0, angle, ji, &Amat[aMatIndex] );
91	gezelter	560
92			// rotate about the x-axis
93			angle = dt2 * ji[0] / dAtom->getIxx();
94	mmeineke	561	this->rotate( 1, 2, angle, ji, &Amat[aMatIndex] );
95	gezelter	560
96			dAtom->setJx( ji[0] );
97			dAtom->setJy( ji[1] );
98			dAtom->setJz( ji[2] );
99			}
100
101			}
102			}
103
104	mmeineke	561	void NVT::moveB( void ){
105	gezelter	560	int i,j,k;
106			int atomIndex;
107			DirectionalAtom* dAtom;
108			double Tb[3];
109			double ji[3];
110	gezelter	565	double instTemp;
111	mmeineke	561
112	gezelter	565	instTemp = tStats->getTemperature();
113			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
114	gezelter	560
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	gezelter	565	vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())eConvert - vel[j]chi);
121	gezelter	560
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	gezelter	565	ji[0] += dt2 * (Tb[0] * eConvert - ji[0]*chi);
142			ji[1] += dt2 * (Tb[1] * eConvert - ji[1]*chi);
143			ji[2] += dt2 * (Tb[2] * eConvert - ji[2]*chi);
144	gezelter	560
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	gezelter	565
166			// We must set tauThermostat.
167
168			if (!have_tau_thermostat) {
169	gezelter	560	sprintf( painCave.errMsg,
170	gezelter	565	"NVT error: If you use the constant temperature\n"
171			" integrator, you must set tauThermostat.\n");
172	gezelter	560	painCave.isFatal = 1;
173			simError();
174			return -1;
175	gezelter	565	}
176	gezelter	560	return 1;
177			}
178