src/integrators/NPT.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 <math.h>

#include "brains/SimInfo.hpp"
#include "brains/Thermo.hpp"
#include "integrators/NPT.hpp"
#include "math/SquareMatrix3.hpp"
#include "primitives/Molecule.hpp"
#include "utils/PhysicalConstants.hpp"
#include "utils/simError.h"

// Basic 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.

namespace OpenMD {

  NPT::NPT(SimInfo* info) :
    VelocityVerletIntegrator(info), etaTolerance(1e-6), chiTolerance(1e-6), 
    maxIterNum_(4) {

      Globals* simParams = info_->getSimParams();
    
      if (!simParams->getUseIntialExtendedSystemState()) {
        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
        currSnapshot->setThermostat(make_pair(0.0, 0.0));
        currSnapshot->setBarostat(Mat3x3d(0.0));
      }
    
      if (!simParams->haveTargetTemp()) {
        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
        painCave.isFatal = 1;
        painCave.severity = OPENMD_ERROR;
        simError();
      } else {
        targetTemp = simParams->getTargetTemp();
      }

      // We must set tauThermostat
      if (!simParams->haveTauThermostat()) {
        sprintf(painCave.errMsg, "If you use the constant temperature\n"
                "\tintegrator, you must set tauThermostat.\n");

        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();
      } else {
        tauThermostat = simParams->getTauThermostat();
      }

      if (!simParams->haveTargetPressure()) {
        sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
                "   without a targetPressure!\n");

        painCave.isFatal = 1;
        simError();
      } else {
        targetPressure = simParams->getTargetPressure();
      }
    
      if (!simParams->haveTauBarostat()) {
        sprintf(painCave.errMsg,
                "If you use the NPT integrator, you must set tauBarostat.\n");
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();
      } else {
        tauBarostat = simParams->getTauBarostat();
      }
    
      tt2 = tauThermostat * tauThermostat;
      tb2 = tauBarostat * tauBarostat;

