src/integrators/NVT.cpp

#include <math.h>

#include "primitives/Atom.hpp"
#include "primitives/SRI.hpp"
#include "primitives/AbstractClasses.hpp"
#include "brains/SimInfo.hpp"
#include "UseTheForce/ForceFields.hpp"
#include "brains/Thermo.hpp"
#include "io/ReadWrite.hpp"
#include "integrators/Integrator.hpp"
#include "utils/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 )
{
  GenericData* data;
  DoubleGenericData * chiValue;
  DoubleGenericData * integralOfChidtValue;

  chiValue = NULL;
  integralOfChidtValue = NULL;

  chi = 0.0;
  have_tau_thermostat = 0;
  have_target_temp = 0;
  have_chi_tolerance = 0;
  integralOfChidt = 0.0;


  if( theInfo->useInitXSstate ){

    // retrieve chi and integralOfChidt from simInfo
    data = info->getProperty(CHIVALUE_ID);
    if(data){
      chiValue = dynamic_cast<DoubleGenericData*>(data);
    }
    
    data = info->getProperty(INTEGRALOFCHIDT_ID);
    if(data){
      integralOfChidtValue = dynamic_cast<DoubleGenericData*>(data);
    }
    
    // chi and integralOfChidt should appear by pair
    if(chiValue && integralOfChidtValue){
      chi = chiValue->getData();
      integralOfChidt = integralOfChidtValue->getData();
    }
  }

  oldVel = new double[3*integrableObjects.size()];
  oldJi = new double[3*integrableObjects.size()];
}

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

template<typename T> void NVT<T>::moveA() {

  int i, j;
  DirectionalAtom* dAtom;
  Vector3d Tb;
  Vector3d ji;
  double mass;
  Vector3d vel;
  Vector3d pos;
  Vector3d frc;

  double instTemp;

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

  instTemp = tStats->getTemperature();

  for( i=0; i < integrableObjects.size(); i++ ){

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

    mass = integrableObjects[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];
    }

    integrableObjects[i]->setVel( vel );
    integrableObjects[i]->setPos( pos );

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

      // get and convert the torque to body frame

      Tb = integrableObjects[i]->getTrq();
      integrableObjects[i]->lab2Body( Tb );

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

      ji = integrableObjects[i]->getJ();

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

      this->rotationPropagation( integrableObjects[i], ji );

      integrableObjects[i]->setJ( ji );
    }
  }
  
  if(nConstrained)
    constrainA();

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

  //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
  
  chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
  integralOfChidt += chi*dt2;

}

template<typename T> void NVT<T>::moveB( void ){
  int i, j, k;
  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 < integrableObjects.size(); i++ ){

    vel = integrableObjects[i]->getVel();

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

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

      ji = integrableObjects[i]->getJ();

      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 < integrableObjects.size(); i++ ){

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

      mass = integrableObjects[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);

      integrableObjects[i]->setVel( vel );

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

        // get and convert the torque to body frame

        Tb = integrableObjects[i]->getTrq();
        integrableObjects[i]->lab2Body( Tb );

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

        integrableObjects[i]->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,
             "You can't use the NVT integrator without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    painCave.severity = OOPSE_ERROR;
    simError();
    return -1;
  }

  // We must set tauThermostat.

  if (!have_tau_thermostat) {
    sprintf( painCave.errMsg,
             "If you use the constant temperature\n"
             "\tintegrator, you must set tauThermostat.\n");
    painCave.severity = OOPSE_ERROR;
    painCave.isFatal = 1;
    simError();
    return -1;
  }

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

  return 1;

}

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

  double conservedQuantity;
  double fkBT;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;

  fkBT = (double)(info->ndf) * kB * targetTemp;

  Energy = tStats->getTotalE();

  thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
    (2.0 * eConvert);

  thermostat_potential = fkBT * integralOfChidt / eConvert;

  conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;

  return conservedQuantity;
}

