OOPSE/libmdtools/NPTi.cpp

#include <cmath>
#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 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.

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

template<typename T> void NPTi<T>::moveA() {
  
  int i, j;
  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;

  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<nAtoms; i++ ){
    atoms[i]->getVel( vel );
    atoms[i]->getPos( pos );
    atoms[i]->getFrc( frc );

    mass = atoms[i]->getMass();

    for (j=0; j < 3; j++) {
      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
      rj[j] = pos[j];
    }

    atoms[i]->setVel( vel );

    info->wrapVector(rj);

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

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

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

  // We must set tauBarostat.
   
  if (!have_tau_barostat) {
    sprintf( painCave.errMsg,
             "NPTi error: If you use the NPTi\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;
}
Revision:	658
Committed:	Thu Jul 31 15:35:07 2003 UTC (20 years, 11 months ago) by tim
File size:	6170 byte(s)
Log Message:	Added Z constraint.
#	User	Rev	Content
1	gezelter	578	#include <cmath>
2	gezelter	574	#include "Atom.hpp"
3			#include "SRI.hpp"
4			#include "AbstractClasses.hpp"
5			#include "SimInfo.hpp"
6			#include "ForceFields.hpp"
7			#include "Thermo.hpp"
8			#include "ReadWrite.hpp"
9			#include "Integrator.hpp"
10			#include "simError.h"
11
12
13			// Basic isotropic thermostating and barostating via the Melchionna
14			// modification of the Hoover algorithm:
15			//
16			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
17			// Molec. Phys., 78, 533.
18			//
19			// and
20			//
21			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
22
23	tim	645	template<typename T> NPTi<T>::NPTi ( SimInfo theInfo, ForceFields the_ff):
24			T( theInfo, the_ff )
25	gezelter	574	{
26			chi = 0.0;
27			eta = 0.0;
28			have_tau_thermostat = 0;
29			have_tau_barostat = 0;
30			have_target_temp = 0;
31			have_target_pressure = 0;
32			}
33
34	tim	645	template<typename T> void NPTi<T>::moveA() {
35	gezelter	574
36	gezelter	600	int i, j;
37	gezelter	574	DirectionalAtom* dAtom;
38	gezelter	600	double Tb[3], ji[3];
39			double A[3][3], I[3][3];
40			double angle, mass;
41			double vel[3], pos[3], frc[3];
42
43	gezelter	574	double rj[3];
44			double instaTemp, instaPress, instaVol;
45	gezelter	611	double tt2, tb2, scaleFactor;
46	gezelter	574
47			tt2 = tauThermostat * tauThermostat;
48			tb2 = tauBarostat * tauBarostat;
49
50			instaTemp = tStats->getTemperature();
51			instaPress = tStats->getPressure();
52			instaVol = tStats->getVolume();
53
54	mmeineke	586	// first evolve chi a half step
55	gezelter	574
56			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
57	mmeineke	586	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
58			(p_convertNkBTtb2));
59	gezelter	574
60			for( i=0; i<nAtoms; i++ ){
61	gezelter	600	atoms[i]->getVel( vel );
62			atoms[i]->getPos( pos );
63			atoms[i]->getFrc( frc );
64	gezelter	574
65	gezelter	600	mass = atoms[i]->getMass();
66	gezelter	574
67	gezelter	600	for (j=0; j < 3; j++) {
68			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
69			rj[j] = pos[j];
70			}
71
72			atoms[i]->setVel( vel );
73
74	gezelter	577	info->wrapVector(rj);
75	gezelter	574
76	gezelter	600	for (j = 0; j < 3; j++)
77			pos[j] += dt * (vel[j] + eta*rj[j]);
78
79			atoms[i]->setPos( pos );
80
81	gezelter	574	if( atoms[i]->isDirectional() ){
82
83			dAtom = (DirectionalAtom *)atoms[i];
84
85			// get and convert the torque to body frame
86
87	gezelter	600	dAtom->getTrq( Tb );
88	gezelter	574	dAtom->lab2Body( Tb );
89
90			// get the angular momentum, and propagate a half step
91
92	gezelter	600	dAtom->getJ( ji );
93
94			for (j=0; j < 3; j++)
95			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
96	gezelter	574
97			// use the angular velocities to propagate the rotation matrix a
98			// full time step
99	gezelter	600
100			dAtom->getA(A);
101			dAtom->getI(I);
102
103	gezelter	574	// rotate about the x-axis
104	gezelter	600	angle = dt2 * ji[0] / I[0][0];
105			this->rotate( 1, 2, angle, ji, A );
106
107	gezelter	574	// rotate about the y-axis
108	gezelter	600	angle = dt2 * ji[1] / I[1][1];
109			this->rotate( 2, 0, angle, ji, A );
110	gezelter	574
111			// rotate about the z-axis
112	gezelter	600	angle = dt * ji[2] / I[2][2];
113			this->rotate( 0, 1, angle, ji, A);
114	gezelter	574
115			// rotate about the y-axis
116	gezelter	600	angle = dt2 * ji[1] / I[1][1];
117			this->rotate( 2, 0, angle, ji, A );
118	gezelter	574
119			// rotate about the x-axis
120	gezelter	600	angle = dt2 * ji[0] / I[0][0];
121			this->rotate( 1, 2, angle, ji, A );
122	gezelter	574
123	gezelter	600	dAtom->setJ( ji );
124			dAtom->setA( A );
125			}
126
127	gezelter	574	}
128	gezelter	611
129	gezelter	577	// Scale the box after all the positions have been moved:
130	gezelter	600
131	gezelter	611	scaleFactor = exp(dt*eta);
132
133	mmeineke	614	if ((scaleFactor > 1.1) \|\| (scaleFactor < 0.9)) {
134	gezelter	611	sprintf( painCave.errMsg,
135			"NPTi error: Attempting a Box scaling of more than 10 percent"
136			" check your tauBarostat, as it is probably too small!\n"
137			" eta = %lf, scaleFactor = %lf\n", eta, scaleFactor
138			);
139			painCave.isFatal = 1;
140			simError();
141			} else {
142			info->scaleBox(exp(dt*eta));
143			}
144	mmeineke	614
145	gezelter	574	}
146
147	tim	645	template<typename T> void NPTi<T>::moveB( void ){
148	gezelter	600
149			int i, j;
150	gezelter	574	DirectionalAtom* dAtom;
151	gezelter	600	double Tb[3], ji[3];
152			double vel[3], frc[3];
153			double mass;
154
155	gezelter	574	double instaTemp, instaPress, instaVol;
156			double tt2, tb2;
157
158			tt2 = tauThermostat * tauThermostat;
159			tb2 = tauBarostat * tauBarostat;
160
161			instaTemp = tStats->getTemperature();
162			instaPress = tStats->getPressure();
163			instaVol = tStats->getVolume();
164
165			chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
166	mmeineke	586	eta += dt2 * ( instaVol * (instaPress - targetPressure) /
167			(p_convertNkBTtb2));
168	gezelter	574
169			for( i=0; i<nAtoms; i++ ){
170	gezelter	600
171			atoms[i]->getVel( vel );
172			atoms[i]->getFrc( frc );
173
174			mass = atoms[i]->getMass();
175
176	gezelter	574	// velocity half step
177	gezelter	600	for (j=0; j < 3; j++)
178			vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*(chi+eta));
179	gezelter	574
180	gezelter	600	atoms[i]->setVel( vel );
181
182	gezelter	574	if( atoms[i]->isDirectional() ){
183	gezelter	600
184	gezelter	574	dAtom = (DirectionalAtom *)atoms[i];
185	gezelter	600
186			// get and convert the torque to body frame
187
188			dAtom->getTrq( Tb );
189	gezelter	574	dAtom->lab2Body( Tb );
190	gezelter	600
191			// get the angular momentum, and propagate a half step
192
193			dAtom->getJ( ji );
194
195			for (j=0; j < 3; j++)
196			ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
197
198			dAtom->setJ( ji );
199	gezelter	574	}
200			}
201			}
202
203	tim	645	template<typename T> int NPTi<T>::readyCheck() {
204	tim	658
205			//check parent's readyCheck() first
206			if (T::readyCheck() == -1)
207			return -1;
208	gezelter	574
209			// First check to see if we have a target temperature.
210			// Not having one is fatal.
211
212			if (!have_target_temp) {
213			sprintf( painCave.errMsg,
214			"NPTi error: You can't use the NPTi integrator\n"
215			" without a targetTemp!\n"
216			);
217			painCave.isFatal = 1;
218			simError();
219			return -1;
220			}
221
222			if (!have_target_pressure) {
223			sprintf( painCave.errMsg,
224			"NPTi error: You can't use the NPTi integrator\n"
225			" without a targetPressure!\n"
226			);
227			painCave.isFatal = 1;
228			simError();
229			return -1;
230			}
231
232			// We must set tauThermostat.
233
234			if (!have_tau_thermostat) {
235			sprintf( painCave.errMsg,
236			"NPTi error: If you use the NPTi\n"
237			" integrator, you must set tauThermostat.\n");
238			painCave.isFatal = 1;
239			simError();
240			return -1;
241			}
242
243			// We must set tauBarostat.
244
245			if (!have_tau_barostat) {
246			sprintf( painCave.errMsg,
247			"NPTi error: If you use the NPTi\n"
248			" integrator, you must set tauBarostat.\n");
249			painCave.isFatal = 1;
250			simError();
251			return -1;
252			}
253
254			// We need NkBT a lot, so just set it here:
255
256			NkBT = (double)info->ndf * kB * targetTemp;
257
258			return 1;
259			}