OOPSE/libmdtools/NPTim.cpp

#include <cmath>
#include "Atom.hpp"
#include "Molecule.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 isotropic thermostating and barostating via the Melchionna
// modification of the Hoover algorithm:
//
//    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
//       Molec. Phys., 78, 533. 
//
//           and
// 
//    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.

//  The NPTim variant scales the molecular center-of-mass coordinates
//  instead of the atomic coordinates

template<typename T> NPTim<T>::NPTim ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  chi = 0.0;
  eta = 0.0;
  integralOfChidt = 0.0;
  have_tau_thermostat = 0;
  have_tau_barostat = 0;
  have_target_temp = 0;
  have_target_pressure = 0;
}

template<typename T> void NPTim<T>::moveA() {
  
  int i, j, k;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double A[3][3], I[3][3];
  double angle, mass;
  double vel[3], pos[3], frc[3];

  double rj[3];
  double instaTemp, instaPress, instaVol;
  double tt2, tb2, scaleFactor;

  int nInMol;
  double rc[3];

  nMols = info->n_mol;
  myMolecules = info->molecules;

  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaVol = tStats->getVolume();
   
   // first evolve chi a half step
  
  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
                 (p_convert*NkBT*tb2));

  for( i = 0; i < nMols; i++) {

    myMolecules[i].getCOM(rc);

    nInMol = myMolecules[i].getNAtoms();
    myAtoms = myMolecules[i].getMyAtoms();

    // find the minimum image coordinates of the molecular centers of mass:

    info->wrapVector(rc);

    for (j = 0; j < nInMol; j++) {

      if(myAtoms[j] != NULL) {

        myAtoms[j]->getVel( vel );
        myAtoms[j]->getPos( pos );
        myAtoms[j]->getFrc( frc );

        mass = myAtoms[j]->getMass();

        for (k=0; k < 3; k++)
          vel[k] += dt2 * ((frc[k] / mass ) * eConvert - vel[k]*(chi+eta));

        myAtoms[j]->setVel( vel );

        for (k = 0; k < 3; k++) 
          pos[k] += dt * (vel[k] + eta*rc[k]);

        myAtoms[j]->setPos( pos );

        if( myAtoms[j]->isDirectional() ){

          dAtom = (DirectionalAtom *)myAtoms[j];
          
          // 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 (k=0; k < 3; k++) 
            ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
          
          // use the angular velocities to propagate the rotation matrix a
          // full time step
          
          dAtom->getA(A);
          dAtom->getI(I);
          
          // rotate about the x-axis      
          angle = dt2 * ji[0] / I[0][0];
          this->rotate( 1, 2, angle, ji, A ); 
          
          // rotate about the y-axis
          angle = dt2 * ji[1] / I[1][1];
          this->rotate( 2, 0, angle, ji, A );
          
          // rotate about the z-axis
          angle = dt * ji[2] / I[2][2];
          this->rotate( 0, 1, angle, ji, A);
          
          // rotate about the y-axis
          angle = dt2 * ji[1] / I[1][1];
          this->rotate( 2, 0, angle, ji, A );
          
          // rotate about the x-axis
          angle = dt2 * ji[0] / I[0][0];
          this->rotate( 1, 2, angle, ji, A );
          
          dAtom->setJ( ji );
          dAtom->setA( A  );    
        }                        
      }
    }
  }

  // Scale the box after all the positions have been moved:
  
  scaleFactor = exp(dt*eta);

  if (scaleFactor > 1.1 || scaleFactor < 0.9) {
    sprintf( painCave.errMsg,
             "NPTi error: Attempting a Box scaling of more than 10 percent"
             " check your tauBarostat, as it is probably too small!\n"
             " eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
             );
    painCave.isFatal = 1;
    simError();
  } else {         
    info->scaleBox(exp(dt*eta));      
  }
}

template<typename T> void NPTim<T>::moveB( void ){
  int i, j;
  DirectionalAtom* dAtom;
  double Tb[3], ji[3];
  double vel[3], frc[3];
  double mass;

  double instaTemp, instaPress, instaVol;
  double tt2, tb2;
 
  tt2 = tauThermostat * tauThermostat;
  tb2 = tauBarostat * tauBarostat;

  instaTemp = tStats->getTemperature();
  instaPress = tStats->getPressure();
  instaVol = tStats->getVolume();

  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
  eta += dt2 * ( instaVol * (instaPress - targetPressure) / 
                 (p_convert*NkBT*tb2));
  
  for( i=0; i<nAtoms; i++ ){
    
    atoms[i]->getVel( vel );
    atoms[i]->getFrc( frc );
    
    mass = atoms[i]->getMass();
    
    // velocity half step
    for (j=0; j < 3; j++) 
      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
    
    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 );
      
      // 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);     
      
      dAtom->setJ( ji );
    }
  }
}