template<typename T> string NVT<T>::getAdditionalParameters(void){
  string parameters;
  const int BUFFERSIZE = 2000; // size of the read buffer
  char buffer[BUFFERSIZE];

  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
  parameters += buffer;

  return parameters;
}
Revision:	1695
Committed:	Mon Nov 1 22:52:57 2004 UTC (19 years, 8 months ago) by tim
File size:	7153 byte(s)
Log Message:	Molecule, Atom, DirectionalAtom, RigidBody and StuntDouble classes get compiled
#	User	Rev	Content
1	tim	1695	#include <math.h>
2
3			#include "primitives/Atom.hpp"
4			#include "primitives/SRI.hpp"
5			#include "primitives/AbstractClasses.hpp"
6			#include "brains/SimInfo.hpp"
7			#include "UseTheForce/ForceFields.hpp"
8			#include "brains/Thermo.hpp"
9			#include "io/ReadWrite.hpp"
10			#include "integrators/Integrator.hpp"
11			#include "utils/simError.h"
12
13
14			// Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
15
16			template<typename T> NVT<T>::NVT ( SimInfo theInfo, ForceFields the_ff):
17			T( theInfo, the_ff )
18			{
19			GenericData* data;
20			DoubleGenericData * chiValue;
21			DoubleGenericData * integralOfChidtValue;
22
23			chiValue = NULL;
24			integralOfChidtValue = NULL;
25
26			chi = 0.0;
27			have_tau_thermostat = 0;
28			have_target_temp = 0;
29			have_chi_tolerance = 0;
30			integralOfChidt = 0.0;
31
32
33			if( theInfo->useInitXSstate ){
34
35			// retrieve chi and integralOfChidt from simInfo
36			data = info->getProperty(CHIVALUE_ID);
37			if(data){
38			chiValue = dynamic_cast<DoubleGenericData*>(data);
39			}
40
41			data = info->getProperty(INTEGRALOFCHIDT_ID);
42			if(data){
43			integralOfChidtValue = dynamic_cast<DoubleGenericData*>(data);
44			}
45
46			// chi and integralOfChidt should appear by pair
47			if(chiValue && integralOfChidtValue){
48			chi = chiValue->getData();
49			integralOfChidt = integralOfChidtValue->getData();
50			}
51			}
52
53			oldVel = new double[3*integrableObjects.size()];
54			oldJi = new double[3*integrableObjects.size()];
55			}
56
57			template<typename T> NVT<T>::~NVT() {
58			delete[] oldVel;
59			delete[] oldJi;
60			}
61
62			template<typename T> void NVT<T>::moveA() {
63
64			int i, j;
65			DirectionalAtom* dAtom;
66			Vector3d Tb;
67			Vector3d ji;
68			double mass;
69			Vector3d vel;
70			Vector3d pos;
71			Vector3d frc;
72
73			double instTemp;
74
75			// We need the temperature at time = t for the chi update below:
76
77			instTemp = tStats->getTemperature();
78
79			for( i=0; i < integrableObjects.size(); i++ ){
80
81			vel = integrableObjects[i]->getVel();
82			pos = integrableObjects[i]->getPos();
83			integrableObjects[i]->getFrc( frc );
84
85			mass = integrableObjects[i]->getMass();
86
87			for (j=0; j < 3; j++) {
88			// velocity half step (use chi from previous step here):
89			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
90			// position whole step
91			pos[j] += dt * vel[j];
92			}
93
94			integrableObjects[i]->setVel( vel );
95			integrableObjects[i]->setPos( pos );
96
97			if( integrableObjects[i]->isDirectional() ){
98
99			// get and convert the torque to body frame
100
101			Tb = integrableObjects[i]->getTrq();
102			integrableObjects[i]->lab2Body( Tb );
103
104			// get the angular momentum, and propagate a half step
105
106			ji = integrableObjects[i]->getJ();
107
108			for (j=0; j < 3; j++)
109			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
110
111			this->rotationPropagation( integrableObjects[i], ji );
112
113			integrableObjects[i]->setJ( ji );
114			}
115			}
116
117			if(nConstrained)
118			constrainA();
119
120			// Finally, evolve chi a half step (just like a velocity) using
121			// temperature at time t, not time t+dt/2
122
123			//std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
124
125			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
126			integralOfChidt += chi*dt2;
127
128			}
129
130			template<typename T> void NVT<T>::moveB( void ){
131			int i, j, k;
132			double Tb[3], ji[3];
133			double vel[3], frc[3];
134			double mass;
135			double instTemp;
136			double oldChi, prevChi;
137
138			// Set things up for the iteration:
139
140			oldChi = chi;
141
142			for( i=0; i < integrableObjects.size(); i++ ){
143
144			vel = integrableObjects[i]->getVel();
145
146			for (j=0; j < 3; j++)
147			oldVel[3*i + j] = vel[j];
148
149			if( integrableObjects[i]->isDirectional() ){
150
151			ji = integrableObjects[i]->getJ();
152
153			for (j=0; j < 3; j++)
154			oldJi[3*i + j] = ji[j];
155
156			}
157			}
158
159			// do the iteration:
160
161			for (k=0; k < 4; k++) {
162
163			instTemp = tStats->getTemperature();
164
165			// evolve chi another half step using the temperature at t + dt/2
166
167			prevChi = chi;
168			chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
169			(tauThermostat*tauThermostat);
170
171			for( i=0; i < integrableObjects.size(); i++ ){
172
173			integrableObjects[i]->getFrc( frc );
174			vel = integrableObjects[i]->getVel();
175
176			mass = integrableObjects[i]->getMass();
177
178			// velocity half step
179			for (j=0; j < 3; j++)
180			vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j]chi);
181
182			integrableObjects[i]->setVel( vel );
183
184			if( integrableObjects[i]->isDirectional() ){
185
186			// get and convert the torque to body frame
187
188			Tb = integrableObjects[i]->getTrq();
189			integrableObjects[i]->lab2Body( Tb );
190
191			for (j=0; j < 3; j++)
192			ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
193
194			integrableObjects[i]->setJ( ji );
195			}
196			}
197
198			if(nConstrained)
199			constrainB();
200
201			if (fabs(prevChi - chi) <= chiTolerance) break;
202			}
203
204			integralOfChidt += dt2*chi;
205			}
206
207			template<typename T> void NVT<T>::resetIntegrator( void ){
208
209			chi = 0.0;
210			integralOfChidt = 0.0;
211			}
212
213			template<typename T> int NVT<T>::readyCheck() {
214
215			//check parent's readyCheck() first
216			if (T::readyCheck() == -1)
217			return -1;
218
219			// 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			"You can't use the NVT integrator without a targetTemp!\n"
225			);
226			painCave.isFatal = 1;
227			painCave.severity = OOPSE_ERROR;
228			simError();
229			return -1;
230			}
231
232			// We must set tauThermostat.
233
234			if (!have_tau_thermostat) {
235			sprintf( painCave.errMsg,
236			"If you use the constant temperature\n"
237			"\tintegrator, you must set tauThermostat.\n");
238			painCave.severity = OOPSE_ERROR;
239			painCave.isFatal = 1;
240			simError();
241			return -1;
242			}
243
244			if (!have_chi_tolerance) {
245			sprintf( painCave.errMsg,
246			"In NVT integrator: setting chi tolerance to 1e-6\n");
247			chiTolerance = 1e-6;
248			have_chi_tolerance = 1;
249			painCave.severity = OOPSE_INFO;
250			painCave.isFatal = 0;
251			simError();
252			}
253
254			return 1;
255
256			}
257
258			template<typename T> double NVT<T>::getConservedQuantity(void){
259
260			double conservedQuantity;
261			double fkBT;
262			double Energy;
263			double thermostat_kinetic;
264			double thermostat_potential;
265
266			fkBT = (double)(info->ndf) * kB * targetTemp;
267
268			Energy = tStats->getTotalE();
269
270			thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
271			(2.0 * eConvert);
272
273			thermostat_potential = fkBT * integralOfChidt / eConvert;
274
275			conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
276
277			return conservedQuantity;
278			}
279
280			template<typename T> string NVT<T>::getAdditionalParameters(void){
281			string parameters;
282			const int BUFFERSIZE = 2000; // size of the read buffer
283			char buffer[BUFFERSIZE];
284
285			sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
286			parameters += buffer;
287
288			return parameters;
289			}