src/integrators/NPTxyz.cpp

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

static IntegratorBuilder<NPTxyz>* NPTxyzCreator = new IntegratorBuilder<NPTxyz>("NPTxyz");
    
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);
  }
}

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

}
Revision:	1837
Committed:	Thu Dec 2 22:15:31 2004 UTC (19 years, 7 months ago) by tim
File size:	3133 byte(s)
Log Message:	refine factory pattern to make it initialized correctly
#	User	Rev	Content
1	tim	1822	#include "brains/SimInfo.hpp"
2			#include "brains/Thermo.hpp"
3	tim	1837	#include "integrators/IntegratorCreator.hpp"
4	tim	1774	#include "integrators/NPTxyz.hpp"
5	tim	1822	#include "primitives/Molecule.hpp"
6			#include "utils/OOPSEConstant.hpp"
7			#include "utils/simError.h"
8	gezelter	1490
9			// Basic non-isotropic thermostating and barostating via the Melchionna
10			// modification of the Hoover algorithm:
11			//
12			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
13			// Molec. Phys., 78, 533.
14			//
15			// and
16			//
17			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
18
19	tim	1774	namespace oopse {
20	gezelter	1490
21	tim	1837	static IntegratorBuilder<NPTxyz>* NPTxyzCreator = new IntegratorBuilder<NPTxyz>("NPTxyz");
22	tim	1774
23			double NPTxyz::calcConservedQuantity(){
24	gezelter	1490
25			double conservedQuantity;
26			double totalEnergy;
27			double thermostat_kinetic;
28			double thermostat_potential;
29			double barostat_kinetic;
30			double barostat_potential;
31			double trEta;
32
33	tim	1822	totalEnergy = thermo.getTotalE();
34	gezelter	1490
35	tim	1822	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
36	gezelter	1490
37	tim	1822	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
38	gezelter	1490
39	tim	1822	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
40			trEta = tmp.trace();
41	gezelter	1490
42	tim	1822	barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
43	gezelter	1490
44	tim	1822	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
45	gezelter	1490
46			conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
47			barostat_kinetic + barostat_potential;
48
49
50			return conservedQuantity;
51
52			}
53
54	tim	1774
55			void NPTxyz::scaleSimBox(){
56	gezelter	1490
57	tim	1774	int i,j,k;
58			Mat3x3d scaleMat;
59			double eta2ij, scaleFactor;
60			double bigScale, smallScale, offDiagMax;
61			Mat3x3d hm;
62			Mat3x3d hmnew;
63	gezelter	1490
64	tim	1774
65
66			// Scale the box after all the positions have been moved:
67
68			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
69			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
70
71			bigScale = 1.0;
72			smallScale = 1.0;
73			offDiagMax = 0.0;
74
75			for(i=0; i<3; i++){
76			for(j=0; j<3; j++){
77			scaleMat(i, j) = 0.0;
78			if(i==j) scaleMat(i, j) = 1.0;
79			}
80	gezelter	1490	}
81
82	tim	1774	for(i=0;i<3;i++){
83	gezelter	1490
84	tim	1774	// calculate the scaleFactors
85
86			scaleFactor = exp(dt*eta(i, i));
87
88			scaleMat(i, i) = scaleFactor;
89
90			if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
91			if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
92			}
93
94			if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
95			sprintf( painCave.errMsg,
96			"NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
97			" Check your tauBarostat, as it is probably too small!\n\n"
98			" scaleMat = [%lf\t%lf\t%lf]\n"
99			" [%lf\t%lf\t%lf]\n"
100			" [%lf\t%lf\t%lf]\n",
101			scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
102			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
103			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
104			painCave.isFatal = 1;
105			simError();
106			} else {
107	tim	1822
108			Mat3x3d hmat = currentSnapshot_->getHmat();
109			hmat = hmat *scaleMat;
110			currentSnapshot_->setHmat(hmat);
111	tim	1774	}
112	gezelter	1490	}
113	tim	1774
114	tim	1837	void NPTf::loadEta() {
115			eta= currentSnapshot_->getEta();
116	tim	1774	}
117	tim	1837
118			}