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

Comparing trunk/src/integrators/NPT.cpp (file contents):
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
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.
# Line 37 | Line 28
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]  Vardeman & Gezelter, in progress (2009).                        
40   */
41  
42   #include <math.h>
# Line 46 | Line 46
46   #include "integrators/NPT.hpp"
47   #include "math/SquareMatrix3.hpp"
48   #include "primitives/Molecule.hpp"
49 < #include "utils/OOPSEConstant.hpp"
49 > #include "utils/PhysicalConstants.hpp"
50   #include "utils/simError.h"
51  
52   // Basic isotropic thermostating and barostating via the Melchionna
# Line 59 | Line 59
59   //
60   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
61  
62 < namespace oopse {
62 > namespace OpenMD {
63  
64    NPT::NPT(SimInfo* info) :
65      VelocityVerletIntegrator(info), chiTolerance(1e-6), etaTolerance(1e-6), maxIterNum_(4) {
# Line 76 | Line 76 | namespace oopse {
76        if (!simParams->haveTargetTemp()) {
77          sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
78          painCave.isFatal = 1;
79 <        painCave.severity = OOPSE_ERROR;
79 >        painCave.severity = OPENMD_ERROR;
80          simError();
81        } else {
82          targetTemp = simParams->getTargetTemp();
# Line 85 | Line 85 | namespace oopse {
85        // We must set tauThermostat
86        if (!simParams->haveTauThermostat()) {
87          sprintf(painCave.errMsg, "If you use the constant temperature\n"
88 <                "\tintegrator, you must set tauThermostat_.\n");
88 >                "\tintegrator, you must set tauThermostat.\n");
89  
90 <        painCave.severity = OOPSE_ERROR;
90 >        painCave.severity = OPENMD_ERROR;
91          painCave.isFatal = 1;
92          simError();
93        } else {
# Line 107 | Line 107 | namespace oopse {
107        if (!simParams->haveTauBarostat()) {
108          sprintf(painCave.errMsg,
109                  "If you use the NPT integrator, you must set tauBarostat.\n");
110 <        painCave.severity = OOPSE_ERROR;
110 >        painCave.severity = OPENMD_ERROR;
111          painCave.isFatal = 1;
112          simError();
113        } else {
# Line 150 | Line 150 | namespace oopse {
150      
151      instaTemp =thermo.getTemperature();
152      press = thermo.getPressureTensor();
153 <    instaPress = OOPSEConstant::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
153 >    instaPress = PhysicalConstants::pressureConvert* (press(0, 0) + press(1, 1) + press(2, 2)) / 3.0;
154      instaVol =thermo.getVolume();
155  
156      Vector3d  COM = info_->getCom();
# Line 171 | Line 171 | namespace oopse {
171          getVelScaleA(sc, vel);
172  
173          // velocity half step  (use chi from previous step here):
174 <        //vel[j] += dt2 * ((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
175 <        vel += dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
174 >        //vel[j] += dt2 * ((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
175 >        vel += dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
176          integrableObject->setVel(vel);
177  
178          if (integrableObject->isDirectional()) {
# Line 185 | Line 185 | namespace oopse {
185  
186            ji = integrableObject->getJ();
187  
188 <          //ji[j] += dt2 * (Tb[j] * OOPSEConstant::energyConvert - ji[j]*chi);
189 <          ji += dt2*OOPSEConstant::energyConvert * Tb - dt2*chi* ji;
188 >          //ji[j] += dt2 * (Tb[j] * PhysicalConstants::energyConvert - ji[j]*chi);
189 >          ji += dt2*PhysicalConstants::energyConvert * Tb - dt2*chi* ji;
190                  
191            rotAlgo->rotate(integrableObject, ji, dt);
192  
# Line 306 | Line 306 | namespace oopse {
306            getVelScaleB(sc, index);
307  
308            // velocity half step
309 <          //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * OOPSEConstant::energyConvert - sc[j]);
310 <          vel = oldVel[index] + dt2*OOPSEConstant::energyConvert/mass* frc - dt2*sc;
309 >          //vel[j] = oldVel[3 * i + j] + dt2 *((frc[j] / mass) * PhysicalConstants::energyConvert - sc[j]);
310 >          vel = oldVel[index] + dt2*PhysicalConstants::energyConvert/mass* frc - dt2*sc;
311            integrableObject->setVel(vel);
312  
313            if (integrableObject->isDirectional()) {
314              // get and convert the torque to body frame
315              Tb = integrableObject->lab2Body(integrableObject->getTrq());
316  
317 <            //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * OOPSEConstant::energyConvert - oldJi[3*i+j]*chi);
318 <            ji = oldJi[index] + dt2*OOPSEConstant::energyConvert*Tb - dt2*chi*oldJi[index];
317 >            //ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * PhysicalConstants::energyConvert - oldJi[3*i+j]*chi);
318 >            ji = oldJi[index] + dt2*PhysicalConstants::energyConvert*Tb - dt2*chi*oldJi[index];
319              integrableObject->setJ(ji);
320            }
321  

Diff Legend

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