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:	246
Committed:	Wed Jan 12 22:41:40 2005 UTC (20 years, 3 months ago) by gezelter
Original Path:	*trunk/src/integrators/NPTf.cpp*
File size:	9431 byte(s)
Log Message:	merging new_design branch into OOPSE-2.0
#	User	Rev	Content
1	gezelter	246	/*
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	3	#include "brains/SimInfo.hpp"
43			#include "brains/Thermo.hpp"
44	gezelter	246	#include "integrators/IntegratorCreator.hpp"
45			#include "integrators/NPTf.hpp"
46			#include "primitives/Molecule.hpp"
47			#include "utils/OOPSEConstant.hpp"
48	tim	3	#include "utils/simError.h"
49	gezelter	2
50	gezelter	246	namespace oopse {
51	gezelter	2
52			// 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	gezelter	246	void NPTf::evolveEtaA() {
63	gezelter	2
64	gezelter	246	int i, j;
65	gezelter	2
66	gezelter	246	for(i = 0; i < 3; i ++){
67			for(j = 0; j < 3; j++){
68			if( i == j) {
69			eta(i, j) += dt2 * instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
70			} else {
71			eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2);
72			}
73			}
74	gezelter	2	}
75	gezelter	246
76			for(i = 0; i < 3; i++) {
77			for (j = 0; j < 3; j++) {
78			oldEta(i, j) = eta(i, j);
79			}
80	gezelter	2	}
81	gezelter	246
82	gezelter	2	}
83
84	gezelter	246	void NPTf::evolveEtaB() {
85	gezelter	2
86	gezelter	246	int i;
87			int j;
88	gezelter	2
89	gezelter	246	for(i = 0; i < 3; i++) {
90			for (j = 0; j < 3; j++) {
91			prevEta(i, j) = eta(i, j);
92			}
93			}
94	gezelter	2
95	gezelter	246	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			(press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
100			} else {
101			eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2);
102			}
103			}
104			}
105	gezelter	2
106
107			}
108
109	gezelter	246	void NPTf::calcVelScale(){
110	gezelter	2
111	gezelter	246	for (int i = 0; i < 3; i++ ) {
112			for (int j = 0; j < 3; j++ ) {
113			vScale(i, j) = eta(i, j);
114	gezelter	2
115			if (i == j) {
116	gezelter	246	vScale(i, j) += chi;
117	gezelter	2	}
118			}
119			}
120			}
121
122	gezelter	246	void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
123			sc = vScale * vel;
124	gezelter	2	}
125
126	gezelter	246	void NPTf::getVelScaleB(Vector3d& sc, int index ) {
127			sc = vScale * oldVel[index];
128	gezelter	2	}
129
130	gezelter	246	void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
131	gezelter	2
132	gezelter	246	/*@todo /
133			Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
134			sc = eta * rj;
135	gezelter	2	}
136
137	gezelter	246	void NPTf::scaleSimBox(){
138	gezelter	2
139	gezelter	246	int i;
140			int j;
141			int k;
142			Mat3x3d scaleMat;
143	gezelter	2	double eta2ij;
144			double bigScale, smallScale, offDiagMax;
145	gezelter	246	Mat3x3d hm;
146			Mat3x3d hmnew;
147	gezelter	2
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	gezelter	246	eta2ij += eta(i, k) * eta(k, j);
168	gezelter	2	}
169
170	gezelter	246	scaleMat(i, j) = 0.0;
171	gezelter	2	// identity matrix (see above):
172	gezelter	246	if (i == j) scaleMat(i, j) = 1.0;
173	gezelter	2	// Taylor expansion for the exponential truncated at second order:
174	gezelter	246	scaleMat(i, j) += dteta(i, j) + 0.5dtdteta2ij;
175	gezelter	2
176
177			if (i != j)
178	gezelter	246	if (fabs(scaleMat(i, j)) > offDiagMax)
179			offDiagMax = fabs(scaleMat(i, j));
180	gezelter	2	}
181
182	gezelter	246	if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
183			if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
184	gezelter	2	}
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	gezelter	246	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	2	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	gezelter	246	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	2	painCave.isFatal = 1;
221			simError();
222			} else {
223	gezelter	246
224			Mat3x3d hmat = currentSnapshot_->getHmat();
225			hmat = hmat *scaleMat;
226			currentSnapshot_->setHmat(hmat);
227
228	gezelter	2	}
229			}
230
231	gezelter	246	bool NPTf::etaConverged() {
232			int i;
233			double diffEta, sumEta;
234	gezelter	2
235	gezelter	246	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	2
242	gezelter	246	return ( diffEta <= etaTolerance );
243	gezelter	2	}
244
245	gezelter	246	double NPTf::calcConservedQuantity(){
246	gezelter	2
247	gezelter	246	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	gezelter	2
256	gezelter	246	// 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			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	2
269	gezelter	246	totalEnergy = thermo.getTotalE();
270	gezelter	2
271	gezelter	246	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
272	gezelter	2
273	gezelter	246	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
274	gezelter	2
275	gezelter	246	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
276			trEta = tmp.trace();
277
278			barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
279	gezelter	2
280	gezelter	246	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
281	gezelter	2
282	gezelter	246	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
283			barostat_kinetic + barostat_potential;
284	gezelter	2
285	gezelter	246	return conservedQuantity;
286	gezelter	2
287			}
288
289	gezelter	246	void NPTf::loadEta() {
290			eta= currentSnapshot_->getEta();
291	gezelter	2
292	gezelter	246	//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			}
299	gezelter	2
300	gezelter	246	void NPTf::saveEta() {
301			currentSnapshot_->setEta(eta);
302			}
303	gezelter	2
304			}