OOPSE/libmdtools/NVT.cpp

#include <math.h>

#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 )
{
  GenericData* data;
  DoubleData * chiValue;
  DoubleData * 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<DoubleData*>(data);
    }
    
    data = info->getProperty(INTEGRALOFCHIDT_ID);
    if(data){
      integralOfChidtValue = dynamic_cast<DoubleData*>(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;
  double Tb[3], ji[3];
  double 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 < integrableObjects.size(); i++ ){

    integrableObjects[i]->getVel( vel );
    integrableObjects[i]->getPos( pos );
    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

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

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

      integrableObjects[i]->getJ( ji );

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

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

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

  consFramework->doConstrainA();

  // 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++ ){

    integrableObjects[i]->getVel( vel );

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

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

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

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

      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

        integrableObjects[i]->getTrq( Tb );
        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 );
      }
    }

    consFramework->doConstrainB();

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