      updateSizes();
    }

  NPT::~NPT() {
  }

  void NPT::doUpdateSizes() {

    oldPos.resize(info_->getNIntegrableObjects());
    oldVel.resize(info_->getNIntegrableObjects());
    oldJi.resize(info_->getNIntegrableObjects());

  }

  void NPT::moveA() {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* sd;
    Vector3d Tb, ji;
    RealType mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d sc;
    int index;

    thermostat = snap->getThermostat();
    loadEta();
    
    instaTemp =thermo.getTemperature();
    press = thermo.getPressureTensor();
    instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
    instaVol =thermo.getVolume();

    Vector3d  COM = thermo.getCom();

    //evolve velocity half step

    calcVelScale();

    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
           sd = mol->nextIntegrableObject(j)) {
                
        vel = sd->getVel();
        frc = sd->getFrc();

        mass = sd->getMass();

        getVelScaleA(sc, vel);

        // velocity half step  (use chi from previous step here):

        vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
        sd->setVel(vel);

        if (sd->isDirectional()) {

          // get and convert the torque to body frame

          Tb = sd->lab2Body(sd->getTrq());

          // get the angular momentum, and propagate a half step

          ji = sd->getJ();

          ji += dt2*PhysicalConstants::energyConvert * Tb 
            - dt2*thermostat.first* ji;
                
          rotAlgo_->rotate(sd, ji, dt);

          sd->setJ(ji);
        }
            
      }
    }
    // evolve chi and eta  half step

    thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
    
    evolveEtaA();

    //calculate the integral of chidt
    thermostat.second += dt2 * thermostat.first;
    
    flucQ_->moveA();


    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
           sd = mol->nextIntegrableObject(j)) {

        oldPos[index++] = sd->getPos();            

      }
    }
    
    //the first estimation of r(t+dt) is equal to  r(t)

    for(int k = 0; k < maxIterNum_; k++) {
      index = 0;
      for (mol = info_->beginMolecule(i); mol != NULL; 
           mol = info_->nextMolecule(i)) {

        for (sd = mol->beginIntegrableObject(j); sd != NULL;
             sd = mol->nextIntegrableObject(j)) {

          vel = sd->getVel();
          pos = sd->getPos();

          this->getPosScale(pos, COM, index, sc);

          pos = oldPos[index] + dt * (vel + sc);
          sd->setPos(pos);     

          ++index;
        }
      }

      rattle_->constraintA();
    }

    // Scale the box after all the positions have been moved:

    this->scaleSimBox();

    snap->setThermostat(thermostat);

    saveEta();
  }

  void NPT::moveB(void) {
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* sd;
    int index;
    Vector3d Tb;
    Vector3d ji;
    Vector3d sc;
    Vector3d vel;
    Vector3d frc;
    RealType mass;

    thermostat = snap->getThermostat();
    RealType oldChi  = thermostat.first;
    RealType prevChi;

    loadEta();
    
    //save velocity and angular momentum
    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; 
         mol = info_->nextMolecule(i)) {

      for (sd = mol->beginIntegrableObject(j); sd != NULL;
           sd = mol->nextIntegrableObject(j)) {
                
        oldVel[index] = sd->getVel();

        if (sd->isDirectional())
           oldJi[index] = sd->getJ();

        ++index;
      }
    }

    // do the iteration:
    instaVol =thermo.getVolume();

    for(int k = 0; k < maxIterNum_; k++) {
      instaTemp =thermo.getTemperature();
      instaPress =thermo.getPressure();

      // evolve chi another half step using the temperature at t + dt/2
      prevChi = thermostat.first;
      thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;

      //evolve eta
      this->evolveEtaB();
      this->calcVelScale();

      index = 0;
      for (mol = info_->beginMolecule(i); mol != NULL; 
           mol = info_->nextMolecule(i)) {

        for (sd = mol->beginIntegrableObject(j); sd != NULL;
             sd = mol->nextIntegrableObject(j)) {            

          frc = sd->getFrc();
          mass = sd->getMass();

          getVelScaleB(sc, index);

          // velocity half step
          vel = oldVel[index] 
            + dt2*PhysicalConstants::energyConvert/mass* frc 
            - dt2*sc;

          sd->setVel(vel);

          if (sd->isDirectional()) {
            // get and convert the torque to body frame
            Tb = sd->lab2Body(sd->getTrq());

            ji = oldJi[index] 
              + dt2*PhysicalConstants::energyConvert*Tb 
              - dt2*thermostat.first*oldJi[index];

            sd->setJ(ji);
          }

          ++index;
        }
      }
        
      rattle_->constraintB();

      if ((fabs(prevChi - thermostat.first) <= chiTolerance) && 
          this->etaConverged())
        break;
    }

    //calculate integral of chidt
    thermostat.second += dt2 * thermostat.first;

    snap->setThermostat(thermostat);

    flucQ_->moveB();
    saveEta();
  }

  void NPT::resetIntegrator(){
    snap->setThermostat(make_pair(0.0, 0.0));
    resetEta();
  }

  void NPT::resetEta() {
    Mat3x3d etaMat(0.0);
    snap->setBarostat(etaMat);    
  }
}
Revision:	2071
Committed:	Sat Mar 7 21:41:51 2015 UTC (10 years, 2 months ago) by gezelter
File size:	9956 byte(s)
Log Message:	Reducing the number of warnings when using g++ to compile.
#	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	1879	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	gezelter	1782	* [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	gezelter	2	#include <math.h>
44
45	tim	3	#include "brains/SimInfo.hpp"
46			#include "brains/Thermo.hpp"
47	gezelter	246	#include "integrators/NPT.hpp"
48			#include "math/SquareMatrix3.hpp"
49			#include "primitives/Molecule.hpp"
50	gezelter	1390	#include "utils/PhysicalConstants.hpp"
51	tim	3	#include "utils/simError.h"
52	gezelter	2
53			// Basic 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
63	gezelter	1390	namespace OpenMD {
64	gezelter	2
65	gezelter	507	NPT::NPT(SimInfo* info) :
66	gezelter	2071	VelocityVerletIntegrator(info), etaTolerance(1e-6), chiTolerance(1e-6),
67			maxIterNum_(4) {
68	gezelter	2
69	gezelter	507	Globals* simParams = info_->getSimParams();
70	gezelter	246
71	tim	665	if (!simParams->getUseIntialExtendedSystemState()) {
72	gezelter	246	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
73	gezelter	1782	currSnapshot->setThermostat(make_pair(0.0, 0.0));
74			currSnapshot->setBarostat(Mat3x3d(0.0));
75	gezelter	507	}
76	gezelter	246
77	gezelter	507	if (!simParams->haveTargetTemp()) {
78	gezelter	246	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
79			painCave.isFatal = 1;
80	gezelter	1390	painCave.severity = OPENMD_ERROR;
81	gezelter	246	simError();
82	gezelter	507	} else {
83	gezelter	246	targetTemp = simParams->getTargetTemp();
84	gezelter	507	}
85	gezelter	2
86	gezelter	507	// We must set tauThermostat
87			if (!simParams->haveTauThermostat()) {
88	gezelter	246	sprintf(painCave.errMsg, "If you use the constant temperature\n"
89	gezelter	1277	"\tintegrator, you must set tauThermostat.\n");
90	gezelter	2
91	gezelter	1390	painCave.severity = OPENMD_ERROR;
92	gezelter	246	painCave.isFatal = 1;
93			simError();
94	gezelter	507	} else {
95	gezelter	246	tauThermostat = simParams->getTauThermostat();
96	gezelter	507	}
97	gezelter	2
98	gezelter	507	if (!simParams->haveTargetPressure()) {
99	gezelter	246	sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
100	gezelter	507	" without a targetPressure!\n");
101	gezelter	2
102	gezelter	246	painCave.isFatal = 1;
103			simError();
104	gezelter	507	} else {
105	gezelter	246	targetPressure = simParams->getTargetPressure();
106	gezelter	507	}
107	gezelter	246
108	gezelter	507	if (!simParams->haveTauBarostat()) {
109	gezelter	246	sprintf(painCave.errMsg,
110			"If you use the NPT integrator, you must set tauBarostat.\n");
111	gezelter	1390	painCave.severity = OPENMD_ERROR;
112	gezelter	246	painCave.isFatal = 1;
113			simError();
114	gezelter	507	} else {
115	gezelter	246	tauBarostat = simParams->getTauBarostat();
116	gezelter	507	}
117	gezelter	246
118	gezelter	507	tt2 = tauThermostat * tauThermostat;
119			tb2 = tauBarostat * tauBarostat;
120	gezelter	2
121	gezelter	1782	updateSizes();
122	gezelter	507	}
123	gezelter	2
124	gezelter	507	NPT::~NPT() {
125			}
126	gezelter	2
127	gezelter	1782	void NPT::doUpdateSizes() {
128	gezelter	2
129	gezelter	246	oldPos.resize(info_->getNIntegrableObjects());
130			oldVel.resize(info_->getNIntegrableObjects());
131			oldJi.resize(info_->getNIntegrableObjects());
132	gezelter	2
133	gezelter	507	}
134	gezelter	2
135	gezelter	507	void NPT::moveA() {
136	gezelter	246	SimInfo::MoleculeIterator i;
137			Molecule::IntegrableObjectIterator j;
138			Molecule* mol;
139	gezelter	1782	StuntDouble* sd;
140	gezelter	246	Vector3d Tb, ji;
141	tim	963	RealType mass;
142	gezelter	246	Vector3d vel;
143			Vector3d pos;
144			Vector3d frc;
145			Vector3d sc;
146			int index;
147	gezelter	2
148	gezelter	1782	thermostat = snap->getThermostat();
149	gezelter	246	loadEta();
150
151			instaTemp =thermo.getTemperature();
152			press = thermo.getPressureTensor();
153	gezelter	1390	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
154	gezelter	246	instaVol =thermo.getVolume();
155	gezelter	2
156	gezelter	1782	Vector3d COM = thermo.getCom();
157	gezelter	2
158	gezelter	246	//evolve velocity half step
159	gezelter	2
160	gezelter	246	calcVelScale();
161	gezelter	2
162	gezelter	1782	for (mol = info_->beginMolecule(i); mol != NULL;
163			mol = info_->nextMolecule(i)) {
164
165			for (sd = mol->beginIntegrableObject(j); sd != NULL;
166			sd = mol->nextIntegrableObject(j)) {
167	gezelter	246
168	gezelter	1782	vel = sd->getVel();
169			frc = sd->getFrc();
170	gezelter	2
171	gezelter	1782	mass = sd->getMass();
172	gezelter	2
173	gezelter	507	getVelScaleA(sc, vel);
174	gezelter	2
175	gezelter	507	// velocity half step (use chi from previous step here):
176	gezelter	1782
177	gezelter	1390	vel += dt2PhysicalConstants::energyConvert/mass frc - dt2*sc;
178	gezelter	1782	sd->setVel(vel);
179	gezelter	2
180	gezelter	1782	if (sd->isDirectional()) {
181	gezelter	2
182	gezelter	507	// get and convert the torque to body frame
183	gezelter	2
184	gezelter	1782	Tb = sd->lab2Body(sd->getTrq());
185	gezelter	2
186	gezelter	507	// get the angular momentum, and propagate a half step
187	gezelter	2
188	gezelter	1782	ji = sd->getJ();
189	gezelter	2
190	gezelter	1782	ji += dt2PhysicalConstants::energyConvert Tb
191			- dt2thermostat.first ji;
192	gezelter	246
193	gezelter	1782	rotAlgo_->rotate(sd, ji, dt);
194	gezelter	2
195	gezelter	1782	sd->setJ(ji);
196	gezelter	507	}
197	gezelter	246
198	gezelter	507	}
199	gezelter	246	}
200			// evolve chi and eta half step
201	gezelter	2
202	gezelter	1782	thermostat.first += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
203	gezelter	246
204			evolveEtaA();
205	gezelter	2
206	gezelter	246	//calculate the integral of chidt
207	gezelter	1782	thermostat.second += dt2 * thermostat.first;
208	gezelter	246
209	gezelter	1782	flucQ_->moveA();
210
211
212	gezelter	246	index = 0;
213	gezelter	1782	for (mol = info_->beginMolecule(i); mol != NULL;
214			mol = info_->nextMolecule(i)) {
215
216			for (sd = mol->beginIntegrableObject(j); sd != NULL;
217			sd = mol->nextIntegrableObject(j)) {
218
219			oldPos[index++] = sd->getPos();
220
221	gezelter	507	}
222	gezelter	2	}
223	gezelter	246
224			//the first estimation of r(t+dt) is equal to r(t)
225	gezelter	2
226	gezelter	246	for(int k = 0; k < maxIterNum_; k++) {
227	gezelter	507	index = 0;
228	gezelter	1782	for (mol = info_->beginMolecule(i); mol != NULL;
229			mol = info_->nextMolecule(i)) {
230	gezelter	2
231	gezelter	1782	for (sd = mol->beginIntegrableObject(j); sd != NULL;
232			sd = mol->nextIntegrableObject(j)) {
233	gezelter	2
234	gezelter	1782	vel = sd->getVel();
235			pos = sd->getPos();
236
237	gezelter	507	this->getPosScale(pos, COM, index, sc);
238	gezelter	2
239	gezelter	507	pos = oldPos[index] + dt * (vel + sc);
240	gezelter	1782	sd->setPos(pos);
241	gezelter	2
242	gezelter	507	++index;
243			}
244			}
245	gezelter	2
246	gezelter	1782	rattle_->constraintA();
247	gezelter	246	}
248	gezelter	2
249	gezelter	246	// Scale the box after all the positions have been moved:
250	gezelter	2
251	gezelter	246	this->scaleSimBox();
252	gezelter	2
253	gezelter	1782	snap->setThermostat(thermostat);
254	gezelter	2
255	gezelter	246	saveEta();
256	gezelter	507	}
257	gezelter	2
258	gezelter	507	void NPT::moveB(void) {
259	gezelter	246	SimInfo::MoleculeIterator i;
260			Molecule::IntegrableObjectIterator j;
261			Molecule* mol;
262	gezelter	1782	StuntDouble* sd;
263	gezelter	246	int index;
264			Vector3d Tb;
265			Vector3d ji;
266			Vector3d sc;
267			Vector3d vel;
268			Vector3d frc;
269	tim	963	RealType mass;
270	gezelter	2
271	gezelter	1782	thermostat = snap->getThermostat();
272			RealType oldChi = thermostat.first;
273	tim	963	RealType prevChi;
274	gezelter	2
275	gezelter	246	loadEta();
276
277			//save velocity and angular momentum
278			index = 0;
279	gezelter	1782	for (mol = info_->beginMolecule(i); mol != NULL;
280			mol = info_->nextMolecule(i)) {
281
282			for (sd = mol->beginIntegrableObject(j); sd != NULL;
283			sd = mol->nextIntegrableObject(j)) {
284	gezelter	246
285	gezelter	1782	oldVel[index] = sd->getVel();
286
287			if (sd->isDirectional())
288			oldJi[index] = sd->getJ();
289
290	gezelter	507	++index;
291			}
292	gezelter	2	}
293
294	gezelter	246	// do the iteration:
295			instaVol =thermo.getVolume();
296	gezelter	2
297	gezelter	246	for(int k = 0; k < maxIterNum_; k++) {
298	gezelter	507	instaTemp =thermo.getTemperature();
299			instaPress =thermo.getPressure();
300	gezelter	2
301	gezelter	507	// evolve chi another half step using the temperature at t + dt/2
302	gezelter	1782	prevChi = thermostat.first;
303			thermostat.first = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
304	gezelter	2
305	gezelter	507	//evolve eta
306			this->evolveEtaB();
307			this->calcVelScale();
308	gezelter	2
309	gezelter	507	index = 0;
310	gezelter	1782	for (mol = info_->beginMolecule(i); mol != NULL;
311			mol = info_->nextMolecule(i)) {
312	gezelter	2
313	gezelter	1782	for (sd = mol->beginIntegrableObject(j); sd != NULL;
314			sd = mol->nextIntegrableObject(j)) {
315	gezelter	2
316	gezelter	1782	frc = sd->getFrc();
317			mass = sd->getMass();
318
319	gezelter	507	getVelScaleB(sc, index);
320	gezelter	2
321	gezelter	507	// velocity half step
322	gezelter	1782	vel = oldVel[index]
323			+ dt2PhysicalConstants::energyConvert/mass frc
324			- dt2*sc;
325	gezelter	2
326	gezelter	1782	sd->setVel(vel);
327
328			if (sd->isDirectional()) {
329	gezelter	507	// get and convert the torque to body frame
330	gezelter	1782	Tb = sd->lab2Body(sd->getTrq());
331	gezelter	2
332	gezelter	1782	ji = oldJi[index]
333			+ dt2PhysicalConstants::energyConvertTb
334			- dt2thermostat.firstoldJi[index];
335
336			sd->setJ(ji);
337	gezelter	507	}
338	gezelter	2
339	gezelter	507	++index;
340			}
341			}
342	gezelter	246
343	gezelter	1782	rattle_->constraintB();
344	gezelter	2
345	gezelter	1782	if ((fabs(prevChi - thermostat.first) <= chiTolerance) &&
346			this->etaConverged())
347	gezelter	507	break;
348	gezelter	2	}
349
350	gezelter	246	//calculate integral of chidt
351	gezelter	1782	thermostat.second += dt2 * thermostat.first;
352	gezelter	2
353	gezelter	1782	snap->setThermostat(thermostat);
354	gezelter	2
355	gezelter	1782	flucQ_->moveB();
356	gezelter	246	saveEta();
357	gezelter	507	}
358	gezelter	2
359	tim	546	void NPT::resetIntegrator(){
360	gezelter	1782	snap->setThermostat(make_pair(0.0, 0.0));
361			resetEta();
362	tim	546	}
363
364	gezelter	1782	void NPT::resetEta() {
365			Mat3x3d etaMat(0.0);
366			snap->setBarostat(etaMat);
367			}
368	gezelter	2	}