ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/integrators/NPT.cpp
(Generate patch)

Comparing trunk/src/integrators/NPT.cpp (file contents):
Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 665 by tim, Thu Oct 13 22:26:47 2005 UTC

# Line 1 | Line 1
1 + /*
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   #include <math.h>
43  
3 – #include "primitives/Atom.hpp"
4 – #include "primitives/SRI.hpp"
5 – #include "primitives/AbstractClasses.hpp"
44   #include "brains/SimInfo.hpp"
7 – #include "UseTheForce/ForceFields.hpp"
45   #include "brains/Thermo.hpp"
46 < #include "io/ReadWrite.hpp"
47 < #include "integrators/Integrator.hpp"
46 > #include "integrators/NPT.hpp"
47 > #include "math/SquareMatrix3.hpp"
48 > #include "primitives/Molecule.hpp"
49 > #include "utils/OOPSEConstant.hpp"
50   #include "utils/simError.h"
51  
13 – #ifdef IS_MPI
14 – #include "brains/mpiSimulation.hpp"
15 – #endif
16 –
17 –
52   // Basic isotropic thermostating and barostating via the Melchionna
53   // modification of the Hoover algorithm:
54   //
# Line 25 | Line 59
59   //
60   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
61  
62 < template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
29 <  T( theInfo, the_ff )
30 < {
31 <  GenericData* data;
32 <  DoubleData * chiValue;
33 <  DoubleData * integralOfChidtValue;
62 > namespace oopse {
63  
64 <  chiValue = NULL;
65 <  integralOfChidtValue = NULL;
64 >  NPT::NPT(SimInfo* info) :
65 >    VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
66  
67 <  chi = 0.0;
68 <  integralOfChidt = 0.0;
69 <  have_tau_thermostat = 0;
70 <  have_tau_barostat = 0;
71 <  have_target_temp = 0;
72 <  have_target_pressure = 0;
73 <  have_chi_tolerance = 0;
74 <  have_eta_tolerance = 0;
75 <  have_pos_iter_tolerance = 0;
67 >      Globals* simParams = info_->getSimParams();
68 >    
69 >      if (!simParams->getUseIntialExtendedSystemState()) {
70 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
71 >        currSnapshot->setChi(0.0);
72 >        currSnapshot->setIntegralOfChiDt(0.0);
73 >        currSnapshot->setEta(Mat3x3d(0.0));
74 >      }
75 >    
76 >      if (!simParams->haveTargetTemp()) {
77 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
78 >        painCave.isFatal = 1;
79 >        painCave.severity = OOPSE_ERROR;
80 >        simError();
81 >      } else {
82 >        targetTemp = simParams->getTargetTemp();
83 >      }
84  
85 <  // retrieve chi and integralOfChidt from simInfo
86 <  data = info->getProperty(CHIVALUE_ID);
87 <  if(data){
88 <    chiValue = dynamic_cast<DoubleData*>(data);
52 <  }
85 >      // We must set tauThermostat
86 >      if (!simParams->haveTauThermostat()) {
87 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
88 >                "\tintegrator, you must set tauThermostat_.\n");
89  
90 <  data = info->getProperty(INTEGRALOFCHIDT_ID);
91 <  if(data){
92 <    integralOfChidtValue = dynamic_cast<DoubleData*>(data);
93 <  }
90 >        painCave.severity = OOPSE_ERROR;
91 >        painCave.isFatal = 1;
92 >        simError();
93 >      } else {
94 >        tauThermostat = simParams->getTauThermostat();
95 >      }
96  
97 <  // chi and integralOfChidt should appear by pair
98 <  if(chiValue && integralOfChidtValue){
99 <    chi = chiValue->getData();
62 <    integralOfChidt = integralOfChidtValue->getData();
63 <  }
97 >      if (!simParams->haveTargetPressure()) {
98 >        sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
99 >                "   without a targetPressure!\n");
100  
101 <  oldPos = new double[3*integrableObjects.size()];
102 <  oldVel = new double[3*integrableObjects.size()];
103 <  oldJi = new double[3*integrableObjects.size()];
101 >        painCave.isFatal = 1;
102 >        simError();
103 >      } else {
104 >        targetPressure = simParams->getTargetPressure();
105 >      }
106 >    
107 >      if (!simParams->haveTauBarostat()) {
108 >        sprintf(painCave.errMsg,
109 >                "If you use the NPT integrator, you must set tauBarostat.\n");
110 >        painCave.severity = OOPSE_ERROR;
111 >        painCave.isFatal = 1;
112 >        simError();
113 >      } else {
114 >        tauBarostat = simParams->getTauBarostat();
115 >      }
116 >    
117 >      tt2 = tauThermostat * tauThermostat;
118 >      tb2 = tauBarostat * tauBarostat;
119  
120 < }
120 >      update();
121 >    }
122  
123 < template<typename T> NPT<T>::~NPT() {
124 <  delete[] oldPos;
73 <  delete[] oldVel;
74 <  delete[] oldJi;
75 < }
123 >  NPT::~NPT() {
124 >  }
125  
126 < template<typename T> void NPT<T>::moveA() {
126 >  void NPT::doUpdate() {
127  
128 <  //new version of NPT
129 <  int i, j, k;
130 <  double Tb[3], ji[3];
82 <  double mass;
83 <  double vel[3], pos[3], frc[3];
84 <  double sc[3];
85 <  double COM[3];
128 >    oldPos.resize(info_->getNIntegrableObjects());
129 >    oldVel.resize(info_->getNIntegrableObjects());
130 >    oldJi.resize(info_->getNIntegrableObjects());
131  
132 <  instaTemp = tStats->getTemperature();
88 <  tStats->getPressureTensor( press );
89 <  instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0;
90 <  instaVol = tStats->getVolume();
132 >  }
133  
134 <  tStats->getCOM(COM);
134 >  void NPT::moveA() {
135 >    SimInfo::MoleculeIterator i;
136 >    Molecule::IntegrableObjectIterator  j;
137 >    Molecule* mol;
138 >    StuntDouble* integrableObject;
139 >    Vector3d Tb, ji;
140 >    double mass;
141 >    Vector3d vel;
142 >    Vector3d pos;
143 >    Vector3d frc;
144 >    Vector3d sc;
145 >    int index;
146  
147 <  //evolve velocity half step
147 >    chi= currentSnapshot_->getChi();
148 >    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
149 >    loadEta();
150 >    
151 >    instaTemp =thermo.getTemperature();
152 >    press = thermo.getPressureTensor();
153 >    instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
154 >    instaVol =thermo.getVolume();
155  
156 <  calcVelScale();
97 <  
98 <  for( i=0; i<integrableObjects.size(); i++ ){
156 >    Vector3d  COM = info_->getCom();
157  
158 <    integrableObjects[i]->getVel( vel );
101 <    integrableObjects[i]->getFrc( frc );
158 >    //evolve velocity half step
159  
160 <    mass = integrableObjects[i]->getMass();
160 >    calcVelScale();
161  
162 <    getVelScaleA( sc, vel );
162 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
163 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
164 >           integrableObject = mol->nextIntegrableObject(j)) {
165 >                
166 >        vel = integrableObject->getVel();
167 >        frc = integrableObject->getFrc();
168  
169 <    for (j=0; j < 3; j++) {
169 >        mass = integrableObject->getMass();
170  
171 <      // velocity half step  (use chi from previous step here):
110 <      vel[j] += dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
171 >        getVelScaleA(sc, vel);
172  
173 <    }
173 >        // velocity half step  (use chi from previous step here):
174 >        //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
175 >        vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
176 >        integrableObject->setVel(vel);
177  
178 <    integrableObjects[i]->setVel( vel );
178 >        if (integrableObject->isDirectional()) {
179  
180 <    if( integrableObjects[i]->isDirectional() ){
180 >          // get and convert the torque to body frame
181  
182 <      // get and convert the torque to body frame
182 >          Tb = integrableObject->lab2Body(integrableObject->getTrq());
183  
184 <      integrableObjects[i]->getTrq( Tb );
121 <      integrableObjects[i]->lab2Body( Tb );
184 >          // get the angular momentum, and propagate a half step
185  
186 <      // get the angular momentum, and propagate a half step
186 >          ji = integrableObject->getJ();
187  
188 <      integrableObjects[i]->getJ( ji );
188 >          //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
189 >          ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
190 >                
191 >          rotAlgo->rotate(integrableObject, ji, dt);
192  
193 <      for (j=0; j < 3; j++)
194 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
193 >          integrableObject->setJ(ji);
194 >        }
195 >            
196 >      }
197 >    }
198 >    // evolve chi and eta  half step
199  
200 <      this->rotationPropagation( integrableObjects[i], ji );
200 >    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
201 >    
202 >    evolveEtaA();
203  
204 <      integrableObjects[i]->setJ( ji );
204 >    //calculate the integral of chidt
205 >    integralOfChidt += dt2 * chi;
206 >    
207 >    index = 0;
208 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
209 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
210 >           integrableObject = mol->nextIntegrableObject(j)) {
211 >        oldPos[index++] = integrableObject->getPos();            
212 >      }
213      }
214 <  }
214 >    
215 >    //the first estimation of r(t+dt) is equal to  r(t)
216  
217 <  // evolve chi and eta  half step
217 >    for(int k = 0; k < maxIterNum_; k++) {
218 >      index = 0;
219 >      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
220 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
221 >             integrableObject = mol->nextIntegrableObject(j)) {
222  
223 <  evolveChiA();
224 <  evolveEtaA();
223 >          vel = integrableObject->getVel();
224 >          pos = integrableObject->getPos();
225  
226 <  //calculate the integral of chidt
142 <  integralOfChidt += dt2*chi;
226 >          this->getPosScale(pos, COM, index, sc);
227  
228 <  //save the old positions
229 <  for(i = 0; i < integrableObjects.size(); i++){
146 <    integrableObjects[i]->getPos(pos);
147 <    for(j = 0; j < 3; j++)
148 <      oldPos[i*3 + j] = pos[j];
149 <  }
228 >          pos = oldPos[index] + dt * (vel + sc);
229 >          integrableObject->setPos(pos);    
230  
231 <  //the first estimation of r(t+dt) is equal to  r(t)
231 >          ++index;
232 >        }
233 >      }
234  
235 <  for(k = 0; k < 5; k ++){
235 >      rattle->constraintA();
236 >    }
237  
238 <    for(i =0 ; i < integrableObjects.size(); i++){
238 >    // Scale the box after all the positions have been moved:
239  
240 <      integrableObjects[i]->getVel(vel);
158 <      integrableObjects[i]->getPos(pos);
240 >    this->scaleSimBox();
241  
242 <      this->getPosScale( pos, COM, i, sc );
243 <
162 <      for(j = 0; j < 3; j++)
163 <        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
242 >    currentSnapshot_->setChi(chi);
243 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
244  
245 <      integrableObjects[i]->setPos( pos );
166 <    }
167 <    
168 <    if(nConstrained)
169 <      constrainA();
245 >    saveEta();
246    }
247  
248 +  void NPT::moveB(void) {
249 +    SimInfo::MoleculeIterator i;
250 +    Molecule::IntegrableObjectIterator  j;
251 +    Molecule* mol;
252 +    StuntDouble* integrableObject;
253 +    int index;
254 +    Vector3d Tb;
255 +    Vector3d ji;
256 +    Vector3d sc;
257 +    Vector3d vel;
258 +    Vector3d frc;
259 +    double mass;
260  
173 –  // Scale the box after all the positions have been moved:
261  
262 <  this->scaleSimBox();
263 < }
262 >    chi= currentSnapshot_->getChi();
263 >    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264 >    double oldChi  = chi;
265 >    double prevChi;
266  
267 < template<typename T> void NPT<T>::moveB( void ){
268 <
269 <  //new version of NPT
270 <  int i, j, k;
271 <  double Tb[3], ji[3], sc[3];
272 <  double vel[3], frc[3];
273 <  double mass;
274 <
275 <  // Set things up for the iteration:
276 <
277 <  for( i=0; i<integrableObjects.size(); i++ ){
278 <
190 <    integrableObjects[i]->getVel( vel );
191 <
192 <    for (j=0; j < 3; j++)
193 <      oldVel[3*i + j]  = vel[j];
194 <
195 <    if( integrableObjects[i]->isDirectional() ){
196 <
197 <      integrableObjects[i]->getJ( ji );
198 <
199 <      for (j=0; j < 3; j++)
200 <        oldJi[3*i + j] = ji[j];
201 <
267 >    loadEta();
268 >    
269 >    //save velocity and angular momentum
270 >    index = 0;
271 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
272 >      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
273 >           integrableObject = mol->nextIntegrableObject(j)) {
274 >                
275 >        oldVel[index] = integrableObject->getVel();
276 >        oldJi[index] = integrableObject->getJ();
277 >        ++index;
278 >      }
279      }
203 –  }
280  
281 <  // do the iteration:
281 >    // do the iteration:
282 >    instaVol =thermo.getVolume();
283  
284 <  instaVol = tStats->getVolume();
284 >    for(int k = 0; k < maxIterNum_; k++) {
285 >      instaTemp =thermo.getTemperature();
286 >      instaPress =thermo.getPressure();
287  
288 <  for (k=0; k < 4; k++) {
288 >      // evolve chi another half step using the temperature at t + dt/2
289 >      prevChi = chi;
290 >      chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
291  
292 <    instaTemp = tStats->getTemperature();
293 <    instaPress = tStats->getPressure();
292 >      //evolve eta
293 >      this->evolveEtaB();
294 >      this->calcVelScale();
295  
296 <    // evolve chi another half step using the temperature at t + dt/2
296 >      index = 0;
297 >      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
298 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
299 >             integrableObject = mol->nextIntegrableObject(j)) {            
300  
301 <    this->evolveChiB();
302 <    this->evolveEtaB();
218 <    this->calcVelScale();
301 >          frc = integrableObject->getFrc();
302 >          vel = integrableObject->getVel();
303  
304 <    for( i=0; i<integrableObjects.size(); i++ ){
304 >          mass = integrableObject->getMass();
305  
306 <      integrableObjects[i]->getFrc( frc );
223 <      integrableObjects[i]->getVel(vel);
306 >          getVelScaleB(sc, index);
307  
308 <      mass = integrableObjects[i]->getMass();
308 >          // velocity half step
309 >          //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
310 >          vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
311 >          integrableObject->setVel(vel);
312  
313 <      getVelScaleB( sc, i );
313 >          if (integrableObject->isDirectional()) {
314 >            // get and convert the torque to body frame
315 >            Tb = integrableObject->lab2Body(integrableObject->getTrq());
316  
317 <      // velocity half step
318 <      for (j=0; j < 3; j++)
319 <        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
317 >            //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
318 >            ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
319 >            integrableObject->setJ(ji);
320 >          }
321  
322 <      integrableObjects[i]->setVel( vel );
323 <
235 <      if( integrableObjects[i]->isDirectional() ){
236 <
237 <        // get and convert the torque to body frame
238 <
239 <        integrableObjects[i]->getTrq( Tb );
240 <        integrableObjects[i]->lab2Body( Tb );
241 <
242 <        for (j=0; j < 3; j++)
243 <          ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
244 <
245 <          integrableObjects[i]->setJ( ji );
322 >          ++index;
323 >        }
324        }
325 +        
326 +      rattle->constraintB();
327 +
328 +      if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
329 +        break;
330      }
331  
332 <    if(nConstrained)
333 <      constrainB();
332 >    //calculate integral of chidt
333 >    integralOfChidt += dt2 * chi;
334  
335 <    if ( this->chiConverged() && this->etaConverged() ) break;
336 <  }
335 >    currentSnapshot_->setChi(chi);
336 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);    
337  
338 <  //calculate integral of chida
256 <  integralOfChidt += dt2*chi;
257 <
258 <
259 < }
260 <
261 < template<typename T> void NPT<T>::resetIntegrator() {
262 <  chi = 0.0;
263 <  T::resetIntegrator();
264 < }
265 <
266 < template<typename T> void NPT<T>::evolveChiA() {
267 <  chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
268 <  oldChi = chi;
269 < }
270 <
271 < template<typename T> void NPT<T>::evolveChiB() {
272 <
273 <  prevChi = chi;
274 <  chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
275 < }
276 <
277 < template<typename T> bool NPT<T>::chiConverged() {
278 <
279 <  return ( fabs( prevChi - chi ) <= chiTolerance );
280 < }
281 <
282 < template<typename T> int NPT<T>::readyCheck() {
283 <
284 <  //check parent's readyCheck() first
285 <  if (T::readyCheck() == -1)
286 <    return -1;
287 <
288 <  // First check to see if we have a target temperature.
289 <  // Not having one is fatal.
290 <
291 <  if (!have_target_temp) {
292 <    sprintf( painCave.errMsg,
293 <             "NPT error: You can't use the NPT integrator\n"
294 <             "   without a targetTemp!\n"
295 <             );
296 <    painCave.isFatal = 1;
297 <    simError();
298 <    return -1;
338 >    saveEta();
339    }
340  
341 <  if (!have_target_pressure) {
342 <    sprintf( painCave.errMsg,
343 <             "NPT error: You can't use the NPT integrator\n"
344 <             "   without a targetPressure!\n"
305 <             );
306 <    painCave.isFatal = 1;
307 <    simError();
308 <    return -1;
341 >  void NPT::resetIntegrator(){
342 >      currentSnapshot_->setChi(0.0);
343 >      currentSnapshot_->setIntegralOfChiDt(0.0);
344 >      resetEta();
345    }
346  
311 –  // We must set tauThermostat.
347  
348 <  if (!have_tau_thermostat) {
349 <    sprintf( painCave.errMsg,
350 <             "NPT error: If you use the NPT\n"
351 <             "   integrator, you must set tauThermostat.\n");
352 <    painCave.isFatal = 1;
318 <    simError();
319 <    return -1;
320 <  }
321 <
322 <  // We must set tauBarostat.
323 <
324 <  if (!have_tau_barostat) {
325 <    sprintf( painCave.errMsg,
326 <             "If you use the NPT integrator, you must set tauBarostat.\n");
327 <    painCave.severity = OOPSE_ERROR;
328 <    painCave.isFatal = 1;
329 <    simError();
330 <    return -1;
331 <  }
332 <
333 <  if (!have_chi_tolerance) {
334 <    sprintf( painCave.errMsg,
335 <             "Setting chi tolerance to 1e-6 in NPT integrator\n");
336 <    chiTolerance = 1e-6;
337 <    have_chi_tolerance = 1;
338 <    painCave.severity = OOPSE_INFO;
339 <    painCave.isFatal = 0;
340 <    simError();
341 <  }
342 <
343 <  if (!have_eta_tolerance) {
344 <    sprintf( painCave.errMsg,
345 <             "Setting eta tolerance to 1e-6 in NPT integrator");
346 <    etaTolerance = 1e-6;
347 <    have_eta_tolerance = 1;
348 <    painCave.severity = OOPSE_INFO;
349 <    painCave.isFatal = 0;
350 <    simError();
351 <  }
352 <
353 <  // We need NkBT a lot, so just set it here: This is the RAW number
354 <  // of integrableObjects, so no subtraction or addition of constraints or
355 <  // orientational degrees of freedom:
356 <
357 <  NkBT = (double)(info->getTotIntegrableObjects()) * kB * targetTemp;
358 <
359 <  // fkBT is used because the thermostat operates on more degrees of freedom
360 <  // than the barostat (when there are particles with orientational degrees
361 <  // of freedom).  
362 <
363 <  fkBT = (double)(info->getNDF()) * kB * targetTemp;
364 <
365 <  tt2 = tauThermostat * tauThermostat;
366 <  tb2 = tauBarostat * tauBarostat;
367 <
368 <  return 1;
348 >    void NPT::resetEta() {
349 >      Mat3x3d etaMat(0.0);
350 >      currentSnapshot_->setEta(etaMat);    
351 >    }
352 >    
353   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines