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 )
{
  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;

  // 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*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 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);

      this->rotationPropagation( dAtom, ji );

      dAtom->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

  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 fkBT;
  double Energy;
  double thermostat_kinetic;
  double thermostat_potential;

  fkBT = (double)(info->getNDF()    ) * 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;

  cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
      "\t" << thermostat_potential << "\t" << conservedQuantity << endl;

  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:	837
Committed:	Wed Oct 29 00:19:10 2003 UTC (22 years ago) by tim
File size:	6286 byte(s)
Log Message:	add chi and eta to the comment line of dump file.
#	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	tim	837	#include "simError.h"
10	mmeineke	561
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	tim	837	GenericData* data;
18			DoubleData * chiValue;
19			DoubleData * integralOfChidtValue;
20
21			chiValue = NULL;
22			integralOfChidtValue = NULL;
23
24	gezelter	565	chi = 0.0;
25	gezelter	560	have_tau_thermostat = 0;
26			have_target_temp = 0;
27	tim	763	have_chi_tolerance = 0;
28			integralOfChidt = 0.0;
29
30	tim	837	// retrieve chi and integralOfChidt from simInfo
31			data = info->getProperty(CHIVALUE_ID);
32			if(data){
33			chiValue = dynamic_cast<DoubleData*>(data);
34			}
35
36			data = info->getProperty(INTEGRALOFCHIDT_ID);
37			if(data){
38			integralOfChidtValue = dynamic_cast<DoubleData*>(data);
39			}
40
41			// chi and integralOfChidt should appear by pair
42			if(chiValue && integralOfChidtValue){
43			chi = chiValue->getData();
44			integralOfChidt = integralOfChidtValue->getData();
45			}
46
47	tim	763	oldVel = new double[3*nAtoms];
48			oldJi = new double[3*nAtoms];
49	gezelter	560	}
50
51	tim	763	template<typename T> NVT<T>::~NVT() {
52			delete[] oldVel;
53			delete[] oldJi;
54			}
55
56	tim	645	template<typename T> void NVT<T>::moveA() {
57	tim	837
58	gezelter	600	int i, j;
59	gezelter	560	DirectionalAtom* dAtom;
60	gezelter	600	double Tb[3], ji[3];
61	mmeineke	778	double mass;
62	gezelter	600	double vel[3], pos[3], frc[3];
63
64	gezelter	565	double instTemp;
65	gezelter	560
66	tim	763	// We need the temperature at time = t for the chi update below:
67
68	gezelter	565	instTemp = tStats->getTemperature();
69	tim	837
70	gezelter	560	for( i=0; i<nAtoms; i++ ){
71
72	gezelter	600	atoms[i]->getVel( vel );
73			atoms[i]->getPos( pos );
74			atoms[i]->getFrc( frc );
75
76			mass = atoms[i]->getMass();
77
78			for (j=0; j < 3; j++) {
79	tim	763	// velocity half step (use chi from previous step here):
80	gezelter	600	vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
81			// position whole step
82	gezelter	560	pos[j] += dt * vel[j];
83	gezelter	600	}
84	gezelter	560
85	gezelter	600	atoms[i]->setVel( vel );
86			atoms[i]->setPos( pos );
87	tim	837
88	gezelter	560	if( atoms[i]->isDirectional() ){
89
90			dAtom = (DirectionalAtom *)atoms[i];
91	tim	837
92	gezelter	560	// get and convert the torque to body frame
93	tim	837
94	gezelter	600	dAtom->getTrq( Tb );
95	gezelter	560	dAtom->lab2Body( Tb );
96	tim	837
97	gezelter	560	// get the angular momentum, and propagate a half step
98
99	gezelter	600	dAtom->getJ( ji );
100
101	tim	837	for (j=0; j < 3; j++)
102	gezelter	600	ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
103	tim	837
104	mmeineke	778	this->rotationPropagation( dAtom, ji );
105	tim	837
106	gezelter	600	dAtom->setJ( ji );
107	tim	837	}
108	gezelter	560	}
109	tim	837
110	mmeineke	768	if (nConstrained){
111			constrainA();
112			}
113	tim	763
114	tim	837	// Finally, evolve chi a half step (just like a velocity) using
115	tim	763	// temperature at time t, not time t+dt/2
116
117			chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
118			integralOfChidt += chi*dt2;
119
120	gezelter	560	}
121
122	tim	645	template<typename T> void NVT<T>::moveB( void ){
123	tim	763	int i, j, k;
124	gezelter	560	DirectionalAtom* dAtom;
125	gezelter	600	double Tb[3], ji[3];
126			double vel[3], frc[3];
127			double mass;
128	tim	763	double instTemp;
129			double oldChi, prevChi;
130	gezelter	600
131	tim	763	// Set things up for the iteration:
132
133			oldChi = chi;
134
135	gezelter	560	for( i=0; i<nAtoms; i++ ){
136	gezelter	600
137			atoms[i]->getVel( vel );
138
139	tim	763	for (j=0; j < 3; j++)
140			oldVel[3*i + j] = vel[j];
141	gezelter	600
142	gezelter	560	if( atoms[i]->isDirectional() ){
143	gezelter	600
144	gezelter	560	dAtom = (DirectionalAtom *)atoms[i];
145	gezelter	600
146	tim	763	dAtom->getJ( ji );
147	gezelter	600
148	tim	763	for (j=0; j < 3; j++)
149			oldJi[3*i + j] = ji[j];
150	gezelter	600
151	tim	763	}
152			}
153	gezelter	600
154	tim	763	// do the iteration:
155	gezelter	600
156	tim	763	for (k=0; k < 4; k++) {
157	tim	837
158	tim	763	instTemp = tStats->getTemperature();
159
160			// evolve chi another half step using the temperature at t + dt/2
161
162			prevChi = chi;
163	tim	837	chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
164	tim	763	(tauThermostat*tauThermostat);
165	tim	837
166	tim	763	for( i=0; i<nAtoms; i++ ){
167
168			atoms[i]->getFrc( frc );
169			atoms[i]->getVel(vel);
170	tim	837
171	tim	763	mass = atoms[i]->getMass();
172	tim	837
173	tim	763	// velocity half step
174	tim	837	for (j=0; j < 3; j++)
175	tim	763	vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j]chi);
176	tim	837
177	tim	763	atoms[i]->setVel( vel );
178	tim	837
179	tim	763	if( atoms[i]->isDirectional() ){
180	tim	837
181	tim	763	dAtom = (DirectionalAtom *)atoms[i];
182	tim	837
183			// get and convert the torque to body frame
184
185	tim	763	dAtom->getTrq( Tb );
186	tim	837	dAtom->lab2Body( Tb );
187
188			for (j=0; j < 3; j++)
189	tim	763	ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
190	tim	837
191	tim	763	dAtom->setJ( ji );
192			}
193			}
194	gezelter	600
195	mmeineke	768	if (nConstrained){
196			constrainB();
197			}
198
199	tim	763	if (fabs(prevChi - chi) <= chiTolerance) break;
200	gezelter	560	}
201	tim	837
202	tim	763	integralOfChidt += dt2*chi;
203	gezelter	560	}
204
205	mmeineke	746	template<typename T> void NVT<T>::resetIntegrator( void ){
206	tim	837
207	mmeineke	746	chi = 0.0;
208	tim	763	integralOfChidt = 0.0;
209	mmeineke	746	}
210
211	tim	645	template<typename T> int NVT<T>::readyCheck() {
212	tim	658
213			//check parent's readyCheck() first
214			if (T::readyCheck() == -1)
215			return -1;
216	tim	837
217			// First check to see if we have a target temperature.
218			// Not having one is fatal.
219
220	gezelter	560	if (!have_target_temp) {
221			sprintf( painCave.errMsg,
222			"NVT error: You can't use the NVT integrator without a targetTemp!\n"
223			);
224			painCave.isFatal = 1;
225			simError();
226			return -1;
227			}
228	tim	837
229	gezelter	565	// We must set tauThermostat.
230	tim	837
231	gezelter	565	if (!have_tau_thermostat) {
232	gezelter	560	sprintf( painCave.errMsg,
233	gezelter	565	"NVT error: If you use the constant temperature\n"
234			" integrator, you must set tauThermostat.\n");
235	gezelter	560	painCave.isFatal = 1;
236			simError();
237			return -1;
238	tim	837	}
239	tim	763
240			if (!have_chi_tolerance) {
241			sprintf( painCave.errMsg,
242			"NVT warning: setting chi tolerance to 1e-6\n");
243			chiTolerance = 1e-6;
244			have_chi_tolerance = 1;
245			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	837	fkBT = (double)(info->getNDF() ) * kB * targetTemp;
262
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			cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
273	tim	769	"\t" << thermostat_potential << "\t" << conservedQuantity << endl;
274
275	tim	837	return conservedQuantity;
276	tim	763	}
277	tim	837
278			template<typename T> string NVT<T>::getAdditionalParameters(void){
279			string parameters;
280			const int BUFFERSIZE = 2000; // size of the read buffer
281			char buffer[BUFFERSIZE];
282
283			sprintf(buffer,"\t%g\t%g;", chi, integralOfChidt);
284			parameters += buffer;
285
286			return parameters;
287			}