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

    // 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 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 NPTxyz::loadEta() {
    eta= currentSnapshot_->getEta();
}

}
Revision:	1867
Committed:	Tue Dec 7 23:08:14 2004 UTC (19 years, 7 months ago) by tim
File size:	3831 byte(s)
Log Message:	NPT in progress
#	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	tim	1867	// We need NkBT a lot, so just set it here: This is the RAW number
26			// of integrableObjects, so no subtraction or addition of constraints or
27			// orientational degrees of freedom:
28			NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
29	gezelter	1490
30	tim	1867	// fkBT is used because the thermostat operates on more degrees of freedom
31			// than the barostat (when there are particles with orientational degrees
32			// of freedom).
33			fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
34	gezelter	1490
35	tim	1867	double conservedQuantity;
36			double totalEnergy;
37			double thermostat_kinetic;
38			double thermostat_potential;
39			double barostat_kinetic;
40			double barostat_potential;
41			double trEta;
42	gezelter	1490
43	tim	1867	totalEnergy = thermo.getTotalE();
44	gezelter	1490
45	tim	1867	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
46
47			thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
48
49	tim	1822	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
50			trEta = tmp.trace();
51	gezelter	1490
52	tim	1867	barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
53	gezelter	1490
54	tim	1867	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
55	gezelter	1490
56	tim	1867	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
57			barostat_kinetic + barostat_potential;
58	gezelter	1490
59
60	tim	1867	return conservedQuantity;
61	gezelter	1490
62			}
63
64	tim	1774
65			void NPTxyz::scaleSimBox(){
66	gezelter	1490
67	tim	1867	int i,j,k;
68			Mat3x3d scaleMat;
69			double eta2ij, scaleFactor;
70			double bigScale, smallScale, offDiagMax;
71			Mat3x3d hm;
72			Mat3x3d hmnew;
73	gezelter	1490
74	tim	1774
75
76			// Scale the box after all the positions have been moved:
77
78			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
79			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
80
81	tim	1867	bigScale = 1.0;
82			smallScale = 1.0;
83			offDiagMax = 0.0;
84	tim	1774
85	tim	1867	for(i=0; i<3; i++){
86			for(j=0; j<3; j++){
87			scaleMat(i, j) = 0.0;
88			if(i==j) {
89			scaleMat(i, j) = 1.0;
90			}
91			}
92	tim	1774	}
93	gezelter	1490
94	tim	1867	for(i=0;i<3;i++){
95	gezelter	1490
96	tim	1774	// calculate the scaleFactors
97
98	tim	1867	scaleFactor = exp(dt*eta(i, i));
99	tim	1774
100	tim	1867	scaleMat(i, i) = scaleFactor;
101	tim	1774
102	tim	1867	if (scaleMat(i, i) > bigScale) {
103			bigScale = scaleMat(i, i);
104			}
105
106			if (scaleMat(i, i) < smallScale) {
107			smallScale = scaleMat(i, i);
108			}
109			}
110	tim	1774
111	tim	1867	if ((bigScale > 1.1) \|\| (smallScale < 0.9)) {
112			sprintf( painCave.errMsg,
113			"NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
114			" Check your tauBarostat, as it is probably too small!\n\n"
115			" scaleMat = [%lf\t%lf\t%lf]\n"
116			" [%lf\t%lf\t%lf]\n"
117			" [%lf\t%lf\t%lf]\n",
118			scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
119			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
120			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
121			painCave.isFatal = 1;
122			simError();
123			} else {
124	tim	1822
125			Mat3x3d hmat = currentSnapshot_->getHmat();
126			hmat = hmat *scaleMat;
127			currentSnapshot_->setHmat(hmat);
128	tim	1867	}
129	gezelter	1490	}
130	tim	1774
131	tim	1852	void NPTxyz::loadEta() {
132	tim	1837	eta= currentSnapshot_->getEta();
133	tim	1774	}
134	tim	1837
135			}