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

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; 
  
}
Revision:	831
Committed:	Tue Oct 28 16:20:29 2003 UTC (22 years ago)
File size:	6845 byte(s)
Log Message:	This commit was manufactured by cvs2svn to create branch 'no-template- branch'.
#	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			#include "simError.h"
11
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			// Molec. Phys., 78, 533.
21			//
22			// and
23			//
24			// 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	mmeineke	780
30			int i,j;
31
32			for(i = 0; i < 3; i++){
33			for (j = 0; j < 3; j++){
34
35	gezelter	588	eta[i][j] = 0.0;
36	mmeineke	780	oldEta[i][j] = 0.0;
37			}
38			}
39			}
40	gezelter	588
41	mmeineke	780	template<typename T> NPTf<T>::~NPTf() {
42	tim	767
43	mmeineke	780	// empty for now
44			}
45	tim	767
46	mmeineke	780	template<typename T> void NPTf<T>::resetIntegrator() {
47
48			int i, j;
49
50			for(i = 0; i < 3; i++)
51			for (j = 0; j < 3; j++)
52			eta[i][j] = 0.0;
53
54			T::resetIntegrator();
55	gezelter	576	}
56
57	mmeineke	780	template<typename T> void NPTf<T>::evolveEtaA() {
58
59			int i, j;
60
61			for(i = 0; i < 3; i ++){
62			for(j = 0; j < 3; j++){
63			if( i == j)
64			eta[i][j] += dt2 * instaVol *
65			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
66			else
67			eta[i][j] += dt2 * instaVol * press[i][j] / (NkBT*tb2);
68			}
69			}
70
71			for(i = 0; i < 3; i++)
72			for (j = 0; j < 3; j++)
73			oldEta[i][j] = eta[i][j];
74	tim	767	}
75
76	mmeineke	780	template<typename T> void NPTf<T>::evolveEtaB() {
77
78			int i,j;
79	gezelter	772
80	mmeineke	780	for(i = 0; i < 3; i++)
81			for (j = 0; j < 3; j++)
82			prevEta[i][j] = eta[i][j];
83	mmeineke	778
84	mmeineke	780	for(i = 0; i < 3; i ++){
85			for(j = 0; j < 3; j++){
86			if( i == j) {
87			eta[i][j] = oldEta[i][j] + dt2 * instaVol *
88			(press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
89			} else {
90			eta[i][j] = oldEta[i][j] + dt2 * instaVol * press[i][j] / (NkBT*tb2);
91			}
92			}
93			}
94			}
95	gezelter	600
96	mmeineke	780	template<typename T> void NPTf<T>::getVelScaleA(double sc[3], double vel[3]) {
97			int i,j;
98			double vScale[3][3];
99	gezelter	576
100	gezelter	588	for (i = 0; i < 3; i++ ) {
101			for (j = 0; j < 3; j++ ) {
102	tim	767	vScale[i][j] = eta[i][j];
103
104	gezelter	588	if (i == j) {
105	tim	767	vScale[i][j] += chi;
106			}
107	gezelter	588	}
108			}
109	tim	767
110	mmeineke	780	info->matVecMul3( vScale, vel, sc );
111			}
112	gezelter	600
113	mmeineke	780	template<typename T> void NPTf<T>::getVelScaleB(double sc[3], int index ){
114			int i,j;
115			double myVel[3];
116			double vScale[3][3];
117	gezelter	600
118	mmeineke	780	for (i = 0; i < 3; i++ ) {
119			for (j = 0; j < 3; j++ ) {
120			vScale[i][j] = eta[i][j];
121	gezelter	576
122	mmeineke	780	if (i == j) {
123			vScale[i][j] += chi;
124			}
125	tim	767	}
126			}
127	gezelter	600
128	mmeineke	780	for (j = 0; j < 3; j++)
129			myVel[j] = oldVel[3*index + j];
130	tim	767
131	mmeineke	780	info->matVecMul3( vScale, myVel, sc );
132			}
133	tim	767
134	mmeineke	780	template<typename T> void NPTf<T>::getPosScale(double pos[3], double COM[3],
135			int index, double sc[3]){
136			int j;
137			double rj[3];
138	tim	767
139	mmeineke	780	for(j=0; j<3; j++)
140			rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
141	tim	767
142	mmeineke	780	info->matVecMul3( eta, rj, sc );
143			}
144	tim	767
145	mmeineke	780	template<typename T> void NPTf<T>::scaleSimBox( void ){
146	tim	767
147	mmeineke	780	int i,j,k;
148			double scaleMat[3][3];
149			double eta2ij;
150			double bigScale, smallScale, offDiagMax;
151			double hm[3][3], hmnew[3][3];
152
153	tim	767
154	mmeineke	780
155	gezelter	577	// Scale the box after all the positions have been moved:
156	gezelter	600
157	gezelter	578	// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
158			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
159	gezelter	600
160	gezelter	617	bigScale = 1.0;
161			smallScale = 1.0;
162			offDiagMax = 0.0;
163	gezelter	600
164	gezelter	578	for(i=0; i<3; i++){
165			for(j=0; j<3; j++){
166	gezelter	600
167	gezelter	588	// Calculate the matrix Product of the eta array (we only need
168			// the ij element right now):
169	gezelter	600
170	gezelter	588	eta2ij = 0.0;
171	gezelter	578	for(k=0; k<3; k++){
172	gezelter	588	eta2ij += eta[i][k] * eta[k][j];
173	gezelter	578	}
174	gezelter	588
175			scaleMat[i][j] = 0.0;
176			// identity matrix (see above):
177			if (i == j) scaleMat[i][j] = 1.0;
178			// Taylor expansion for the exponential truncated at second order:
179			scaleMat[i][j] += dteta[i][j] + 0.5dtdteta2ij;
180	gezelter	617
181			if (i != j)
182			if (fabs(scaleMat[i][j]) > offDiagMax)
183			offDiagMax = fabs(scaleMat[i][j]);
184	gezelter	578	}
185	gezelter	617
186			if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
187			if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
188	gezelter	578	}
189	gezelter	600
190	gezelter	617	if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
191			sprintf( painCave.errMsg,
192			"NPTf error: Attempting a Box scaling of more than 10 percent.\n"
193			" Check your tauBarostat, as it is probably too small!\n\n"
194			" scaleMat = [%lf\t%lf\t%lf]\n"
195			" [%lf\t%lf\t%lf]\n"
196			" [%lf\t%lf\t%lf]\n",
197			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
198			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
199			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
200			painCave.isFatal = 1;
201			simError();
202			} else if (offDiagMax > 0.1) {
203			sprintf( painCave.errMsg,
204			"NPTf error: Attempting an off-diagonal Box scaling of more than 10 percent.\n"
205			" Check your tauBarostat, as it is probably too small!\n\n"
206			" scaleMat = [%lf\t%lf\t%lf]\n"
207			" [%lf\t%lf\t%lf]\n"
208			" [%lf\t%lf\t%lf]\n",
209			scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
210			scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
211			scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
212			painCave.isFatal = 1;
213			simError();
214			} else {
215			info->getBoxM(hm);
216			info->matMul3(hm, scaleMat, hmnew);
217			info->setBoxM(hmnew);
218			}
219	gezelter	576	}
220
221	mmeineke	780	template<typename T> bool NPTf<T>::etaConverged() {
222			int i;
223			double diffEta, sumEta;
224	gezelter	600
225	mmeineke	780	sumEta = 0;
226	tim	767	for(i = 0; i < 3; i++)
227	mmeineke	780	sumEta += pow(prevEta[i][i] - eta[i][i], 2);
228	tim	767
229	mmeineke	780	diffEta = sqrt( sumEta / 3.0 );
230	tim	767
231	mmeineke	780	return ( diffEta <= etaTolerance );
232	gezelter	576	}
233
234	mmeineke	780	template<typename T> double NPTf<T>::getConservedQuantity(void){
235	mmeineke	746
236	tim	763	double conservedQuantity;
237	mmeineke	782	double totalEnergy;
238	gezelter	772	double thermostat_kinetic;
239			double thermostat_potential;
240			double barostat_kinetic;
241			double barostat_potential;
242			double trEta;
243			double a[3][3], b[3][3];
244	tim	763
245	mmeineke	782	totalEnergy = tStats->getTotalE();
246	tim	763
247	mmeineke	787	thermostat_kinetic = fkBT * tt2 * chi * chi /
248	gezelter	772	(2.0 * eConvert);
249	tim	763
250	gezelter	772	thermostat_potential = fkBT* integralOfChidt / eConvert;
251	tim	763
252	gezelter	772	info->transposeMat3(eta, a);
253			info->matMul3(a, eta, b);
254			trEta = info->matTrace3(b);
255	tim	767
256	mmeineke	780	barostat_kinetic = NkBT * tb2 * trEta /
257	gezelter	772	(2.0 * eConvert);
258
259			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
260			eConvert;
261	tim	767
262	mmeineke	782	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
263	gezelter	772	barostat_kinetic + barostat_potential;
264
265	mmeineke	780	// cout.width(8);
266			// cout.precision(8);
267	tim	767
268	mmeineke	780	// cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
269			// "\t" << thermostat_potential << "\t" << barostat_kinetic <<
270			// "\t" << barostat_potential << "\t" << conservedQuantity << endl;
271	gezelter	772
272			return conservedQuantity;
273	mmeineke	780
274	tim	763	}