src/integrators/NPTxyz.cpp

#include "integrators/NPTxyz.hpp"

// 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.

namespace oopse {
/*
NPTxyz::NPTxyz (SimInfo* info): NPT(info) {
  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->getPropertyByName(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);
          }
        }
      }
    }
  }
}


void NPTxyz::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);
        }
    }
    
}

void NPTxyz::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;
      }
    }
  }
}

void NPTxyz::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;
      }
    }
  }
}


void NPTxyz::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
    sc = vScale * vel;
}

void NPTxyz::getVelScaleB(Vector3d& sc, int index ) {
  sc = vScale * oldVel[index];
}

void NPTxyz::getPosScale(const Vector3d& pos, const Vector3d& COM,
                           int index, Vector3d& sc) {
                           
    Vector3d rj = (oldPos[index] + pos[j])/2.0 -COM;
    sc = eta * rj;
}

bool NPTxyz::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 );
}

*/
    
double NPTxyz::calcConservedQuantity(){

  double conservedQuantity;
  double totalEnergy;
  double thermostat_kinetic;
  double thermostat_potential;
  double barostat_kinetic;
  double barostat_potential;
  double trEta;
  Mat3x3d a;
  Mat3x3d b;

  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;

}

    
void NPTxyz::scaleSimBox(){

  int i,j,k;
  Mat3x3d scaleMat;
  double eta2ij, scaleFactor;
  double bigScale, smallScale, offDiagMax;
  Mat3x3d hm;
  Mat3x3d hmnew;


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

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

}
Revision:	1774
Committed:	Tue Nov 23 23:12:23 2004 UTC (19 years, 9 months ago) by tim
File size:	4897 byte(s)
Log Message:	refactory NPT integrator
#	User	Rev	Content
1	tim	1774	#include "integrators/NPTxyz.hpp"
2	gezelter	1490
3			// Basic non-isotropic thermostating and barostating via the Melchionna
4			// modification of the Hoover algorithm:
5			//
6			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
7			// Molec. Phys., 78, 533.
8			//
9			// and
10			//
11			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
12
13	tim	1774	namespace oopse {
14			/*
15			NPTxyz::NPTxyz (SimInfo* info): NPT(info) {
16	gezelter	1490	GenericData* data;
17	tim	1625	DoubleVectorGenericData * etaValue;
18	gezelter	1490	int i,j;
19
20			for(i = 0; i < 3; i++){
21			for (j = 0; j < 3; j++){
22
23	tim	1774	eta(i, j) = 0.0;
24			oldEta(i, j) = 0.0;
25	gezelter	1490	}
26			}
27
28
29			if( theInfo->useInitXSstate ){
30
31			// retrieve eta array from simInfo if it exists
32	tim	1701	data = info->getPropertyByName(ETAVALUE_ID);
33	gezelter	1490	if(data){
34	tim	1625	etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
35	gezelter	1490
36			if(etaValue){
37
38			for(i = 0; i < 3; i++){
39			for (j = 0; j < 3; j++){
40	tim	1774	eta(i, j) = (etaValue)[3i+j];
41			oldEta(i, j) = eta(i, j);
42	gezelter	1490	}
43			}
44			}
45			}
46			}
47			}
48
49
50
51	tim	1774	void NPTxyz::evolveEtaA() {
52	gezelter	1490
53			int i, j;
54
55	tim	1774	for(i = 0; i < 3; i ++){
56			for(j = 0; j < 3; j++){
57			if( i == j) {
58			eta(i, j) += dt2 * instaVol (press(i, j) - targetPressure/p_convert) / (NkBTtb2);
59			} else {
60			eta(i, j) = 0.0;
61			}
62			}
63			}
64	gezelter	1490
65	tim	1774	for(i = 0; i < 3; i++) {
66			for (j = 0; j < 3; j++) {
67			oldEta(i, j) = eta(i, j);
68			}
69	gezelter	1490	}
70	tim	1774
71	gezelter	1490	}
72
73	tim	1774	void NPTxyz::evolveEtaB() {
74	gezelter	1490
75			int i,j;
76
77			for(i = 0; i < 3; i++)
78			for (j = 0; j < 3; j++)
79	tim	1774	prevEta(i, j) = eta(i, j);
80	gezelter	1490
81			for(i = 0; i < 3; i ++){
82			for(j = 0; j < 3; j++){
83			if( i == j) {
84	tim	1774	eta(i, j) = oldEta(i, j) + dt2 * instaVol *
85			(press(i, j) - targetPressure/p_convert) / (NkBT*tb2);
86	gezelter	1490	} else {
87	tim	1774	eta(i, j) = 0.0;
88	gezelter	1490	}
89			}
90			}
91			}
92
93	tim	1774	void NPTxyz::calcVelScale(void) {
94	gezelter	1490	int i,j;
95
96			for (i = 0; i < 3; i++ ) {
97			for (j = 0; j < 3; j++ ) {
98	tim	1774	vScale(i, j) = eta(i, j);
99	gezelter	1490
100			if (i == j) {
101	tim	1774	vScale(i, j) += chi;
102	gezelter	1490	}
103			}
104			}
105			}
106
107
108	tim	1774	void NPTxyz::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
109			sc = vScale * vel;
110	gezelter	1490	}
111
112	tim	1774	void NPTxyz::getVelScaleB(Vector3d& sc, int index ) {
113			sc = vScale * oldVel[index];
114	gezelter	1490	}
115
116	tim	1774	void NPTxyz::getPosScale(const Vector3d& pos, const Vector3d& COM,
117			int index, Vector3d& sc) {
118
119			Vector3d rj = (oldPos[index] + pos[j])/2.0 -COM;
120			sc = eta * rj;
121	gezelter	1490	}
122
123	tim	1774	bool NPTxyz::etaConverged() {
124	gezelter	1490	int i;
125			double diffEta, sumEta;
126
127			sumEta = 0;
128			for(i = 0; i < 3; i++)
129	tim	1774	sumEta += pow(prevEta(i, i) - eta(i, i), 2);
130	gezelter	1490
131			diffEta = sqrt( sumEta / 3.0 );
132
133			return ( diffEta <= etaTolerance );
134			}
135
136	tim	1774	*/
137
138			double NPTxyz::calcConservedQuantity(){
139	gezelter	1490
140			double conservedQuantity;
141			double totalEnergy;
142			double thermostat_kinetic;
143			double thermostat_potential;
144			double barostat_kinetic;
145			double barostat_potential;
146			double trEta;
147	tim	1774	Mat3x3d a;
148			Mat3x3d b;
149	gezelter	1490
150			totalEnergy = tStats->getTotalE();
151
152			thermostat_kinetic = fkBT * tt2 * chi * chi /
153			(2.0 * eConvert);
154
155			thermostat_potential = fkBT* integralOfChidt / eConvert;
156
157			transposeMat3(eta, a);
158			matMul3(a, eta, b);
159			trEta = matTrace3(b);
160
161			barostat_kinetic = NkBT * tb2 * trEta /
162			(2.0 * eConvert);
163
164			barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
165			eConvert;
166
167			conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
168			barostat_kinetic + barostat_potential;
169
170
171			return conservedQuantity;
172
173			}
174
175	tim	1774
176			void NPTxyz::scaleSimBox(){
177	gezelter	1490
178	tim	1774	int i,j,k;
179			Mat3x3d scaleMat;
180			double eta2ij, scaleFactor;
181			double bigScale, smallScale, offDiagMax;
182			Mat3x3d hm;
183			Mat3x3d hmnew;
184	gezelter	1490
185	tim	1774
186
187			// Scale the box after all the positions have been moved:
188
189			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
190			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
191
192			bigScale = 1.0;
193			smallScale = 1.0;
194			offDiagMax = 0.0;
195
196			for(i=0; i<3; i++){
197			for(j=0; j<3; j++){
198			scaleMat(i, j) = 0.0;
199			if(i==j) scaleMat(i, j) = 1.0;
200			}
201	gezelter	1490	}
202
203	tim	1774	for(i=0;i<3;i++){
204	gezelter	1490
205	tim	1774	// calculate the scaleFactors
206
207			scaleFactor = exp(dt*eta(i, i));
208
209			scaleMat(i, i) = scaleFactor;
210
211			if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
212			if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
213			}
214
215			if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
216			sprintf( painCave.errMsg,
217			"NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
218			" Check your tauBarostat, as it is probably too small!\n\n"
219			" scaleMat = [%lf\t%lf\t%lf]\n"
220			" [%lf\t%lf\t%lf]\n"
221			" [%lf\t%lf\t%lf]\n",
222			scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
223			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
224			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
225			painCave.isFatal = 1;
226			simError();
227			} else {
228			info->getBoxM(hm);
229			matMul3(hm, scaleMat, hmnew);
230			info->setBoxM(hmnew);
231			}
232	gezelter	1490	}
233	tim	1774
234			}