OOPSE/libmdtools/NPTf.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> NPTf<T>::NPTf ( 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> NPTf<T>::~NPTf() {

  // empty for now
}

template<typename T> void NPTf<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 NPTf<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] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
    }
  }

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

template<typename T> void NPTf<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] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
      }
    }
  }
}

template<typename T> void NPTf<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 NPTf<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 NPTf<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 NPTf<T>::scaleSimBox( void ){

  int i,j,k;
  double scaleMat[3][3];
  double eta2ij;
  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++){

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