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

template<typename T> NVT<T>::NVT ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  chi = 0.0;
  have_tau_thermostat = 0;
  have_target_temp = 0;
}

template<typename T> void NVT<T>::moveA() {
  
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double A[3][3], I[3][3];
  double angle, mass;
  double vel[3], pos[3], frc[3];

  double instTemp;

  instTemp = tStats->getTemperature();

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

  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );
    atoms[i]->getPos( pos );
    atoms[i]->getFrc( frc );

    mass = atoms[i]->getMass();

    for (j=0; j < 3; j++) {
      // velocity half step
      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
      // position whole step
      pos[j] += dt * vel[j];
    }

    atoms[i]->setVel( vel );
    atoms[i]->setPos( pos );
   
    if( atoms[i]->isDirectional() ){

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

      dAtom->getJ( ji );

      for (j=0; j < 3; j++) 
        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
      
      // use the angular velocities to propagate the rotation matrix a
      // full time step

      dAtom->getA(A);
      dAtom->getI(I);
    
      // rotate about the x-axis      
      angle = dt2 * ji[0] / I[0][0];
      this->rotate( 1, 2, angle, ji, A ); 

      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate( 2, 0, angle, ji, A );
      
      // rotate about the z-axis
      angle = dt * ji[2] / I[2][2];
      this->rotate( 0, 1, angle, ji, A);
      
      // rotate about the y-axis
      angle = dt2 * ji[1] / I[1][1];
      this->rotate( 2, 0, angle, ji, A );
      
       // rotate about the x-axis
      angle = dt2 * ji[0] / I[0][0];
      this->rotate( 1, 2, angle, ji, A );
      
      dAtom->setJ( ji );
      dAtom->setA( A  );    
    }    
  }
}

template<typename T> void NVT<T>::moveB( void ){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  double instTemp;
  
  instTemp = tStats->getTemperature();
  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
  
  for( i=0; i<nAtoms; i++ ){

    atoms[i]->getVel( vel );
    atoms[i]->getFrc( frc );

    mass = atoms[i]->getMass();

    // velocity half step
    for (j=0; j < 3; j++) 
      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
    
    atoms[i]->setVel( vel );

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

      dAtom = (DirectionalAtom *)atoms[i];

      // get and convert the torque to body frame      

      dAtom->getTrq( Tb );
      dAtom->lab2Body( Tb );

      // get the angular momentum, and propagate a half step

      dAtom->getJ( ji );

      for (j=0; j < 3; j++) 
        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
      

      dAtom->setJ( ji );
    }
  }
}

template<typename T> int NVT<T>::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:	645
Committed:	Tue Jul 22 19:54:52 2003 UTC (20 years, 11 months ago) by tim
File size:	4071 byte(s)
Log Message:	* empty log message *
#	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	tim	645	template<typename T> NVT<T>::NVT ( SimInfo theInfo, ForceFields the_ff):
15			T( theInfo, the_ff )
16	mmeineke	561	{
17	gezelter	565	chi = 0.0;
18	gezelter	560	have_tau_thermostat = 0;
19			have_target_temp = 0;
20			}
21
22	tim	645	template<typename T> void NVT<T>::moveA() {
23	gezelter	560
24	gezelter	600	int i, j;
25	gezelter	560	DirectionalAtom* dAtom;
26	gezelter	600	double Tb[3], ji[3];
27			double A[3][3], I[3][3];
28			double angle, mass;
29			double vel[3], pos[3], frc[3];
30
31	gezelter	565	double instTemp;
32	gezelter	560
33	gezelter	565	instTemp = tStats->getTemperature();
34	mmeineke	561
35	gezelter	565	// first evolve chi a half step
36
37			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
38	gezelter	560
39			for( i=0; i<nAtoms; i++ ){
40
41	gezelter	600	atoms[i]->getVel( vel );
42			atoms[i]->getPos( pos );
43			atoms[i]->getFrc( frc );
44
45			mass = atoms[i]->getMass();
46
47			for (j=0; j < 3; j++) {
48			// velocity half step
49			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
50			// position whole step
51	gezelter	560	pos[j] += dt * vel[j];
52	gezelter	600	}
53	gezelter	560
54	gezelter	600	atoms[i]->setVel( vel );
55			atoms[i]->setPos( pos );
56	gezelter	560
57			if( atoms[i]->isDirectional() ){
58
59			dAtom = (DirectionalAtom *)atoms[i];
60
61			// get and convert the torque to body frame
62
63	gezelter	600	dAtom->getTrq( Tb );
64	gezelter	560	dAtom->lab2Body( Tb );
65
66			// get the angular momentum, and propagate a half step
67
68	gezelter	600	dAtom->getJ( ji );
69
70			for (j=0; j < 3; j++)
71			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
72	gezelter	560
73			// use the angular velocities to propagate the rotation matrix a
74			// full time step
75	gezelter	600
76			dAtom->getA(A);
77			dAtom->getI(I);
78
79	gezelter	560	// rotate about the x-axis
80	gezelter	600	angle = dt2 * ji[0] / I[0][0];
81			this->rotate( 1, 2, angle, ji, A );
82
83	gezelter	560	// rotate about the y-axis
84	gezelter	600	angle = dt2 * ji[1] / I[1][1];
85			this->rotate( 2, 0, angle, ji, A );
86	gezelter	560
87			// rotate about the z-axis
88	gezelter	600	angle = dt * ji[2] / I[2][2];
89			this->rotate( 0, 1, angle, ji, A);
90	gezelter	560
91			// rotate about the y-axis
92	gezelter	600	angle = dt2 * ji[1] / I[1][1];
93			this->rotate( 2, 0, angle, ji, A );
94	gezelter	560
95			// rotate about the x-axis
96	gezelter	600	angle = dt2 * ji[0] / I[0][0];
97			this->rotate( 1, 2, angle, ji, A );
98	gezelter	560
99	gezelter	600	dAtom->setJ( ji );
100			dAtom->setA( A );
101			}
102	gezelter	560	}
103			}
104
105	tim	645	template<typename T> void NVT<T>::moveB( void ){
106	gezelter	600	int i, j;
107	gezelter	560	DirectionalAtom* dAtom;
108	gezelter	600	double Tb[3], ji[3];
109			double vel[3], frc[3];
110			double mass;
111
112	gezelter	565	double instTemp;
113	mmeineke	561
114	gezelter	565	instTemp = tStats->getTemperature();
115			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
116	gezelter	560
117			for( i=0; i<nAtoms; i++ ){
118	gezelter	600
119			atoms[i]->getVel( vel );
120			atoms[i]->getFrc( frc );
121
122			mass = atoms[i]->getMass();
123
124	gezelter	560	// velocity half step
125	gezelter	600	for (j=0; j < 3; j++)
126			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
127	gezelter	560
128	gezelter	600	atoms[i]->setVel( vel );
129
130	gezelter	560	if( atoms[i]->isDirectional() ){
131	gezelter	600
132	gezelter	560	dAtom = (DirectionalAtom *)atoms[i];
133	gezelter	600
134			// get and convert the torque to body frame
135
136			dAtom->getTrq( Tb );
137	gezelter	560	dAtom->lab2Body( Tb );
138	gezelter	600
139			// get the angular momentum, and propagate a half step
140
141			dAtom->getJ( ji );
142
143			for (j=0; j < 3; j++)
144			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
145	gezelter	560
146	gezelter	600
147			dAtom->setJ( ji );
148	gezelter	560	}
149			}
150			}
151
152	tim	645	template<typename T> int NVT<T>::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			}