src/integrators/NPTf.cpp

#include <math.h>

#include "math/MatVec3.h"
#include "primitives/Atom.hpp"
#include "primitives/SRI.hpp"
#include "primitives/AbstractClasses.hpp"
#include "brains/SimInfo.hpp"
#include "UseTheForce/ForceFields.hpp"
#include "brains/Thermo.hpp"
#include "io/ReadWrite.hpp"
#include "integrators/Integrator.hpp"
#include "utils/simError.h"

#ifdef IS_MPI
#include "brains/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.

NPTf::NPTf (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);
          }
        }
      }
    }
  }

}

NPTf::~NPTf() {

  // empty for now
}

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

void NPTf::evolveEtaB() {

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

  
}

void NPTf::calcVelScale(){

  for (int i = 0; i < 3; i++ ) {
    for (int j = 0; j < 3; j++ ) {
      vScale(i, j) = eta(i, j);

      if (i == j) {
        vScale(i, j) += chi;
      }
    }
  }
}

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

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

void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM,
                           int index, Vector3d& sc) {

    /**@todo */
    Vector3d rj = (oldPos[index] + pos[j])/2.0 -COM;
    sc = eta * rj;
}

void NPTf::scaleSimBox(){

  int i;
  int j;
  int k;
  Mat3x3d scaleMat;
  double eta2ij;
  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++){

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

bool NPTf::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 NPTf::calcConservedQuantity(){

    double conservedQuantity;
    double totalEnergy;
    double thermostat_kinetic;
    double thermostat_potential;
    double barostat_kinetic;
    double barostat_potential;
    double trEta;

    totalEnergy = tStats->getTotalE();

    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * eConvert);

    thermostat_potential = fkBT* integralOfChidt / eConvert;

    trEta = (eta.transpose() * eta).trace();
    
    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;

}

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