src/integrators/NPTf.cpp

#include "brains/SimInfo.hpp"
#include "brains/Thermo.hpp"
#include "integrators/NPTf.hpp"
#include "primitives/Molecule.hpp"
#include "utils/OOPSEConstant.hpp"
#include "utils/simError.h"

// 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){

}


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/OOPSEConstant::pressureConvert) / (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/OOPSEConstant::pressureConvert) / (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)/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 {

        Mat3x3d hmat = currentSnapshot_->getHmat();
        hmat = hmat *scaleMat;
        currentSnapshot_->setHmat(hmat);
        
  }
}

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 = thermo.getTotalE();

    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);

    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;

    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
    trEta = tmp.trace();
    
    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);

    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;

    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
        barostat_kinetic + barostat_potential;

    return conservedQuantity;

}

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