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, 7 months ago) by tim
File size:	4897 byte(s)
Log Message:	refactory NPT integrator
#	Content
1	#include "integrators/NPTxyz.hpp"
2
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	namespace oopse {
14	/*
15	NPTxyz::NPTxyz (SimInfo* info): NPT(info) {
16	GenericData* data;
17	DoubleVectorGenericData * etaValue;
18	int i,j;
19
20	for(i = 0; i < 3; i++){
21	for (j = 0; j < 3; j++){
22
23	eta(i, j) = 0.0;
24	oldEta(i, j) = 0.0;
25	}
26	}
27
28
29	if( theInfo->useInitXSstate ){
30
31	// retrieve eta array from simInfo if it exists
32	data = info->getPropertyByName(ETAVALUE_ID);
33	if(data){
34	etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
35
36	if(etaValue){
37
38	for(i = 0; i < 3; i++){
39	for (j = 0; j < 3; j++){
40	eta(i, j) = (etaValue)[3i+j];
41	oldEta(i, j) = eta(i, j);
42	}
43	}
44	}
45	}
46	}
47	}
48
49
50
51	void NPTxyz::evolveEtaA() {
52
53	int i, j;
54
55	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
65	for(i = 0; i < 3; i++) {
66	for (j = 0; j < 3; j++) {
67	oldEta(i, j) = eta(i, j);
68	}
69	}
70
71	}
72
73	void NPTxyz::evolveEtaB() {
74
75	int i,j;
76
77	for(i = 0; i < 3; i++)
78	for (j = 0; j < 3; j++)
79	prevEta(i, j) = eta(i, j);
80
81	for(i = 0; i < 3; i ++){
82	for(j = 0; j < 3; j++){
83	if( i == j) {
84	eta(i, j) = oldEta(i, j) + dt2 * instaVol *
85	(press(i, j) - targetPressure/p_convert) / (NkBT*tb2);
86	} else {
87	eta(i, j) = 0.0;
88	}
89	}
90	}
91	}
92
93	void NPTxyz::calcVelScale(void) {
94	int i,j;
95
96	for (i = 0; i < 3; i++ ) {
97	for (j = 0; j < 3; j++ ) {
98	vScale(i, j) = eta(i, j);
99
100	if (i == j) {
101	vScale(i, j) += chi;
102	}
103	}
104	}
105	}
106
107
108	void NPTxyz::getVelScaleA(Vector3d& sc, const Vector3d& vel) {
109	sc = vScale * vel;
110	}
111
112	void NPTxyz::getVelScaleB(Vector3d& sc, int index ) {
113	sc = vScale * oldVel[index];
114	}
115
116	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	}
122
123	bool NPTxyz::etaConverged() {
124	int i;
125	double diffEta, sumEta;
126
127	sumEta = 0;
128	for(i = 0; i < 3; i++)
129	sumEta += pow(prevEta(i, i) - eta(i, i), 2);
130
131	diffEta = sqrt( sumEta / 3.0 );
132
133	return ( diffEta <= etaTolerance );
134	}
135
136	*/
137
138	double NPTxyz::calcConservedQuantity(){
139
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	Mat3x3d a;
148	Mat3x3d b;
149
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
176	void NPTxyz::scaleSimBox(){
177
178	int i,j,k;
179	Mat3x3d scaleMat;
180	double eta2ij, scaleFactor;
181	double bigScale, smallScale, offDiagMax;
182	Mat3x3d hm;
183	Mat3x3d hmnew;
184
185
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	}
202
203	for(i=0;i<3;i++){
204
205	// 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	}
233
234	}