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:	1930
Committed:	Wed Jan 12 22:41:40 2005 UTC (19 years, 6 months ago) by gezelter
File size:	9431 byte(s)
Log Message:	merging new_design branch into OOPSE-2.0
#	User	Rev	Content
1	gezelter	1930	/*
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	gezelter	1930	#include "integrators/IntegratorCreator.hpp"
45			#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	gezelter	1930	namespace oopse {
51	gezelter	1490
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	1930	void NPTf::evolveEtaA() {
63	gezelter	1490
64	gezelter	1930	int i, j;
65	gezelter	1490
66	gezelter	1930	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	1490	}
75	gezelter	1930
76			for(i = 0; i < 3; i++) {
77			for (j = 0; j < 3; j++) {
78			oldEta(i, j) = eta(i, j);
79			}
80	gezelter	1490	}
81	gezelter	1930
82	gezelter	1490	}
83
84	gezelter	1930	void NPTf::evolveEtaB() {
85	gezelter	1490
86	gezelter	1930	int i;
87			int j;
88	gezelter	1490
89	gezelter	1930	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	gezelter	1930	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	1490
106
107			}
108
109	gezelter	1930	void NPTf::calcVelScale(){
110	gezelter	1490
111	gezelter	1930	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	gezelter	1930	vScale(i, j) += chi;
117	gezelter	1490	}
118			}
119			}
120			}
121
122	gezelter	1930	void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
123			sc = vScale * vel;
124	gezelter	1490	}
125
126	gezelter	1930	void NPTf::getVelScaleB(Vector3d& sc, int index ) {
127			sc = vScale * oldVel[index];
128	gezelter	1490	}
129
130	gezelter	1930	void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
131	gezelter	1490
132	gezelter	1930	/*@todo /
133			Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
134			sc = eta * rj;
135	gezelter	1490	}
136
137	gezelter	1930	void NPTf::scaleSimBox(){
138	gezelter	1490
139	gezelter	1930	int i;
140			int j;
141			int k;
142			Mat3x3d scaleMat;
143	gezelter	1490	double eta2ij;
144			double bigScale, smallScale, offDiagMax;
145	gezelter	1930	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	gezelter	1930	eta2ij += eta(i, k) * eta(k, j);
168	gezelter	1490	}
169
170	gezelter	1930	scaleMat(i, j) = 0.0;
171	gezelter	1490	// identity matrix (see above):
172	gezelter	1930	if (i == j) scaleMat(i, j) = 1.0;
173	gezelter	1490	// Taylor expansion for the exponential truncated at second order:
174	gezelter	1930	scaleMat(i, j) += dteta(i, j) + 0.5dtdteta2ij;
175	gezelter	1490
176
177			if (i != j)
178	gezelter	1930	if (fabs(scaleMat(i, j)) > offDiagMax)
179			offDiagMax = fabs(scaleMat(i, j));
180	gezelter	1490	}
181
182	gezelter	1930	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	gezelter	1930	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	gezelter	1930	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	gezelter	1930
224			Mat3x3d hmat = currentSnapshot_->getHmat();
225			hmat = hmat *scaleMat;
226			currentSnapshot_->setHmat(hmat);
227
228	gezelter	1490	}
229			}
230
231	gezelter	1930	bool NPTf::etaConverged() {
232			int i;
233			double diffEta, sumEta;
234	gezelter	1490
235	gezelter	1930	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	gezelter	1930	return ( diffEta <= etaTolerance );
243	gezelter	1490	}
244
245	gezelter	1930	double NPTf::calcConservedQuantity(){
246	gezelter	1490
247	gezelter	1930	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	1490
256	gezelter	1930	// 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	1490
269	gezelter	1930	totalEnergy = thermo.getTotalE();
270	gezelter	1490
271	gezelter	1930	thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
272	gezelter	1490
273	gezelter	1930	thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
274	gezelter	1490
275	gezelter	1930	SquareMatrix<double, 3> tmp = eta.transpose() * eta;
276			trEta = tmp.trace();
277
278			barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
279	gezelter	1490
280	gezelter	1930	barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
281	gezelter	1490
282	gezelter	1930	conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
283			barostat_kinetic + barostat_potential;
284	gezelter	1490
285	gezelter	1930	return conservedQuantity;
286	gezelter	1490
287			}
288
289	gezelter	1930	void NPTf::loadEta() {
290			eta= currentSnapshot_->getEta();
291	gezelter	1490
292	gezelter	1930	//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	1490
300	gezelter	1930	void NPTf::saveEta() {
301			currentSnapshot_->setEta(eta);
302			}
303	gezelter	1490
304			}