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. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. 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.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
#include "brains/SimInfo.hpp"
#include "brains/Thermo.hpp"
#include "integrators/IntegratorCreator.hpp"
#include "integrators/NPTf.hpp"
#include "primitives/Molecule.hpp"
#include "utils/PhysicalConstants.hpp"
#include "utils/simError.h"

namespace OpenMD {

  // 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/PhysicalConstants::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/PhysicalConstants::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) += thermostat.first;
        }
      }
    }
  }

  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)/(RealType)2.0 -COM;
    sc = eta * rj;
  }

  void NPTf::scaleSimBox(){

    int i;
    int j;
    int k;
    Mat3x3d scaleMat;
    RealType eta2ij;
    RealType 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 = snap->getHmat();
      hmat = hmat *scaleMat;
      snap->setHmat(hmat);
        
    }
  }

  bool NPTf::etaConverged() {
    int i;
    RealType 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 );
  }

  RealType NPTf::calcConservedQuantity(){
    
    thermostat = snap->getThermostat();
    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()*PhysicalConstants::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()*PhysicalConstants::kB *targetTemp;    
    
    RealType conservedQuantity;
    RealType totalEnergy;
    RealType thermostat_kinetic;
    RealType thermostat_potential;
    RealType barostat_kinetic;
    RealType barostat_potential;
    RealType trEta;

    totalEnergy = thermo.getTotalEnergy();
    
    thermostat_kinetic = fkBT * tt2 * thermostat.first * 
      thermostat.first /(2.0 * PhysicalConstants::energyConvert);

    thermostat_potential = fkBT* thermostat.second / PhysicalConstants::energyConvert;

    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
    trEta = tmp.trace();
    
    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * PhysicalConstants::energyConvert);

    barostat_potential = (targetPressure * thermo.getVolume() / PhysicalConstants::pressureConvert) /PhysicalConstants::energyConvert;

    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
      barostat_kinetic + barostat_potential;

    return conservedQuantity;

  }

  void NPTf::loadEta() {
    eta= snap->getBarostat();

    //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() {
    snap->setBarostat(eta);
  }

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