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 {

// 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,
    //             "NPTf error: the diagonal elements of eta matrix are not the same or etaMat is not a diagonal matrix");
    //    painCave.isFatal = 1;
    //    simError();
    //}
}

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

}
Revision:	1883
Committed:	Mon Dec 13 22:30:27 2004 UTC (19 years, 9 months ago) by tim
File size:	7354 byte(s)
Log Message:	MPI version is built
#	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	gezelter	1490	// Basic non-isotropic thermostating and barostating via the Melchionna
12			// modification of the Hoover algorithm:
13			//
14			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
15			// Molec. Phys., 78, 533.
16			//
17			// and
18			//
19			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
20
21	tim	1774	void NPTf::evolveEtaA() {
22	gezelter	1490
23			int i, j;
24
25	tim	1774	for(i = 0; i < 3; i ++){
26			for(j = 0; j < 3; j++){
27			if( i == j) {
28	tim	1822	eta(i, j) += dt2 * instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
29	tim	1774	} else {
30			eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2);
31			}
32			}
33	gezelter	1490	}
34
35	tim	1774	for(i = 0; i < 3; i++) {
36			for (j = 0; j < 3; j++) {
37			oldEta(i, j) = eta(i, j);
38			}
39			}
40
41	gezelter	1490	}
42
43	tim	1774	void NPTf::evolveEtaB() {
44	gezelter	1490
45	tim	1774	int i;
46			int j;
47	gezelter	1490
48	tim	1774	for(i = 0; i < 3; i++) {
49			for (j = 0; j < 3; j++) {
50			prevEta(i, j) = eta(i, j);
51			}
52			}
53	gezelter	1490
54	tim	1774	for(i = 0; i < 3; i ++){
55			for(j = 0; j < 3; j++){
56			if( i == j) {
57			eta(i, j) = oldEta(i, j) + dt2 * instaVol *
58	tim	1822	(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
59	tim	1774	} else {
60			eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2);
61			}
62			}
63	gezelter	1490	}
64	tim	1774
65
66	gezelter	1490	}
67
68	tim	1774	void NPTf::calcVelScale(){
69	gezelter	1490
70	tim	1774	for (int i = 0; i < 3; i++ ) {
71			for (int j = 0; j < 3; j++ ) {
72			vScale(i, j) = eta(i, j);
73	gezelter	1490
74			if (i == j) {
75	tim	1774	vScale(i, j) += chi;
76	gezelter	1490	}
77			}
78			}
79			}
80
81	tim	1822	void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
82	tim	1774	sc = vScale * vel;
83	gezelter	1490	}
84
85	tim	1774	void NPTf::getVelScaleB(Vector3d& sc, int index ) {
86			sc = vScale * oldVel[index];
87	gezelter	1490	}
88
89	tim	1822	void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
90	gezelter	1490
91	tim	1774	/*@todo /
92	tim	1822	Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
93	tim	1774	sc = eta * rj;
94	gezelter	1490	}
95
96	tim	1774	void NPTf::scaleSimBox(){
97	gezelter	1490
98	tim	1774	int i;
99			int j;
100			int k;
101			Mat3x3d scaleMat;
102	gezelter	1490	double eta2ij;
103			double bigScale, smallScale, offDiagMax;
104	tim	1774	Mat3x3d hm;
105			Mat3x3d hmnew;
106	gezelter	1490
107
108
109			// Scale the box after all the positions have been moved:
110
111			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
112			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
113
114			bigScale = 1.0;
115			smallScale = 1.0;
116			offDiagMax = 0.0;
117
118			for(i=0; i<3; i++){
119			for(j=0; j<3; j++){
120
121			// Calculate the matrix Product of the eta array (we only need
122			// the ij element right now):
123
124			eta2ij = 0.0;
125			for(k=0; k<3; k++){
126	tim	1774	eta2ij += eta(i, k) * eta(k, j);
127	gezelter	1490	}
128
129	tim	1774	scaleMat(i, j) = 0.0;
130	gezelter	1490	// identity matrix (see above):
131	tim	1774	if (i == j) scaleMat(i, j) = 1.0;
132	gezelter	1490	// Taylor expansion for the exponential truncated at second order:
133	tim	1774	scaleMat(i, j) += dteta(i, j) + 0.5dtdteta2ij;
134	gezelter	1490
135
136			if (i != j)
137	tim	1774	if (fabs(scaleMat(i, j)) > offDiagMax)
138			offDiagMax = fabs(scaleMat(i, j));
139	gezelter	1490	}
140
141	tim	1774	if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
142			if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
143	gezelter	1490	}
144
145			if ((bigScale > 1.01) \|\| (smallScale < 0.99)) {
146			sprintf( painCave.errMsg,
147			"NPTf error: Attempting a Box scaling of more than 1 percent.\n"
148			" Check your tauBarostat, as it is probably too small!\n\n"
149			" scaleMat = [%lf\t%lf\t%lf]\n"
150			" [%lf\t%lf\t%lf]\n"
151			" [%lf\t%lf\t%lf]\n"
152			" eta = [%lf\t%lf\t%lf]\n"
153			" [%lf\t%lf\t%lf]\n"
154			" [%lf\t%lf\t%lf]\n",
155	tim	1774	scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
156			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
157			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
158			eta(0, 0),eta(0, 1),eta(0, 2),
159			eta(1, 0),eta(1, 1),eta(1, 2),
160			eta(2, 0),eta(2, 1),eta(2, 2));
161	gezelter	1490	painCave.isFatal = 1;
162			simError();
163			} else if (offDiagMax > 0.01) {
164			sprintf( painCave.errMsg,
165			"NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
166			" Check your tauBarostat, as it is probably too small!\n\n"
167			" scaleMat = [%lf\t%lf\t%lf]\n"
168			" [%lf\t%lf\t%lf]\n"
169			" [%lf\t%lf\t%lf]\n"
170			" eta = [%lf\t%lf\t%lf]\n"
171			" [%lf\t%lf\t%lf]\n"
172			" [%lf\t%lf\t%lf]\n",
173	tim	1774	scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
174			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
175			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
176			eta(0, 0),eta(0, 1),eta(0, 2),
177			eta(1, 0),eta(1, 1),eta(1, 2),
178			eta(2, 0),eta(2, 1),eta(2, 2));
179	gezelter	1490	painCave.isFatal = 1;
180			simError();
181			} else {
182	tim	1822
183			Mat3x3d hmat = currentSnapshot_->getHmat();
184			hmat = hmat *scaleMat;
185			currentSnapshot_->setHmat(hmat);
186
187	gezelter	1490	}
188			}
189
190	tim	1774	bool NPTf::etaConverged() {
191			int i;
192			double diffEta, sumEta;
193	gezelter	1490
194	tim	1774	sumEta = 0;
195			for(i = 0; i < 3; i++) {
196			sumEta += pow(prevEta(i, i) - eta(i, i), 2);
197			}
198
199			diffEta = sqrt( sumEta / 3.0 );
200	gezelter	1490
201	tim	1774	return ( diffEta <= etaTolerance );
202	gezelter	1490	}
203
204	tim	1774	double NPTf::calcConservedQuantity(){
205	gezelter	1490
206	tim	1867	chi= currentSnapshot_->getChi();
207			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
208			loadEta();
209
210			// We need NkBT a lot, so just set it here: This is the RAW number
211			// of integrableObjects, so no subtraction or addition of constraints or
212			// orientational degrees of freedom:
213			NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
214
215			// fkBT is used because the thermostat operates on more degrees of freedom
216			// than the barostat (when there are particles with orientational degrees
217			// of freedom).
218			fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
219
220	tim	1774	double conservedQuantity;
221			double totalEnergy;
222			double thermostat_kinetic;
223			double thermostat_potential;
224			double barostat_kinetic;
225			double barostat_potential;
226			double trEta;
227	gezelter	1490
228	tim	1822	totalEnergy = thermo.getTotalE();
229	gezelter	1490
230	tim	1822	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
231	gezelter	1490
232	tim	1822	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
233	gezelter	1490
234	tim	1822	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
235			trEta = tmp.trace();
236	tim	1774
237	tim	1822	barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
238	gezelter	1490
239	tim	1822	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
240	gezelter	1490
241	tim	1774	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
242			barostat_kinetic + barostat_potential;
243	gezelter	1490
244	tim	1774	return conservedQuantity;
245	gezelter	1490
246			}
247
248	tim	1837	void NPTf::loadEta() {
249			eta= currentSnapshot_->getEta();
250
251	tim	1883	//if (!eta.isDiagonal()) {
252			// sprintf( painCave.errMsg,
253			// "NPTf error: the diagonal elements of eta matrix are not the same or etaMat is not a diagonal matrix");
254			// painCave.isFatal = 1;
255			// simError();
256			//}
257	tim	1837	}
258
259			void NPTf::saveEta() {
260			currentSnapshot_->setEta(eta);
261			}
262
263			}