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

Comparing trunk/OOPSE-2.0/src/integrators/NPT.cpp (file contents):
Revision 1625 by tim, Thu Oct 21 16:22:01 2004 UTC vs.
Revision 1930 by gezelter, Wed Jan 12 22:41:40 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 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
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 <  DoubleGenericData * chiValue;
33 <  DoubleGenericData * 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->getUseInitXSstate()) {
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<DoubleGenericData*>(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<DoubleGenericData*>(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 +    update();
121   }
122  
123 < template<typename T> NPT<T>::~NPT() {
72 <  delete[] oldPos;
73 <  delete[] oldVel;
74 <  delete[] oldJi;
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 );
242 >    currentSnapshot_->setChi(chi);
243 >    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
244  
245 <      for(j = 0; j < 3; j++)
163 <        pos[j] = oldPos[i*3 + j] + dt*(vel[j] + sc[j]);
164 <
165 <      integrableObjects[i]->setPos( pos );
166 <    }
167 <    
168 <    if(nConstrained)
169 <      constrainA();
170 <  }
171 <
172 <
173 <  // Scale the box after all the positions have been moved:
174 <
175 <  this->scaleSimBox();
245 >    saveEta();
246   }
247  
248 < template<typename T> void NPT<T>::moveB( void ){
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  
180 –  //new version of NPT
181 –  int i, j, k;
182 –  double Tb[3], ji[3], sc[3];
183 –  double vel[3], frc[3];
184 –  double mass;
261  
262 <  // Set things up for the iteration:
262 >    chi= currentSnapshot_->getChi();
263 >    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
264 >    double oldChi  = chi;
265 >    double prevChi;
266  
267 <  for( i=0; i<integrableObjects.size(); i++ ){
268 <
269 <    integrableObjects[i]->getVel( vel );
270 <
271 <    for (j=0; j < 3; j++)
272 <      oldVel[3*i + j]  = vel[j];
273 <
274 <    if( integrableObjects[i]->isDirectional() ){
275 <
276 <      integrableObjects[i]->getJ( ji );
277 <
278 <      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 );
322 >                ++index;
323 >            }
324 >        }
325 >        
326 >        rattle->constraintB();
327  
328 <      if( integrableObjects[i]->isDirectional() ){
329 <
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 );
246 <      }
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 <
338 >    saveEta();
339   }
340  
261 – template<typename T> void NPT<T>::resetIntegrator() {
262 –  chi = 0.0;
263 –  T::resetIntegrator();
341   }
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;
299 –  }
300 –
301 –  if (!have_target_pressure) {
302 –    sprintf( painCave.errMsg,
303 –             "NPT error: You can't use the NPT integrator\n"
304 –             "   without a targetPressure!\n"
305 –             );
306 –    painCave.isFatal = 1;
307 –    simError();
308 –    return -1;
309 –  }
310 –
311 –  // We must set tauThermostat.
312 –
313 –  if (!have_tau_thermostat) {
314 –    sprintf( painCave.errMsg,
315 –             "NPT error: If you use the NPT\n"
316 –             "   integrator, you must set tauThermostat.\n");
317 –    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;
369 – }

Diff Legend

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