src/integrators/NPTf.cpp

 /*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#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:	1927
Committed:	Wed Jan 12 17:14:35 2005 UTC (21 years, 4 months ago) by tim
File size:	9431 byte(s)
Log Message:	change license
#	User	Rev	Content
1	tim	1927	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42	tim	1492	#include "brains/SimInfo.hpp"
43			#include "brains/Thermo.hpp"
44	tim	1837	#include "integrators/IntegratorCreator.hpp"
45	tim	1822	#include "integrators/NPTf.hpp"
46			#include "primitives/Molecule.hpp"
47			#include "utils/OOPSEConstant.hpp"
48	tim	1492	#include "utils/simError.h"
49	gezelter	1490
50	tim	1837	namespace oopse {
51
52	gezelter	1490	// Basic non-isotropic thermostating and barostating via the Melchionna
53			// modification of the Hoover algorithm:
54			//
55			// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
56			// Molec. Phys., 78, 533.
57			//
58			// and
59			//
60			// Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
61
62	tim	1774	void NPTf::evolveEtaA() {
63	gezelter	1490
64			int i, j;
65
66	tim	1774	for(i = 0; i < 3; i ++){
67			for(j = 0; j < 3; j++){
68			if( i == j) {
69	tim	1822	eta(i, j) += dt2 * instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
70	tim	1774	} else {
71			eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2);
72			}
73			}
74	gezelter	1490	}
75
76	tim	1774	for(i = 0; i < 3; i++) {
77			for (j = 0; j < 3; j++) {
78			oldEta(i, j) = eta(i, j);
79			}
80			}
81
82	gezelter	1490	}
83
84	tim	1774	void NPTf::evolveEtaB() {
85	gezelter	1490
86	tim	1774	int i;
87			int j;
88	gezelter	1490
89	tim	1774	for(i = 0; i < 3; i++) {
90			for (j = 0; j < 3; j++) {
91			prevEta(i, j) = eta(i, j);
92			}
93			}
94	gezelter	1490
95	tim	1774	for(i = 0; i < 3; i ++){
96			for(j = 0; j < 3; j++){
97			if( i == j) {
98			eta(i, j) = oldEta(i, j) + dt2 * instaVol *
99	tim	1822	(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
100	tim	1774	} else {
101			eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2);
102			}
103			}
104	gezelter	1490	}
105	tim	1774
106
107	gezelter	1490	}
108
109	tim	1774	void NPTf::calcVelScale(){
110	gezelter	1490
111	tim	1774	for (int i = 0; i < 3; i++ ) {
112			for (int j = 0; j < 3; j++ ) {
113			vScale(i, j) = eta(i, j);
114	gezelter	1490
115			if (i == j) {
116	tim	1774	vScale(i, j) += chi;
117	gezelter	1490	}
118			}
119			}
120			}
121
122	tim	1822	void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
123	tim	1774	sc = vScale * vel;
124	gezelter	1490	}
125
126	tim	1774	void NPTf::getVelScaleB(Vector3d& sc, int index ) {
127			sc = vScale * oldVel[index];
128	gezelter	1490	}
129
130	tim	1822	void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
131	gezelter	1490
132	tim	1774	/*@todo /
133	tim	1822	Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
134	tim	1774	sc = eta * rj;
135	gezelter	1490	}
136
137	tim	1774	void NPTf::scaleSimBox(){
138	gezelter	1490
139	tim	1774	int i;
140			int j;
141			int k;
142			Mat3x3d scaleMat;
143	gezelter	1490	double eta2ij;
144			double bigScale, smallScale, offDiagMax;
145	tim	1774	Mat3x3d hm;
146			Mat3x3d hmnew;
147	gezelter	1490
148
149
150			// Scale the box after all the positions have been moved:
151
152			// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat)
153			// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2)
154
155			bigScale = 1.0;
156			smallScale = 1.0;
157			offDiagMax = 0.0;
158
159			for(i=0; i<3; i++){
160			for(j=0; j<3; j++){
161
162			// Calculate the matrix Product of the eta array (we only need
163			// the ij element right now):
164
165			eta2ij = 0.0;
166			for(k=0; k<3; k++){
167	tim	1774	eta2ij += eta(i, k) * eta(k, j);
168	gezelter	1490	}
169
170	tim	1774	scaleMat(i, j) = 0.0;
171	gezelter	1490	// identity matrix (see above):
172	tim	1774	if (i == j) scaleMat(i, j) = 1.0;
173	gezelter	1490	// Taylor expansion for the exponential truncated at second order:
174	tim	1774	scaleMat(i, j) += dteta(i, j) + 0.5dtdteta2ij;
175	gezelter	1490
176
177			if (i != j)
178	tim	1774	if (fabs(scaleMat(i, j)) > offDiagMax)
179			offDiagMax = fabs(scaleMat(i, j));
180	gezelter	1490	}
181
182	tim	1774	if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
183			if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
184	gezelter	1490	}
185
186			if ((bigScale > 1.01) \|\| (smallScale < 0.99)) {
187			sprintf( painCave.errMsg,
188			"NPTf error: Attempting a Box scaling of more than 1 percent.\n"
189			" Check your tauBarostat, as it is probably too small!\n\n"
190			" scaleMat = [%lf\t%lf\t%lf]\n"
191			" [%lf\t%lf\t%lf]\n"
192			" [%lf\t%lf\t%lf]\n"
193			" eta = [%lf\t%lf\t%lf]\n"
194			" [%lf\t%lf\t%lf]\n"
195			" [%lf\t%lf\t%lf]\n",
196	tim	1774	scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
197			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
198			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
199			eta(0, 0),eta(0, 1),eta(0, 2),
200			eta(1, 0),eta(1, 1),eta(1, 2),
201			eta(2, 0),eta(2, 1),eta(2, 2));
202	gezelter	1490	painCave.isFatal = 1;
203			simError();
204			} else if (offDiagMax > 0.01) {
205			sprintf( painCave.errMsg,
206			"NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
207			" Check your tauBarostat, as it is probably too small!\n\n"
208			" scaleMat = [%lf\t%lf\t%lf]\n"
209			" [%lf\t%lf\t%lf]\n"
210			" [%lf\t%lf\t%lf]\n"
211			" eta = [%lf\t%lf\t%lf]\n"
212			" [%lf\t%lf\t%lf]\n"
213			" [%lf\t%lf\t%lf]\n",
214	tim	1774	scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
215			scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
216			scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
217			eta(0, 0),eta(0, 1),eta(0, 2),
218			eta(1, 0),eta(1, 1),eta(1, 2),
219			eta(2, 0),eta(2, 1),eta(2, 2));
220	gezelter	1490	painCave.isFatal = 1;
221			simError();
222			} else {
223	tim	1822
224			Mat3x3d hmat = currentSnapshot_->getHmat();
225			hmat = hmat *scaleMat;
226			currentSnapshot_->setHmat(hmat);
227
228	gezelter	1490	}
229			}
230
231	tim	1774	bool NPTf::etaConverged() {
232			int i;
233			double diffEta, sumEta;
234	gezelter	1490
235	tim	1774	sumEta = 0;
236			for(i = 0; i < 3; i++) {
237			sumEta += pow(prevEta(i, i) - eta(i, i), 2);
238			}
239
240			diffEta = sqrt( sumEta / 3.0 );
241	gezelter	1490
242	tim	1774	return ( diffEta <= etaTolerance );
243	gezelter	1490	}
244
245	tim	1774	double NPTf::calcConservedQuantity(){
246	gezelter	1490
247	tim	1867	chi= currentSnapshot_->getChi();
248			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
249			loadEta();
250
251			// We need NkBT a lot, so just set it here: This is the RAW number
252			// of integrableObjects, so no subtraction or addition of constraints or
253			// orientational degrees of freedom:
254			NkBT = info_->getNGlobalIntegrableObjects()OOPSEConstant::kB targetTemp;
255
256			// fkBT is used because the thermostat operates on more degrees of freedom
257			// than the barostat (when there are particles with orientational degrees
258			// of freedom).
259			fkBT = info_->getNdf()OOPSEConstant::kB targetTemp;
260
261	tim	1774	double conservedQuantity;
262			double totalEnergy;
263			double thermostat_kinetic;
264			double thermostat_potential;
265			double barostat_kinetic;
266			double barostat_potential;
267			double trEta;
268	gezelter	1490
269	tim	1822	totalEnergy = thermo.getTotalE();
270	gezelter	1490
271	tim	1822	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
272	gezelter	1490
273	tim	1822	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
274	gezelter	1490
275	tim	1822	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
276			trEta = tmp.trace();
277	tim	1774
278	tim	1822	barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
279	gezelter	1490
280	tim	1822	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
281	gezelter	1490
282	tim	1774	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
283			barostat_kinetic + barostat_potential;
284	gezelter	1490
285	tim	1774	return conservedQuantity;
286	gezelter	1490
287			}
288
289	tim	1837	void NPTf::loadEta() {
290			eta= currentSnapshot_->getEta();
291
292	tim	1883	//if (!eta.isDiagonal()) {
293			// sprintf( painCave.errMsg,
294			// "NPTf error: the diagonal elements of eta matrix are not the same or etaMat is not a diagonal matrix");
295			// painCave.isFatal = 1;
296			// simError();
297			//}
298	tim	1837	}
299
300			void NPTf::saveEta() {
301			currentSnapshot_->setEta(eta);
302			}
303
304			}