src/integrators/NPTxyz.cpp

#include <math.h>
#include "math/MatVec3.h"
#include "primitives/Atom.hpp"
#include "primitives/SRI.hpp"
#include "primitives/AbstractClasses.hpp"
#include "brains/SimInfo.hpp"
#include "UseTheForce/ForceFields.hpp"
#include "brains/Thermo.hpp"
#include "io/ReadWrite.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"

#ifdef IS_MPI
#include "brains/mpiSimulation.hpp"
#endif

// Basic non-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> NPTxyz<T>::NPTxyz ( SimInfo *theInfo, ForceFields* the_ff):
  T( theInfo, the_ff )
{
  GenericData* data;
  DoubleVectorGenericData * etaValue;
  int i,j;

  for(i = 0; i < 3; i++){
    for (j = 0; j < 3; j++){

      eta[i][j] = 0.0;
      oldEta[i][j] = 0.0;
    }
  }


  if( theInfo->useInitXSstate ){

    // retrieve eta array from simInfo if it exists
    data = info->getProperty(ETAVALUE_ID);
    if(data){
      etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
      
      if(etaValue){
        
        for(i = 0; i < 3; i++){
          for (j = 0; j < 3; j++){
            eta[i][j] = (*etaValue)[3*i+j];
            oldEta[i][j] = eta[i][j];
          }
        }
      }
    }
  }
}

template<typename T> NPTxyz<T>::~NPTxyz() {

  // empty for now
}

template<typename T> void NPTxyz<T>::resetIntegrator() {

  int i, j;

  for(i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      eta[i][j] = 0.0;

  T::resetIntegrator();
}

template<typename T> void NPTxyz<T>::evolveEtaA() {

  int i, j;

  for(i = 0; i < 3; i ++){
    for(j = 0; j < 3; j++){
      if( i == j)
        eta[i][j] += dt2 *  instaVol *
          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
      else
        eta[i][j] = 0.0;
    }
  }

  for(i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      oldEta[i][j] = eta[i][j];
}

template<typename T> void NPTxyz<T>::evolveEtaB() {

  int i,j;

  for(i = 0; i < 3; i++)
    for (j = 0; j < 3; j++)
      prevEta[i][j] = eta[i][j];

  for(i = 0; i < 3; i ++){
    for(j = 0; j < 3; j++){
      if( i == j) {
        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
      } else {
        eta[i][j] = 0.0;
      }
    }
  }
}

template<typename T> void NPTxyz<T>::calcVelScale(void) {
  int i,j;

  for (i = 0; i < 3; i++ ) {
    for (j = 0; j < 3; j++ ) {
      vScale[i][j] = eta[i][j];

      if (i == j) {
        vScale[i][j] += chi;
      }
    }
  }
}

template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
  matVecMul3( vScale, vel, sc );
}

template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
  int j;
  double myVel[3];

  for (j = 0; j < 3; j++)
    myVel[j] = oldVel[3*index + j];

  matVecMul3( vScale, myVel, sc );
}

template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
                                               int index, double sc[3]){
  int j;
  double rj[3];

  for(j=0; j<3; j++)
    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];

  matVecMul3( eta, rj, sc );
}

template<typename T> void NPTxyz<T>::scaleSimBox( void ){

  int i,j,k;
  double scaleMat[3][3];
  double eta2ij, scaleFactor;
  double bigScale, smallScale, offDiagMax;
  double hm[3][3], hmnew[3][3];


  // Scale the box after all the positions have been moved:

  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)

  bigScale = 1.0;
  smallScale = 1.0;
  offDiagMax = 0.0;

  for(i=0; i<3; i++){
    for(j=0; j<3; j++){
      scaleMat[i][j] = 0.0;
      if(i==j) scaleMat[i][j] = 1.0;
    }
  }

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

    // calculate the scaleFactors

    scaleFactor = exp(dt*eta[i][i]);

    scaleMat[i][i] = scaleFactor;

    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
  }

//   for(i=0; i<3; i++){
//     for(j=0; j<3; j++){

//       // Calculate the matrix Product of the eta array (we only need
//       // the ij element right now):

//       eta2ij = 0.0;
//       for(k=0; k<3; k++){
//         eta2ij += eta[i][k] * eta[k][j];
//       }

//       scaleMat[i][j] = 0.0;
//       // identity matrix (see above):
//       if (i == j) scaleMat[i][j] = 1.0;
//       // Taylor expansion for the exponential truncated at second order:
//       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;

//       if (i != j)
//         if (fabs(scaleMat[i][j]) > offDiagMax)
//           offDiagMax = fabs(scaleMat[i][j]);
//     }

//     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
//     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
//   }

  if ((bigScale > 1.1) || (smallScale < 0.9)) {
    sprintf( painCave.errMsg,
             "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
             " Check your tauBarostat, as it is probably too small!\n\n"
             " scaleMat = [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n"
             "            [%lf\t%lf\t%lf]\n",
             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
    painCave.isFatal = 1;
    simError();
  } else {
    info->getBoxM(hm);
    matMul3(hm, scaleMat, hmnew);
    info->setBoxM(hmnew);
  }
}

