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, 24107 (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->setChi(0.0);
        currSnapshot->setIntegralOfChiDt(0.0);
        currSnapshot->setEta(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* integrableObject;
    Vector3d Tb, ji;
    RealType mass;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d sc;
    int index;

    chi= currentSnapshot_->getChi();
    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
    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 = info_->getCom();

    //evolve velocity half step

    calcVelScale();

    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
                
        vel = integrableObject->getVel();
        frc = integrableObject->getFrc();

        mass = integrableObject->getMass();

        getVelScaleA(sc, vel);

        // velocity half step  (use chi from previous step here):
        //vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
        vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
        integrableObject->setVel(vel);

        if (integrableObject->isDirectional()) {

          // get and convert the torque to body frame

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

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

          ji = integrableObject->getJ();

          //ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
          ji += dt2*PhysicalConstants::energyConvert * Tb - dt2*chi* ji;
                
          rotAlgo_->rotate(integrableObject, ji, dt);

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

    chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
    
    evolveEtaA();

    //calculate the integral of chidt
    integralOfChidt += dt2 * chi;
    
    flucQ_->moveA();


    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
        oldPos[index++] = integrableObject->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 (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
             integrableObject = mol->nextIntegrableObject(j)) {

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

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

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

          ++index;
        }
      }

      rattle_->constraintA();
    }

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

    this->scaleSimBox();

    currentSnapshot_->setChi(chi);
    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);

    saveEta();
  }

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


    chi= currentSnapshot_->getChi();
    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
    RealType oldChi  = chi;
    RealType prevChi;

    loadEta();
    
    //save velocity and angular momentum
    index = 0;
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
           integrableObject = mol->nextIntegrableObject(j)) {
                
        oldVel[index] = integrableObject->getVel();

        if (integrableObject->isDirectional())
           oldJi[index] = integrableObject->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 = chi;
      chi = 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 (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
             integrableObject = mol->nextIntegrableObject(j)) {            

          frc = integrableObject->getFrc();
          vel = integrableObject->getVel();

          mass = integrableObject->getMass();

          getVelScaleB(sc, index);

          // velocity half step
          //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
          vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
          integrableObject->setVel(vel);

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

            //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
            ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
            integrableObject->setJ(ji);
          }

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

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

    //calculate integral of chidt
    integralOfChidt += dt2 * chi;

    currentSnapshot_->setChi(chi);
    currentSnapshot_->setIntegralOfChiDt(integralOfChidt);    

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

  void NPT::resetIntegrator(){
      currentSnapshot_->setChi(0.0);
      currentSnapshot_->setIntegralOfChiDt(0.0);
      resetEta();
  }


    void NPT::resetEta() {
      Mat3x3d etaMat(0.0);
      currentSnapshot_->setEta(etaMat);    
    }
    
}
Revision:	1759
Committed:	Thu Jun 21 18:55:33 2012 UTC (12 years, 10 months ago) by gezelter
File size:	10975 byte(s)
Log Message:	fixed a directional StuntDouble bug
#	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			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	gezelter	1665	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40			* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	gezelter	246	*/
42
43	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			currSnapshot->setChi(0.0);
73			currSnapshot->setIntegralOfChiDt(0.0);
74			currSnapshot->setEta(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	1715	updateSizes();
122	gezelter	507	}
123	gezelter	2
124	gezelter	507	NPT::~NPT() {
125			}
126	gezelter	2
127	gezelter	1715	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			StuntDouble* integrableObject;
140			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	246	chi= currentSnapshot_->getChi();
149			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
150			loadEta();
151
152			instaTemp =thermo.getTemperature();
153			press = thermo.getPressureTensor();
154	gezelter	1390	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
155	gezelter	246	instaVol =thermo.getVolume();
156	gezelter	2
157	gezelter	246	Vector3d COM = info_->getCom();
158	gezelter	2
159	gezelter	246	//evolve velocity half step
160	gezelter	2
161	gezelter	246	calcVelScale();
162	gezelter	2
163	gezelter	246	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164	gezelter	507	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165			integrableObject = mol->nextIntegrableObject(j)) {
166	gezelter	246
167	gezelter	507	vel = integrableObject->getVel();
168			frc = integrableObject->getFrc();
169	gezelter	2
170	gezelter	507	mass = integrableObject->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	1390	//vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
176			vel += dt2PhysicalConstants::energyConvert/mass frc - dt2*sc;
177	gezelter	507	integrableObject->setVel(vel);
178	gezelter	2
179	gezelter	507	if (integrableObject->isDirectional()) {
180	gezelter	2
181	gezelter	507	// get and convert the torque to body frame
182	gezelter	2
183	gezelter	507	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184	gezelter	2
185	gezelter	507	// get the angular momentum, and propagate a half step
186	gezelter	2
187	gezelter	507	ji = integrableObject->getJ();
188	gezelter	2
189	gezelter	1390	//ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
190			ji += dt2PhysicalConstants::energyConvert Tb - dt2chi ji;
191	gezelter	246
192	gezelter	1715	rotAlgo_->rotate(integrableObject, ji, dt);
193	gezelter	2
194	gezelter	507	integrableObject->setJ(ji);
195			}
196	gezelter	246
197	gezelter	507	}
198	gezelter	246	}
199			// evolve chi and eta half step
200	gezelter	2
201	gezelter	246	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202
203			evolveEtaA();
204	gezelter	2
205	gezelter	246	//calculate the integral of chidt
206			integralOfChidt += dt2 * chi;
207
208	gezelter	1715	flucQ_->moveA();
209
210
211	gezelter	246	index = 0;
212			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
213	gezelter	507	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
214			integrableObject = mol->nextIntegrableObject(j)) {
215			oldPos[index++] = integrableObject->getPos();
216			}
217	gezelter	2	}
218	gezelter	246
219			//the first estimation of r(t+dt) is equal to r(t)
220	gezelter	2
221	gezelter	246	for(int k = 0; k < maxIterNum_; k++) {
222	gezelter	507	index = 0;
223			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
224			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
225			integrableObject = mol->nextIntegrableObject(j)) {
226	gezelter	2
227	gezelter	507	vel = integrableObject->getVel();
228			pos = integrableObject->getPos();
229	gezelter	2
230	gezelter	507	this->getPosScale(pos, COM, index, sc);
231	gezelter	2
232	gezelter	507	pos = oldPos[index] + dt * (vel + sc);
233			integrableObject->setPos(pos);
234	gezelter	2
235	gezelter	507	++index;
236			}
237			}
238	gezelter	2
239	gezelter	1715	rattle_->constraintA();
240	gezelter	246	}
241	gezelter	2
242	gezelter	246	// Scale the box after all the positions have been moved:
243	gezelter	2
244	gezelter	246	this->scaleSimBox();
245	gezelter	2
246	gezelter	246	currentSnapshot_->setChi(chi);
247			currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
248	gezelter	2
249	gezelter	246	saveEta();
250	gezelter	507	}
251	gezelter	2
252	gezelter	507	void NPT::moveB(void) {
253	gezelter	246	SimInfo::MoleculeIterator i;
254			Molecule::IntegrableObjectIterator j;
255			Molecule* mol;
256			StuntDouble* integrableObject;
257			int index;
258			Vector3d Tb;
259			Vector3d ji;
260			Vector3d sc;
261			Vector3d vel;
262			Vector3d frc;
263	tim	963	RealType mass;
264	gezelter	2
265
266	gezelter	246	chi= currentSnapshot_->getChi();
267			integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
268	tim	963	RealType oldChi = chi;
269			RealType prevChi;
270	gezelter	2
271	gezelter	246	loadEta();
272
273			//save velocity and angular momentum
274			index = 0;
275			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
276	gezelter	507	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
277			integrableObject = mol->nextIntegrableObject(j)) {
278	gezelter	246
279	gezelter	507	oldVel[index] = integrableObject->getVel();
280	gezelter	1759
281			if (integrableObject->isDirectional())
282			oldJi[index] = integrableObject->getJ();
283
284	gezelter	507	++index;
285			}
286	gezelter	2	}
287
288	gezelter	246	// do the iteration:
289			instaVol =thermo.getVolume();
290	gezelter	2
291	gezelter	246	for(int k = 0; k < maxIterNum_; k++) {
292	gezelter	507	instaTemp =thermo.getTemperature();
293			instaPress =thermo.getPressure();
294	gezelter	2
295	gezelter	507	// evolve chi another half step using the temperature at t + dt/2
296			prevChi = chi;
297			chi = oldChi + dt2 * (instaTemp / targetTemp - 1.0) / tt2;
298	gezelter	2
299	gezelter	507	//evolve eta
300			this->evolveEtaB();
301			this->calcVelScale();
302	gezelter	2
303	gezelter	507	index = 0;
304			for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
305			for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
306			integrableObject = mol->nextIntegrableObject(j)) {
307	gezelter	2
308	gezelter	507	frc = integrableObject->getFrc();
309			vel = integrableObject->getVel();
310	gezelter	2
311	gezelter	507	mass = integrableObject->getMass();
312	gezelter	2
313	gezelter	507	getVelScaleB(sc, index);
314	gezelter	2
315	gezelter	507	// velocity half step
316	gezelter	1390	//vel[j] = oldVel[3 * i + j] + dt2 ((frc[j] / mass) PhysicalConstants::energyConvert - sc[j]);
317			vel = oldVel[index] + dt2PhysicalConstants::energyConvert/mass frc - dt2*sc;
318	gezelter	507	integrableObject->setVel(vel);
319	gezelter	2
320	gezelter	507	if (integrableObject->isDirectional()) {
321			// get and convert the torque to body frame
322			Tb = integrableObject->lab2Body(integrableObject->getTrq());
323	gezelter	2
324	gezelter	1390	//ji[j] = oldJi[3i + j] + dt2 (Tb[j] * PhysicalConstants::energyConvert - oldJi[3i+j]chi);
325			ji = oldJi[index] + dt2PhysicalConstants::energyConvertTb - dt2chioldJi[index];
326	gezelter	507	integrableObject->setJ(ji);
327			}
328	gezelter	2
329	gezelter	507	++index;
330			}
331			}
332	gezelter	246
333	gezelter	1715	rattle_->constraintB();
334	gezelter	2
335	gezelter	507	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
336			break;
337	gezelter	2	}
338
339	gezelter	246	//calculate integral of chidt
340			integralOfChidt += dt2 * chi;
341	gezelter	2
342	gezelter	246	currentSnapshot_->setChi(chi);
343			currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
344	gezelter	2
345	gezelter	1715	flucQ_->moveB();
346	gezelter	246	saveEta();
347	gezelter	507	}
348	gezelter	2
349	tim	546	void NPT::resetIntegrator(){
350			currentSnapshot_->setChi(0.0);
351			currentSnapshot_->setIntegralOfChiDt(0.0);
352			resetEta();
353			}
354
355
356			void NPT::resetEta() {
357			Mat3x3d etaMat(0.0);
358			currentSnapshot_->setEta(etaMat);
359			}
360
361	gezelter	2	}