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

NVT::NVT ( SimInfo *theInfo, ForceFields* the_ff):
  Integrator( theInfo, the_ff )
{
  GenericData* data;

  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 != NULL ){
    chi = data->getDval();
  }

  data = info->getProperty(INTEGRALOFCHIDT_ID);
  if(data != NULL ){
    integralOfChidt = data->getDval();
  }

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

NVT::~NVT() {
  delete[] oldVel;
  delete[] oldJi;
}

void NVT::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;

}

void NVT::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;
}

void NVT::resetIntegrator( void ){

  chi = 0.0;
  integralOfChidt = 0.0;
}

int NVT::readyCheck() {

  //check parent's readyCheck() first
  if (Integrator::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;

}

double NVT::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;
}

char* NVT::getAdditionalParameters(void){

  sprintf(addParamBuffer,
          "\t%G\t%G;",
          chi, integralOfChidt
          );

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