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

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

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

    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;

 std::cerr << "ndf = " << info->getNDF() << " fkbt = " << fkBT << "\n";

  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:	1127
Committed:	Tue Apr 20 16:56:40 2004 UTC (21 years, 6 months ago) by tim
File size:	6719 byte(s)
Log Message:	fixed getCOMVel and velocitize at thermo
#	Content
1	#include <math.h>
2
3	#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	#include "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	DoubleData * chiValue;
21	DoubleData * 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<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	}
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	double Tb[3], ji[3];
67	double mass;
68	double vel[3], pos[3], frc[3];
69
70	double instTemp;
71
72	// We need the temperature at time = t for the chi update below:
73
74	instTemp = tStats->getTemperature();
75
76	for( i=0; i < integrableObjects.size(); i++ ){
77
78	integrableObjects[i]->getVel( vel );
79	integrableObjects[i]->getPos( pos );
80	integrableObjects[i]->getFrc( frc );
81
82	mass = integrableObjects[i]->getMass();
83
84	for (j=0; j < 3; j++) {
85	// velocity half step (use chi from previous step here):
86	vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
87	// position whole step
88	pos[j] += dt * vel[j];
89	}
90
91	integrableObjects[i]->setVel( vel );
92	integrableObjects[i]->setPos( pos );
93
94	if( integrableObjects[i]->isDirectional() ){
95
96	// get and convert the torque to body frame
97
98	integrableObjects[i]->getTrq( Tb );
99	integrableObjects[i]->lab2Body( Tb );
100
101	// get the angular momentum, and propagate a half step
102
103	integrableObjects[i]->getJ( ji );
104
105	for (j=0; j < 3; j++)
106	ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
107
108	this->rotationPropagation( integrableObjects[i], ji );
109
110	integrableObjects[i]->setJ( ji );
111	}
112	}
113
114	if (nConstrained){
115	constrainA();
116	}
117
118	// Finally, evolve chi a half step (just like a velocity) using
119	// temperature at time t, not time t+dt/2
120
121	std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
122
123	chi += dt2 * ( instTemp / targetTemp - 1.0) / (tauThermostat*tauThermostat);
124	integralOfChidt += chi*dt2;
125
126	}
127
128	template<typename T> void NVT<T>::moveB( void ){
129	int i, j, k;
130	double Tb[3], ji[3];
131	double vel[3], frc[3];
132	double mass;
133	double instTemp;
134	double oldChi, prevChi;
135
136	// Set things up for the iteration:
137
138	oldChi = chi;
139
140	for( i=0; i < integrableObjects.size(); i++ ){
141
142	integrableObjects[i]->getVel( vel );
143
144	for (j=0; j < 3; j++)
145	oldVel[3*i + j] = vel[j];
146
147	if( integrableObjects[i]->isDirectional() ){
148
149	integrableObjects[i]->getJ( ji );
150
151	for (j=0; j < 3; j++)
152	oldJi[3*i + j] = ji[j];
153
154	}
155	}
156
157	// do the iteration:
158
159	for (k=0; k < 4; k++) {
160
161	instTemp = tStats->getTemperature();
162
163	// evolve chi another half step using the temperature at t + dt/2
164
165	prevChi = chi;
166	chi = oldChi + dt2 * ( instTemp / targetTemp - 1.0) /
167	(tauThermostat*tauThermostat);
168
169	for( i=0; i < integrableObjects.size(); i++ ){
170
171	integrableObjects[i]->getFrc( frc );
172	integrableObjects[i]->getVel(vel);
173
174	mass = integrableObjects[i]->getMass();
175
176	// velocity half step
177	for (j=0; j < 3; j++)
178	vel[j] = oldVel[3i+j] + dt2 ((frc[j] / mass ) * eConvert - oldVel[3i + j]chi);
179
180	integrableObjects[i]->setVel( vel );
181
182	if( integrableObjects[i]->isDirectional() ){
183
184	// get and convert the torque to body frame
185
186	integrableObjects[i]->getTrq( Tb );
187	integrableObjects[i]->lab2Body( Tb );
188
189	for (j=0; j < 3; j++)
190	ji[j] = oldJi[3i + j] + dt2 (Tb[j] * eConvert - oldJi[3i+j]chi);
191
192	integrableObjects[i]->setJ( ji );
193	}
194	}
195
196	if (nConstrained){
197	constrainB();
198	}
199
200	if (fabs(prevChi - chi) <= chiTolerance) break;
201	}
202
203	integralOfChidt += dt2*chi;
204	}
205
206	template<typename T> void NVT<T>::resetIntegrator( void ){
207
208	chi = 0.0;
209	integralOfChidt = 0.0;
210	}
211
212	template<typename T> int NVT<T>::readyCheck() {
213
214	//check parent's readyCheck() first
215	if (T::readyCheck() == -1)
216	return -1;
217
218	// First check to see if we have a target temperature.
219	// Not having one is fatal.
220
221	if (!have_target_temp) {
222	sprintf( painCave.errMsg,
223	"NVT error: You can't use the NVT integrator without a targetTemp!\n"
224	);
225	painCave.isFatal = 1;
226	simError();
227	return -1;
228	}
229
230	// We must set tauThermostat.
231
232	if (!have_tau_thermostat) {
233	sprintf( painCave.errMsg,
234	"NVT error: If you use the constant temperature\n"
235	" integrator, you must set tauThermostat.\n");
236	painCave.isFatal = 1;
237	simError();
238	return -1;
239	}
240
241	if (!have_chi_tolerance) {
242	sprintf( painCave.errMsg,
243	"NVT warning: setting chi tolerance to 1e-6\n");
244	chiTolerance = 1e-6;
245	have_chi_tolerance = 1;
246	painCave.isFatal = 0;
247	simError();
248	}
249
250	return 1;
251
252	}
253
254	template<typename T> double NVT<T>::getConservedQuantity(void){
255
256	double conservedQuantity;
257	double fkBT;
258	double Energy;
259	double thermostat_kinetic;
260	double thermostat_potential;
261
262	fkBT = (double)(info->getNDF() ) * kB * targetTemp;
263
264	std::cerr << "ndf = " << info->getNDF() << " fkbt = " << fkBT << "\n";
265
266	Energy = tStats->getTotalE();
267
268	thermostat_kinetic = fkBT* tauThermostat * tauThermostat * chi * chi /
269	(2.0 * eConvert);
270
271	thermostat_potential = fkBT * integralOfChidt / eConvert;
272
273	conservedQuantity = Energy + thermostat_kinetic + thermostat_potential;
274
275	return conservedQuantity;
276	}
277
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	}