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;
  have_chi_tolerance = 0;
  integralOfChidt = 0.0;

  oldVel = new double[3*nAtoms];
  oldJi = new double[3*nAtoms];
}

template<typename T> NVT<T>::~NVT() {
  delete[] oldVel;
  delete[] oldJi;
}

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;

  // We need the temperature at time = t for the chi update below:

  instTemp = tStats->getTemperature();
  
  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  (use chi from previous step here):
      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  );    
    }    
  }
  
  if (nConstrained){
    constrainA();
  }

  // Finally, evolve chi a half step (just like a velocity) using 
  // temperature at time t, not time t+dt/2

  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
  integralOfChidt += chi*dt2;

}

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

  // Set things up for the iteration:

  oldChi = chi;

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

    atoms[i]->getVel( vel );

    for (j=0; j < 3; j++)
      oldVel[3*i + j]  = vel[j];

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

      dAtom = (DirectionalAtom *)atoms[i];

      dAtom->getJ( ji );

      for (j=0; j < 3; j++)
        oldJi[3*i + j] = ji[j];

    }
  }

  // do the iteration:

  for (k=0; k < 4; k++) {
    
    instTemp = tStats->getTemperature();

    // evolve chi another half step using the temperature at t + dt/2

    prevChi = chi;
    chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) / 
      (tauThermostat*tauThermostat);
  
    for( i=0; i<nAtoms; i++ ){

      atoms[i]->getFrc( frc );
      atoms[i]->getVel(vel);
      
      mass = atoms[i]->getMass();
      
      // velocity half step
      for (j=0; j < 3; j++) 
        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + 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 );       
             
        for (j=0; j < 3; j++) 
          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
      
        dAtom->setJ( ji );
      }
    }

    if (nConstrained){
      constrainB();
    }

    if (fabs(prevChi - chi) <= chiTolerance) break;
  }
  
  integralOfChidt += dt2*chi;
}

template<typename T> void NVT<T>::resetIntegrator( void ){
  
  chi = 0.0;
  integralOfChidt = 0.0;
}

template<typename T> int NVT<T>::readyCheck() {

  //check parent's readyCheck() first
  if (T::readyCheck() == -1)
    return -1;
   
  // 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;
  }    

  if (!have_chi_tolerance) {
    sprintf( painCave.errMsg,
             "NVT warning: setting chi tolerance to 1e-6\n");
    chiTolerance = 1e-6;
    have_chi_tolerance = 1;
    painCave.isFatal = 0;
    simError();
  }    

  return 1;    

}

