src/integrators/NPTf.cpp

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

namespace oopse {

static IntegratorBuilder<NPTf>* NPTfCreator = new IntegratorBuilder<NPTf>("NPTf");

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

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

    chi= currentSnapshot_->getChi();
    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
    loadEta();
    
    // We need NkBT a lot, so just set it here: This is the RAW number
    // of integrableObjects, so no subtraction or addition of constraints or
    // orientational degrees of freedom:
    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;

    // fkBT is used because the thermostat operates on more degrees of freedom
    // than the barostat (when there are particles with orientational degrees
    // of freedom).  
    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
    
    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;

}

void NPTf::loadEta() {
    eta= currentSnapshot_->getEta();

    if (!eta.isDiagonal()) {
        sprintf( painCave.errMsg,
                 "NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
        painCave.isFatal = 1;
        simError();
    }
}

void NPTf::saveEta() {
    currentSnapshot_->setEta(eta);
}

}
Revision:	1867
Committed:	Tue Dec 7 23:08:14 2004 UTC (19 years, 7 months ago) by tim
File size:	7425 byte(s)
Log Message:	NPT in progress
#	User	Rev	Content
1	tim	1492	#include "brains/SimInfo.hpp"
2			#include "brains/Thermo.hpp"
3	tim	1837	#include "integrators/IntegratorCreator.hpp"
4	tim	1822	#include "integrators/NPTf.hpp"
5			#include "primitives/Molecule.hpp"
6			#include "utils/OOPSEConstant.hpp"
7	tim	1492	#include "utils/simError.h"
8	gezelter	1490
9	tim	1837	namespace oopse {
10
11			static IntegratorBuilder<NPTf>* NPTfCreator = new IntegratorBuilder<NPTf>("NPTf");
12
13	gezelter	1490	// Basic non-isotropic thermostating and barostating via the Melchionna
14			// modification of the Hoover algorithm:
15			//
16			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
17			// Molec. Phys., 78, 533.
18			//
19			// and
20			//
21			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
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	1867	chi= currentSnapshot_->getChi();
209			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
210			loadEta();
211
212			// We need NkBT a lot, so just set it here: This is the RAW number
213			// of integrableObjects, so no subtraction or addition of constraints or
214			// orientational degrees of freedom:
215			NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
216
217			// fkBT is used because the thermostat operates on more degrees of freedom
218			// than the barostat (when there are particles with orientational degrees
219			// of freedom).
220			fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
221
222	tim	1774	double conservedQuantity;
223			double totalEnergy;
224			double thermostat_kinetic;
225			double thermostat_potential;
226			double barostat_kinetic;
227			double barostat_potential;
228			double trEta;
229	gezelter	1490
230	tim	1822	totalEnergy = thermo.getTotalE();
231	gezelter	1490
232	tim	1822	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
233	gezelter	1490
234	tim	1822	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
235	gezelter	1490
236	tim	1822	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
237			trEta = tmp.trace();
238	tim	1774
239	tim	1822	barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
240	gezelter	1490
241	tim	1822	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
242	gezelter	1490
243	tim	1774	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
244			barostat_kinetic + barostat_potential;
245	gezelter	1490
246	tim	1774	return conservedQuantity;
247	gezelter	1490
248			}
249
250	tim	1837	void NPTf::loadEta() {
251			eta= currentSnapshot_->getEta();
252
253			if (!eta.isDiagonal()) {
254			sprintf( painCave.errMsg,
255			"NPTi error: the diagonal elements of are eta matrix is not same or etaMat is not a diagonal matrix");
256			painCave.isFatal = 1;
257			simError();
258			}
259			}
260
261			void NPTf::saveEta() {
262			currentSnapshot_->setEta(eta);
263			}
264
265			}