OOPSE/libmdtools/NPTxyz.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"

#ifdef IS_MPI
#include "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;
  DoubleArrayData * etaValue;
  vector<double> etaArray;
  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;
    }
  }

  // retrieve eta array from simInfo if it exists
  data = info->getProperty(ETAVALUE_ID);
  if(data){
    etaValue = dynamic_cast<DoubleArrayData*>(data);

    if(etaValue){
      etaArray = etaValue->getData();

      for(i = 0; i < 3; i++){
        for (j = 0; j < 3; j++){
          eta[i][j] = etaArray[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>::getVelScaleA(double sc[3], double vel[3]) {
  int i,j;
  double vScale[3][3];

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

  info->matVecMul3( vScale, vel, sc );
}

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

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

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

  info->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];

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

  info->transposeMat3(eta, a);
  info->matMul3(a, eta, b);
  trEta = info->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%f\t%f;", chi, integralOfChidt);
  parameters += buffer;

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

  return parameters;

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