OOPSE/libmdtools/NPTf.cpp

#include <math.h>

#include "MatVec3.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;
    }
  }


  if( theInfo->useInitXSstate ){
    // 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>::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 NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
 
  matVecMul3( vScale, vel, sc );
}

template<typename T> void NPTf<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 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];

  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.01) || (smallScale < 0.99)) {
    sprintf( painCave.errMsg,
             "NPTf error: Attempting a Box scaling of more than 1 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"
             "      eta = [%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],
             eta[0][0],eta[0][1],eta[0][2],
             eta[1][0],eta[1][1],eta[1][2],
             eta[2][0],eta[2][1],eta[2][2]);
    painCave.isFatal = 1;
    simError();
  } else if (offDiagMax > 0.01) {
    sprintf( painCave.errMsg,
             "NPTf error: Attempting an off-diagonal Box scaling of more than 1 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"
             "      eta = [%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],
             eta[0][0],eta[0][1],eta[0][2],
             eta[1][0],eta[1][1],eta[1][2],
             eta[2][0],eta[2][1],eta[2][2]);
    painCave.isFatal = 1;
    simError();
  } else {
    info->getBoxM(hm);
    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;

  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;

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