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