template<typename T> void NPTim<T>::resetIntegrator() {
  chi = 0.0;
  eta = 0.0;
}

template<typename T> int NPTim<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,
             "NPTim error: You can't use the NPTim integrator\n"
             "   without a targetTemp!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }

  if (!have_target_pressure) {
    sprintf( painCave.errMsg,
             "NPTim error: You can't use the NPTim integrator\n"
             "   without a targetPressure!\n"
             );
    painCave.isFatal = 1;
    simError();
    return -1;
  }
  
  // We must set tauThermostat.
   
  if (!have_tau_thermostat) {
    sprintf( painCave.errMsg,
             "NPTim error: If you use the NPTim\n"
             "   integrator, you must set tauThermostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPTim error: If you use the NPTim\n"
             "   integrator, you must set tauBarostat.\n");
    painCave.isFatal = 1;
    simError();
    return -1;
  }    

  // We need NkBT a lot, so just set it here:

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

  return 1;
}

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

  double conservedQuantity;
  double tb2;
  double eta2;  


  //HNVE
  conservedQuantity = tStats->getTotalE();

  //HNVT
  conservedQuantity += (info->getNDF() * kB * targetTemp * 
                        (integralOfChidt + tauThermostat * tauThermostat * chi * chi /2))/ eConvert ;

  //HNPT
  tb2 = tauBarostat *tauBarostat;
  eta2 = eta * eta;
  
  conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
                        3*NkBT/2 * tb2 * eta2) / eConvert;
  
  return conservedQuantity; 
}
Revision:	763
Committed:	Mon Sep 15 16:52:02 2003 UTC (21 years ago) by tim
File size:	7606 byte(s)
Log Message:	add conserved quantity to statWriter fix bug of vector wrapping at NPTi
#	User	Rev	Content
1	gezelter	596	#include <cmath>
2			#include "Atom.hpp"
3	gezelter	604	#include "Molecule.hpp"
4	gezelter	596	#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 isotropic thermostating and barostating via the Melchionna
15			// modification of the Hoover algorithm:
16			//
17			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
18			// Molec. Phys., 78, 533.
19			//
20			// and
21			//
22			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
23
24			// The NPTim variant scales the molecular center-of-mass coordinates
25			// instead of the atomic coordinates
26
27	tim	645	template<typename T> NPTim<T>::NPTim ( SimInfo theInfo, ForceFields the_ff):
28			T( theInfo, the_ff )
29	gezelter	596	{
30			chi = 0.0;
31			eta = 0.0;
32	tim	763	integralOfChidt = 0.0;
33	gezelter	596	have_tau_thermostat = 0;
34			have_tau_barostat = 0;
35			have_target_temp = 0;
36			have_target_pressure = 0;
37			}
38
39	tim	645	template<typename T> void NPTim<T>::moveA() {
40	gezelter	596
41	gezelter	605	int i, j, k;
42	gezelter	596	DirectionalAtom* dAtom;
43	gezelter	604	double Tb[3], ji[3];
44			double A[3][3], I[3][3];
45			double angle, mass;
46			double vel[3], pos[3], frc[3];
47
48			double rj[3];
49	gezelter	596	double instaTemp, instaPress, instaVol;
50	gezelter	611	double tt2, tb2, scaleFactor;
51	gezelter	596
52	gezelter	604	int nInMol;
53			double rc[3];
54
55	gezelter	596	nMols = info->n_mol;
56	gezelter	604	myMolecules = info->molecules;
57	gezelter	596
58			tt2 = tauThermostat * tauThermostat;
59			tb2 = tauBarostat * tauBarostat;
60
61			instaTemp = tStats->getTemperature();
62			instaPress = tStats->getPressure();
63			instaVol = tStats->getVolume();
64
65			// first evolve chi a half step
66
67			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
68			eta += dt2 * ( instaVol * (instaPress - targetPressure) /
69			(p_convertNkBTtb2));
70
71			for( i = 0; i < nMols; i++) {
72
73	gezelter	604	myMolecules[i].getCOM(rc);
74	gezelter	596
75	gezelter	604	nInMol = myMolecules[i].getNAtoms();
76			myAtoms = myMolecules[i].getMyAtoms();
77	gezelter	596
78			// find the minimum image coordinates of the molecular centers of mass:
79
80			info->wrapVector(rc);
81
82			for (j = 0; j < nInMol; j++) {
83
84	gezelter	604	if(myAtoms[j] != NULL) {
85	gezelter	596
86	gezelter	605	myAtoms[j]->getVel( vel );
87			myAtoms[j]->getPos( pos );
88			myAtoms[j]->getFrc( frc );
89	gezelter	596
90	gezelter	605	mass = myAtoms[j]->getMass();
91	gezelter	596
92	gezelter	605	for (k=0; k < 3; k++)
93			vel[k] += dt2 * ((frc[k] / mass ) * eConvert - vel[k]*(chi+eta));
94	gezelter	596
95	gezelter	605	myAtoms[j]->setVel( vel );
96	gezelter	596
97	gezelter	605	for (k = 0; k < 3; k++)
98			pos[k] += dt * (vel[k] + eta*rc[k]);
99	gezelter	596
100	gezelter	605	myAtoms[j]->setPos( pos );
101	gezelter	596
102	gezelter	604	if( myAtoms[j]->isDirectional() ){
103	gezelter	596
104	gezelter	604	dAtom = (DirectionalAtom *)myAtoms[j];
105	gezelter	596
106			// get and convert the torque to body frame
107	gezelter	604
108			dAtom->getTrq( Tb );
109	gezelter	596	dAtom->lab2Body( Tb );
110	gezelter	604
111	gezelter	596	// get the angular momentum, and propagate a half step
112
113	gezelter	604	dAtom->getJ( ji );
114
115	gezelter	605	for (k=0; k < 3; k++)
116			ji[k] += dt2 * (Tb[k] * eConvert - ji[k]*chi);
117	gezelter	596
118			// use the angular velocities to propagate the rotation matrix a
119			// full time step
120
121	gezelter	604	dAtom->getA(A);
122			dAtom->getI(I);
123
124	gezelter	596	// rotate about the x-axis
125	gezelter	604	angle = dt2 * ji[0] / I[0][0];
126			this->rotate( 1, 2, angle, ji, A );
127	gezelter	596
128			// rotate about the y-axis
129	gezelter	604	angle = dt2 * ji[1] / I[1][1];
130			this->rotate( 2, 0, angle, ji, A );
131	gezelter	596
132			// rotate about the z-axis
133	gezelter	604	angle = dt * ji[2] / I[2][2];
134			this->rotate( 0, 1, angle, ji, A);
135	gezelter	596
136			// rotate about the y-axis
137	gezelter	604	angle = dt2 * ji[1] / I[1][1];
138			this->rotate( 2, 0, angle, ji, A );
139	gezelter	596
140			// rotate about the x-axis
141	gezelter	604	angle = dt2 * ji[0] / I[0][0];
142			this->rotate( 1, 2, angle, ji, A );
143	gezelter	596
144	gezelter	604	dAtom->setJ( ji );
145			dAtom->setA( A );
146			}
147	gezelter	596	}
148			}
149			}
150	gezelter	611
151	gezelter	596	// Scale the box after all the positions have been moved:
152
153	gezelter	611	scaleFactor = exp(dt*eta);
154
155			if (scaleFactor > 1.1 \|\| scaleFactor < 0.9) {
156			sprintf( painCave.errMsg,
157			"NPTi error: Attempting a Box scaling of more than 10 percent"
158			" check your tauBarostat, as it is probably too small!\n"
159			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
160			);
161			painCave.isFatal = 1;
162			simError();
163			} else {
164			info->scaleBox(exp(dt*eta));
165			}
166	gezelter	596	}
167
168	tim	645	template<typename T> void NPTim<T>::moveB( void ){
169	gezelter	604	int i, j;
170	gezelter	596	DirectionalAtom* dAtom;
171	gezelter	604	double Tb[3], ji[3];
172			double vel[3], frc[3];
173			double mass;
174
175	gezelter	596	double instaTemp, instaPress, instaVol;
176			double tt2, tb2;
177	gezelter	604
178	gezelter	596	tt2 = tauThermostat * tauThermostat;
179			tb2 = tauBarostat * tauBarostat;
180
181			instaTemp = tStats->getTemperature();
182			instaPress = tStats->getPressure();
183			instaVol = tStats->getVolume();
184
185			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
186			eta += dt2 * ( instaVol * (instaPress - targetPressure) /
187			(p_convertNkBTtb2));
188
189			for( i=0; i<nAtoms; i++ ){
190
191	gezelter	604	atoms[i]->getVel( vel );
192			atoms[i]->getFrc( frc );
193
194			mass = atoms[i]->getMass();
195
196	gezelter	596	// velocity half step
197	gezelter	604	for (j=0; j < 3; j++)
198			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
199	gezelter	596
200	gezelter	604	atoms[i]->setVel( vel );
201
202	gezelter	596	if( atoms[i]->isDirectional() ){
203
204			dAtom = (DirectionalAtom *)atoms[i];
205
206	gezelter	604	// get and convert the torque to body frame
207	gezelter	596
208	gezelter	604	dAtom->getTrq( Tb );
209	gezelter	596	dAtom->lab2Body( Tb );
210
211	gezelter	604	// get the angular momentum, and propagate a half step
212	gezelter	596
213	gezelter	604	dAtom->getJ( ji );
214	gezelter	596
215	gezelter	604	for (j=0; j < 3; j++)
216			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
217	gezelter	596
218	gezelter	604	dAtom->setJ( ji );
219	gezelter	596	}
220			}
221			}
222
223	mmeineke	746	template<typename T> void NPTim<T>::resetIntegrator() {
224			chi = 0.0;
225			eta = 0.0;
226			}
227
228	tim	645	template<typename T> int NPTim<T>::readyCheck() {
229	tim	658
230			//check parent's readyCheck() first
231			if (T::readyCheck() == -1)
232			return -1;
233	gezelter	596
234			// First check to see if we have a target temperature.
235			// Not having one is fatal.
236
237			if (!have_target_temp) {
238			sprintf( painCave.errMsg,
239	gezelter	604	"NPTim error: You can't use the NPTim integrator\n"
240	gezelter	596	" without a targetTemp!\n"
241			);
242			painCave.isFatal = 1;
243			simError();
244			return -1;
245			}
246
247			if (!have_target_pressure) {
248			sprintf( painCave.errMsg,
249	gezelter	604	"NPTim error: You can't use the NPTim integrator\n"
250	gezelter	596	" without a targetPressure!\n"
251			);
252			painCave.isFatal = 1;
253			simError();
254			return -1;
255			}
256
257			// We must set tauThermostat.
258
259			if (!have_tau_thermostat) {
260			sprintf( painCave.errMsg,
261	gezelter	604	"NPTim error: If you use the NPTim\n"
262	gezelter	596	" integrator, you must set tauThermostat.\n");
263			painCave.isFatal = 1;
264			simError();
265			return -1;
266			}
267
268			// We must set tauBarostat.
269
270			if (!have_tau_barostat) {
271			sprintf( painCave.errMsg,
272	gezelter	604	"NPTim error: If you use the NPTim\n"
273	gezelter	596	" integrator, you must set tauBarostat.\n");
274			painCave.isFatal = 1;
275			simError();
276			return -1;
277			}
278
279			// We need NkBT a lot, so just set it here:
280
281			NkBT = (double)info->ndf * kB * targetTemp;
282
283			return 1;
284			}
285	tim	763
286			template<typename T> double NPTim<T>::getConservedQuantity(void){
287
288			double conservedQuantity;
289			double tb2;
290			double eta2;
291
292
293			//HNVE
294			conservedQuantity = tStats->getTotalE();
295
296			//HNVT
297			conservedQuantity += (info->getNDF() * kB * targetTemp *
298			(integralOfChidt + tauThermostat * tauThermostat * chi * chi /2))/ eConvert ;
299
300			//HNPT
301			tb2 = tauBarostat *tauBarostat;
302			eta2 = eta * eta;
303
304			conservedQuantity += (targetPressure * tStats->getVolume() / p_convert +
305			3NkBT/2 tb2 * eta2) / eConvert;
306
307			return conservedQuantity;
308			}