template<typename T> double NVT<T>::getConservedQuantity(void){

  double conservedQuantity;
  double E_NVT;

  //HNVE
  conservedQuantity = tStats->getTotalE();
  //HNVE
  
  E_NVT =  (info->getNDF() * kB * targetTemp * 
                (integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0 )) / eConvert;

  conservedQuantity += E_NVT;

  //cerr << info->getTime() << "\t" << chi << "\t" << integralOfChidt << "\t" << E_NVT << endl;

  return conservedQuantity; 
}
Revision:	768
Committed:	Wed Sep 17 14:22:15 2003 UTC (20 years, 9 months ago) by mmeineke
File size:	6090 byte(s)
Log Message:	fixed NPTi to now work with constraints.
#	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	tim	763	have_chi_tolerance = 0;
21			integralOfChidt = 0.0;
22
23			oldVel = new double[3*nAtoms];
24			oldJi = new double[3*nAtoms];
25	gezelter	560	}
26
27	tim	763	template<typename T> NVT<T>::~NVT() {
28			delete[] oldVel;
29			delete[] oldJi;
30			}
31
32	tim	645	template<typename T> void NVT<T>::moveA() {
33	gezelter	560
34	gezelter	600	int i, j;
35	gezelter	560	DirectionalAtom* dAtom;
36	gezelter	600	double Tb[3], ji[3];
37			double A[3][3], I[3][3];
38			double angle, mass;
39			double vel[3], pos[3], frc[3];
40
41	gezelter	565	double instTemp;
42	gezelter	560
43	tim	763	// We need the temperature at time = t for the chi update below:
44
45	gezelter	565	instTemp = tStats->getTemperature();
46
47	gezelter	560	for( i=0; i<nAtoms; i++ ){
48
49	gezelter	600	atoms[i]->getVel( vel );
50			atoms[i]->getPos( pos );
51			atoms[i]->getFrc( frc );
52
53			mass = atoms[i]->getMass();
54
55			for (j=0; j < 3; j++) {
56	tim	763	// velocity half step (use chi from previous step here):
57	gezelter	600	vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
58			// position whole step
59	gezelter	560	pos[j] += dt * vel[j];
60	gezelter	600	}
61	gezelter	560
62	gezelter	600	atoms[i]->setVel( vel );
63			atoms[i]->setPos( pos );
64	gezelter	560
65			if( atoms[i]->isDirectional() ){
66
67			dAtom = (DirectionalAtom *)atoms[i];
68
69			// get and convert the torque to body frame
70
71	gezelter	600	dAtom->getTrq( Tb );
72	gezelter	560	dAtom->lab2Body( Tb );
73
74			// get the angular momentum, and propagate a half step
75
76	gezelter	600	dAtom->getJ( ji );
77
78			for (j=0; j < 3; j++)
79			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
80	gezelter	560
81			// use the angular velocities to propagate the rotation matrix a
82			// full time step
83	gezelter	600
84			dAtom->getA(A);
85			dAtom->getI(I);
86
87	gezelter	560	// rotate about the x-axis
88	gezelter	600	angle = dt2 * ji[0] / I[0][0];
89			this->rotate( 1, 2, angle, ji, A );
90
91	gezelter	560	// 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 z-axis
96	gezelter	600	angle = dt * ji[2] / I[2][2];
97			this->rotate( 0, 1, angle, ji, A);
98	gezelter	560
99			// rotate about the y-axis
100	gezelter	600	angle = dt2 * ji[1] / I[1][1];
101			this->rotate( 2, 0, angle, ji, A );
102	gezelter	560
103			// rotate about the x-axis
104	gezelter	600	angle = dt2 * ji[0] / I[0][0];
105			this->rotate( 1, 2, angle, ji, A );
106	gezelter	560
107	gezelter	600	dAtom->setJ( ji );
108			dAtom->setA( A );
109			}
110	gezelter	560	}
111	mmeineke	768
112			if (nConstrained){
113			constrainA();
114			}
115	tim	763
116			// Finally, evolve chi a half step (just like a velocity) using
117			// temperature at time t, not time t+dt/2
118
119			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
120			integralOfChidt += chi*dt2;
121
122	gezelter	560	}
123
124	tim	645	template<typename T> void NVT<T>::moveB( void ){
125	tim	763	int i, j, k;
126	gezelter	560	DirectionalAtom* dAtom;
127	gezelter	600	double Tb[3], ji[3];
128			double vel[3], frc[3];
129			double mass;
130	tim	763	double instTemp;
131			double oldChi, prevChi;
132	gezelter	600
133	tim	763	// Set things up for the iteration:
134
135			oldChi = chi;
136
137	gezelter	560	for( i=0; i<nAtoms; i++ ){
138	gezelter	600
139			atoms[i]->getVel( vel );
140
141	tim	763	for (j=0; j < 3; j++)
142			oldVel[3*i + j] = vel[j];
143	gezelter	600
144	gezelter	560	if( atoms[i]->isDirectional() ){
145	gezelter	600
146	gezelter	560	dAtom = (DirectionalAtom *)atoms[i];
147	gezelter	600
148	tim	763	dAtom->getJ( ji );
149	gezelter	600
150	tim	763	for (j=0; j < 3; j++)
151			oldJi[3*i + j] = ji[j];
152	gezelter	600
153	tim	763	}
154			}
155	gezelter	600
156	tim	763	// do the iteration:
157	gezelter	600
158	tim	763	for (k=0; k < 4; k++) {
159
160			instTemp = tStats->getTemperature();
161
162			// evolve chi another half step using the temperature at t + dt/2
163
164			prevChi = chi;
165			chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
166			(tauThermostat*tauThermostat);
167
168			for( i=0; i<nAtoms; i++ ){
169
170			atoms[i]->getFrc( frc );
171			atoms[i]->getVel(vel);
172
173			mass = atoms[i]->getMass();
174
175			// velocity half step
176	gezelter	600	for (j=0; j < 3; j++)
177	tim	763	vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j]chi);
178	gezelter	560
179	tim	763	atoms[i]->setVel( vel );
180
181			if( atoms[i]->isDirectional() ){
182
183			dAtom = (DirectionalAtom *)atoms[i];
184
185			// get and convert the torque to body frame
186
187			dAtom->getTrq( Tb );
188			dAtom->lab2Body( Tb );
189
190			for (j=0; j < 3; j++)
191			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
192
193			dAtom->setJ( ji );
194			}
195			}
196	gezelter	600
197	mmeineke	768	if (nConstrained){
198			constrainB();
199			}
200
201	tim	763	if (fabs(prevChi - chi) <= chiTolerance) break;
202	gezelter	560	}
203	tim	763
204			integralOfChidt += dt2*chi;
205	gezelter	560	}
206
207	mmeineke	746	template<typename T> void NVT<T>::resetIntegrator( void ){
208
209			chi = 0.0;
210	tim	763	integralOfChidt = 0.0;
211	mmeineke	746	}
212
213	tim	645	template<typename T> int NVT<T>::readyCheck() {
214	tim	658
215			//check parent's readyCheck() first
216			if (T::readyCheck() == -1)
217			return -1;
218
219	gezelter	560	// First check to see if we have a target temperature.
220			// Not having one is fatal.
221
222			if (!have_target_temp) {
223			sprintf( painCave.errMsg,
224			"NVT error: You can't use the NVT integrator without a targetTemp!\n"
225			);
226			painCave.isFatal = 1;
227			simError();
228			return -1;
229			}
230	gezelter	565
231			// We must set tauThermostat.
232
233			if (!have_tau_thermostat) {
234	gezelter	560	sprintf( painCave.errMsg,
235	gezelter	565	"NVT error: If you use the constant temperature\n"
236			" integrator, you must set tauThermostat.\n");
237	gezelter	560	painCave.isFatal = 1;
238			simError();
239			return -1;
240	gezelter	565	}
241	tim	763
242			if (!have_chi_tolerance) {
243			sprintf( painCave.errMsg,
244			"NVT warning: setting chi tolerance to 1e-6\n");
245			chiTolerance = 1e-6;
246			have_chi_tolerance = 1;
247			painCave.isFatal = 0;
248			simError();
249			}
250
251			return 1;
252
253	gezelter	560	}
254	tim	763
255			template<typename T> double NVT<T>::getConservedQuantity(void){
256
257			double conservedQuantity;
258			double E_NVT;
259
260			//HNVE
261			conservedQuantity = tStats->getTotalE();
262			//HNVE
263
264			E_NVT = (info->getNDF() * kB * targetTemp *
265			(integralOfChidt + tauThermostat * tauThermostat * chi * chi / 2.0 )) / eConvert;
266
267			conservedQuantity += E_NVT;
268
269			//cerr << info->getTime() << "\t" << chi << "\t" << integralOfChidt << "\t" << E_NVT << endl;
270
271			return conservedQuantity;
272			}