template<typename T> bool NPTxyz<T>::etaConverged() {
  int i;
  double diffEta, sumEta;

  sumEta = 0;
  for(i = 0; i < 3; i++)
    sumEta += pow(prevEta[i][i] - eta[i][i], 2);

  diffEta = sqrt( sumEta / 3.0 );

  return ( diffEta <= etaTolerance );
}

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

  double conservedQuantity;
  double totalEnergy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;
  double trEta;
  double a[3][3], b[3][3];

  totalEnergy = tStats->getTotalE();

  thermostat_kinetic = fkBT * tt2 * chi * chi /
    (2.0 * eConvert);

  thermostat_potential = fkBT* integralOfChidt / eConvert;

  transposeMat3(eta, a);
  matMul3(a, eta, b);
  trEta = matTrace3(b);

  barostat_kinetic = NkBT * tb2 * trEta /
    (2.0 * eConvert);

  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
    eConvert;

  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
    barostat_kinetic + barostat_potential;

//   cout.width(8);
//   cout.precision(8);

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

  return conservedQuantity;

}

template<typename T> string NPTxyz<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;

  for(int i = 0; i < 3; i++){
    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
    parameters += buffer;
  }

  return parameters;

}
Revision:	132
Committed:	Thu Oct 21 16:22:01 2004 UTC (20 years, 6 months ago) by tim
Original Path:	*trunk/src/integrators/NPTxyz.cpp*
File size:	7281 byte(s)
Log Message:	replace old GebericData with new GenericData
#	User	Rev	Content
1	gezelter	2	#include <math.h>
2	tim	3	#include "math/MatVec3.h"
3			#include "primitives/Atom.hpp"
4			#include "primitives/SRI.hpp"
5			#include "primitives/AbstractClasses.hpp"
6			#include "brains/SimInfo.hpp"
7			#include "UseTheForce/ForceFields.hpp"
8			#include "brains/Thermo.hpp"
9			#include "io/ReadWrite.hpp"
10			#include "integrators/Integrator.hpp"
11			#include "utils/simError.h"
12	gezelter	2
13			#ifdef IS_MPI
14	tim	3	#include "brains/mpiSimulation.hpp"
15	gezelter	2	#endif
16
17			// Basic non-isotropic thermostating and barostating via the Melchionna
18			// modification of the Hoover algorithm:
19			//
20			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
21			// Molec. Phys., 78, 533.
22			//
23			// and
24			//
25			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
26
27			template<typename T> NPTxyz<T>::NPTxyz ( SimInfo theInfo, ForceFields the_ff):
28			T( theInfo, the_ff )
29			{
30			GenericData* data;
31	tim	132	DoubleVectorGenericData * etaValue;
32	gezelter	2	int i,j;
33
34			for(i = 0; i < 3; i++){
35			for (j = 0; j < 3; j++){
36
37			eta[i][j] = 0.0;
38			oldEta[i][j] = 0.0;
39			}
40			}
41
42
43			if( theInfo->useInitXSstate ){
44
45			// retrieve eta array from simInfo if it exists
46			data = info->getProperty(ETAVALUE_ID);
47			if(data){
48	tim	132	etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
49	gezelter	2
50			if(etaValue){
51
52			for(i = 0; i < 3; i++){
53			for (j = 0; j < 3; j++){
54	tim	132	eta[i][j] = (etaValue)[3i+j];
55	gezelter	2	oldEta[i][j] = eta[i][j];
56			}
57			}
58			}
59			}
60			}
61			}
62
63			template<typename T> NPTxyz<T>::~NPTxyz() {
64
65			// empty for now
66			}
67
68			template<typename T> void NPTxyz<T>::resetIntegrator() {
69
70			int i, j;
71
72			for(i = 0; i < 3; i++)
73			for (j = 0; j < 3; j++)
74			eta[i][j] = 0.0;
75
76			T::resetIntegrator();
77			}
78
79			template<typename T> void NPTxyz<T>::evolveEtaA() {
80
81			int i, j;
82
83			for(i = 0; i < 3; i ++){
84			for(j = 0; j < 3; j++){
85			if( i == j)
86			eta[i][j] += dt2 * instaVol *
87			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
88			else
89			eta[i][j] = 0.0;
90			}
91			}
92
93			for(i = 0; i < 3; i++)
94			for (j = 0; j < 3; j++)
95			oldEta[i][j] = eta[i][j];
96			}
97
98			template<typename T> void NPTxyz<T>::evolveEtaB() {
99
100			int i,j;
101
102			for(i = 0; i < 3; i++)
103			for (j = 0; j < 3; j++)
104			prevEta[i][j] = eta[i][j];
105
106			for(i = 0; i < 3; i ++){
107			for(j = 0; j < 3; j++){
108			if( i == j) {
109			eta[i][j] = oldEta[i][j] + dt2 * instaVol *
110			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
111			} else {
112			eta[i][j] = 0.0;
113			}
114			}
115			}
116			}
117
118			template<typename T> void NPTxyz<T>::calcVelScale(void) {
119			int i,j;
120
121			for (i = 0; i < 3; i++ ) {
122			for (j = 0; j < 3; j++ ) {
123			vScale[i][j] = eta[i][j];
124
125			if (i == j) {
126			vScale[i][j] += chi;
127			}
128			}
129			}
130			}
131
132			template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
133			matVecMul3( vScale, vel, sc );
134			}
135
136			template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
137			int j;
138			double myVel[3];
139
140			for (j = 0; j < 3; j++)
141			myVel[j] = oldVel[3*index + j];
142
143			matVecMul3( vScale, myVel, sc );
144			}
145
146			template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
147			int index, double sc[3]){
148			int j;
149			double rj[3];
150
151			for(j=0; j<3; j++)
152			rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
153
154			matVecMul3( eta, rj, sc );
155			}
156
157			template<typename T> void NPTxyz<T>::scaleSimBox( void ){
158
159			int i,j,k;
160			double scaleMat[3][3];
161			double eta2ij, scaleFactor;
162			double bigScale, smallScale, offDiagMax;
163			double hm[3][3], hmnew[3][3];
164
165
166
167			// Scale the box after all the positions have been moved:
168
169			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
170			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
171
172			bigScale = 1.0;
173			smallScale = 1.0;
174			offDiagMax = 0.0;
175
176			for(i=0; i<3; i++){
177			for(j=0; j<3; j++){
178			scaleMat[i][j] = 0.0;
179			if(i==j) scaleMat[i][j] = 1.0;
180			}
181			}
182
183			for(i=0;i<3;i++){
184
185			// calculate the scaleFactors
186
187			scaleFactor = exp(dt*eta[i][i]);
188
189			scaleMat[i][i] = scaleFactor;
190
191			if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
192			if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
193			}
194
195			// for(i=0; i<3; i++){
196			// for(j=0; j<3; j++){
197
198			// // Calculate the matrix Product of the eta array (we only need
199			// // the ij element right now):
200
201			// eta2ij = 0.0;
202			// for(k=0; k<3; k++){
203			// eta2ij += eta[i][k] * eta[k][j];
204			// }
205
206			// scaleMat[i][j] = 0.0;
207			// // identity matrix (see above):
208			// if (i == j) scaleMat[i][j] = 1.0;
209			// // Taylor expansion for the exponential truncated at second order:
210			// scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
211
212			// if (i != j)
213			// if (fabs(scaleMat[i][j]) > offDiagMax)
214			// offDiagMax = fabs(scaleMat[i][j]);
215			// }
216
217			// if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
218			// if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
219			// }
220
221			if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
222			sprintf( painCave.errMsg,
223			"NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
224			" Check your tauBarostat, as it is probably too small!\n\n"
225			" scaleMat = [%lf\t%lf\t%lf]\n"
226			" [%lf\t%lf\t%lf]\n"
227			" [%lf\t%lf\t%lf]\n",
228			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
229			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
230			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
231			painCave.isFatal = 1;
232			simError();
233			} else {
234			info->getBoxM(hm);
235			matMul3(hm, scaleMat, hmnew);
236			info->setBoxM(hmnew);
237			}
238			}
239
240			template<typename T> bool NPTxyz<T>::etaConverged() {
241			int i;
242			double diffEta, sumEta;
243
244			sumEta = 0;
245			for(i = 0; i < 3; i++)
246			sumEta += pow(prevEta[i][i] - eta[i][i], 2);
247
248			diffEta = sqrt( sumEta / 3.0 );
249
250			return ( diffEta <= etaTolerance );
251			}
252
253			template<typename T> double NPTxyz<T>::getConservedQuantity(void){
254
255			double conservedQuantity;
256			double totalEnergy;
257			double thermostat_kinetic;
258			double thermostat_potential;
259			double barostat_kinetic;
260			double barostat_potential;
261			double trEta;
262			double a[3][3], b[3][3];
263
264			totalEnergy = tStats->getTotalE();
265
266			thermostat_kinetic = fkBT * tt2 * chi * chi /
267			(2.0 * eConvert);
268
269			thermostat_potential = fkBT* integralOfChidt / eConvert;
270
271			transposeMat3(eta, a);
272			matMul3(a, eta, b);
273			trEta = matTrace3(b);
274
275			barostat_kinetic = NkBT * tb2 * trEta /
276			(2.0 * eConvert);
277
278			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
279			eConvert;
280
281			conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
282			barostat_kinetic + barostat_potential;
283
284			// cout.width(8);
285			// cout.precision(8);
286
287			// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
288			// "\t" << thermostat_potential << "\t" << barostat_kinetic <<
289			// "\t" << barostat_potential << "\t" << conservedQuantity << endl;
290
291			return conservedQuantity;
292
293			}
294
295			template<typename T> string NPTxyz<T>::getAdditionalParameters(void){
296			string parameters;
297			const int BUFFERSIZE = 2000; // size of the read buffer
298			char buffer[BUFFERSIZE];
299
300			sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
301			parameters += buffer;
302
303			for(int i = 0; i < 3; i++){
304			sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
305			parameters += buffer;
306			}
307
308			return parameters;
309
310			}