src/integrators/NPTxyz.cpp

#include "brains/SimInfo.hpp"
#include "brains/Thermo.hpp"
#include "integrators/NPTxyz.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.

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

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

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

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