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
#	Content
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. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* 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	*
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	* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40	* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	*/
42
43	#include <math.h>
44
45	#include "brains/SimInfo.hpp"
46	#include "brains/Thermo.hpp"
47	#include "integrators/NPT.hpp"
48	#include "math/SquareMatrix3.hpp"
49	#include "primitives/Molecule.hpp"
50	#include "utils/PhysicalConstants.hpp"
51	#include "utils/simError.h"
52
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	namespace OpenMD {
64
65	NPT::NPT(SimInfo* info) :
66	VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
67
68	Globals* simParams = info_->getSimParams();
69
70	if (!simParams->getUseIntialExtendedSystemState()) {
71	Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
72	currSnapshot->setChi(0.0);
73	currSnapshot->setIntegralOfChiDt(0.0);
74	currSnapshot->setEta(Mat3x3d(0.0));
75	}
76
77	if (!simParams->haveTargetTemp()) {
78	sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
79	painCave.isFatal = 1;
80	painCave.severity = OPENMD_ERROR;
81	simError();
82	} else {
83	targetTemp = simParams->getTargetTemp();
84	}
85
86	// We must set tauThermostat
87	if (!simParams->haveTauThermostat()) {
88	sprintf(painCave.errMsg, "If you use the constant temperature\n"
89	"\tintegrator, you must set tauThermostat.\n");
90
91	painCave.severity = OPENMD_ERROR;
92	painCave.isFatal = 1;
93	simError();
94	} else {
95	tauThermostat = simParams->getTauThermostat();
96	}
97
98	if (!simParams->haveTargetPressure()) {
99	sprintf(painCave.errMsg, "NPT error: You can't use the NPT integrator\n"
100	" without a targetPressure!\n");
101
102	painCave.isFatal = 1;
103	simError();
104	} else {
105	targetPressure = simParams->getTargetPressure();
106	}
107
108	if (!simParams->haveTauBarostat()) {
109	sprintf(painCave.errMsg,
110	"If you use the NPT integrator, you must set tauBarostat.\n");
111	painCave.severity = OPENMD_ERROR;
112	painCave.isFatal = 1;
113	simError();
114	} else {
115	tauBarostat = simParams->getTauBarostat();
116	}
117
118	tt2 = tauThermostat * tauThermostat;
119	tb2 = tauBarostat * tauBarostat;
120
121	updateSizes();
122	}
123
124	NPT::~NPT() {
125	}
126
127	void NPT::doUpdateSizes() {
128
129	oldPos.resize(info_->getNIntegrableObjects());
130	oldVel.resize(info_->getNIntegrableObjects());
131	oldJi.resize(info_->getNIntegrableObjects());
132
133	}
134
135	void NPT::moveA() {
136	SimInfo::MoleculeIterator i;
137	Molecule::IntegrableObjectIterator j;
138	Molecule* mol;
139	StuntDouble* integrableObject;
140	Vector3d Tb, ji;
141	RealType mass;
142	Vector3d vel;
143	Vector3d pos;
144	Vector3d frc;
145	Vector3d sc;
146	int index;
147
148	chi= currentSnapshot_->getChi();
149	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
150	loadEta();
151
152	instaTemp =thermo.getTemperature();
153	press = thermo.getPressureTensor();
154	instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
155	instaVol =thermo.getVolume();
156
157	Vector3d COM = info_->getCom();
158
159	//evolve velocity half step
160
161	calcVelScale();
162
163	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
164	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
165	integrableObject = mol->nextIntegrableObject(j)) {
166
167	vel = integrableObject->getVel();
168	frc = integrableObject->getFrc();
169
170	mass = integrableObject->getMass();
171
172	getVelScaleA(sc, vel);
173
174	// velocity half step (use chi from previous step here):
175	//vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
176	vel += dt2PhysicalConstants::energyConvert/mass frc - dt2*sc;
177	integrableObject->setVel(vel);
178
179	if (integrableObject->isDirectional()) {
180
181	// get and convert the torque to body frame
182
183	Tb = integrableObject->lab2Body(integrableObject->getTrq());
184
185	// get the angular momentum, and propagate a half step
186
187	ji = integrableObject->getJ();
188
189	//ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
190	ji += dt2PhysicalConstants::energyConvert Tb - dt2chi ji;
191
192	rotAlgo_->rotate(integrableObject, ji, dt);
193
194	integrableObject->setJ(ji);
195	}
196
197	}
198	}
199	// evolve chi and eta half step
200
201	chi += dt2 * (instaTemp / targetTemp - 1.0) / tt2;
202
203	evolveEtaA();
204
205	//calculate the integral of chidt
206	integralOfChidt += dt2 * chi;
207
208	flucQ_->moveA();
209
210
211	index = 0;
212	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
213	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
214	integrableObject = mol->nextIntegrableObject(j)) {
215	oldPos[index++] = integrableObject->getPos();
216	}
217	}
218
219	//the first estimation of r(t+dt) is equal to r(t)
220
221	for(int k = 0; k < maxIterNum_; k++) {
222	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
227	vel = integrableObject->getVel();
228	pos = integrableObject->getPos();
229
230	this->getPosScale(pos, COM, index, sc);
231
232	pos = oldPos[index] + dt * (vel + sc);
233	integrableObject->setPos(pos);
234
235	++index;
236	}
237	}
238
239	rattle_->constraintA();
240	}
241
242	// Scale the box after all the positions have been moved:
243
244	this->scaleSimBox();
245
246	currentSnapshot_->setChi(chi);
247	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
248
249	saveEta();
250	}
251
252	void NPT::moveB(void) {
253	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	RealType mass;
264
265
266	chi= currentSnapshot_->getChi();
267	integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
268	RealType oldChi = chi;
269	RealType prevChi;
270
271	loadEta();
272
273	//save velocity and angular momentum
274	index = 0;
275	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
276	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
277	integrableObject = mol->nextIntegrableObject(j)) {
278
279	oldVel[index] = integrableObject->getVel();
280
281	if (integrableObject->isDirectional())
282	oldJi[index] = integrableObject->getJ();
283
284	++index;
285	}
286	}
287
288	// do the iteration:
289	instaVol =thermo.getVolume();
290
291	for(int k = 0; k < maxIterNum_; k++) {
292	instaTemp =thermo.getTemperature();
293	instaPress =thermo.getPressure();
294
295	// 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
299	//evolve eta
300	this->evolveEtaB();
301	this->calcVelScale();
302
303	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
308	frc = integrableObject->getFrc();
309	vel = integrableObject->getVel();
310
311	mass = integrableObject->getMass();
312
313	getVelScaleB(sc, index);
314
315	// velocity half step
316	//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	integrableObject->setVel(vel);
319
320	if (integrableObject->isDirectional()) {
321	// get and convert the torque to body frame
322	Tb = integrableObject->lab2Body(integrableObject->getTrq());
323
324	//ji[j] = oldJi[3i + j] + dt2 (Tb[j] * PhysicalConstants::energyConvert - oldJi[3i+j]chi);
325	ji = oldJi[index] + dt2PhysicalConstants::energyConvertTb - dt2chioldJi[index];
326	integrableObject->setJ(ji);
327	}
328
329	++index;
330	}
331	}
332
333	rattle_->constraintB();
334
335	if ((fabs(prevChi - chi) <= chiTolerance) && this->etaConverged())
336	break;
337	}
338
339	//calculate integral of chidt
340	integralOfChidt += dt2 * chi;
341
342	currentSnapshot_->setChi(chi);
343	currentSnapshot_->setIntegralOfChiDt(integralOfChidt);
344
345	flucQ_->moveB();
346	saveEta();
347	}
348
349	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	}