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

Comparing trunk/src/integrators/NPTxyz.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# 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 "brains/SimInfo.hpp"
# Line 44 | Line 44
44   #include "integrators/IntegratorCreator.hpp"
45   #include "integrators/NPTxyz.hpp"
46   #include "primitives/Molecule.hpp"
47 < #include "utils/OOPSEConstant.hpp"
47 > #include "utils/PhysicalConstants.hpp"
48   #include "utils/simError.h"
49  
50   // Basic non-isotropic thermostating and barostating via the Melchionna
# Line 57 | Line 57
57   //
58   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
59  
60 < namespace oopse {
60 > namespace OpenMD {
61  
62      
63 < double NPTxyz::calcConservedQuantity(){
63 >  RealType NPTxyz::calcConservedQuantity(){
64  
65      // We need NkBT a lot, so just set it here: This is the RAW number
66      // of integrableObjects, so no subtraction or addition of constraints or
67      // orientational degrees of freedom:
68 <    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;
68 >    NkBT = info_->getNGlobalIntegrableObjects()*PhysicalConstants::kB *targetTemp;
69  
70      // fkBT is used because the thermostat operates on more degrees of freedom
71      // than the barostat (when there are particles with orientational degrees
72      // of freedom).  
73 <    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;        
73 >    fkBT = info_->getNdf()*PhysicalConstants::kB *targetTemp;        
74  
75 <    double conservedQuantity;
76 <    double totalEnergy;
77 <    double thermostat_kinetic;
78 <    double thermostat_potential;
79 <    double barostat_kinetic;
80 <    double barostat_potential;
81 <    double trEta;
75 >    RealType conservedQuantity;
76 >    RealType totalEnergy;
77 >    RealType thermostat_kinetic;
78 >    RealType thermostat_potential;
79 >    RealType barostat_kinetic;
80 >    RealType barostat_potential;
81 >    RealType trEta;
82  
83      totalEnergy = thermo.getTotalE();
84  
85 <    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
85 >    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * PhysicalConstants::energyConvert);
86  
87 <    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
87 >    thermostat_potential = fkBT* integralOfChidt / PhysicalConstants::energyConvert;
88  
89 <    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
89 >    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
90      trEta = tmp.trace();
91  
92 <    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
92 >    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * PhysicalConstants::energyConvert);
93  
94 <    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
94 >    barostat_potential = (targetPressure * thermo.getVolume() / PhysicalConstants::pressureConvert) /PhysicalConstants::energyConvert;
95  
96      conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
97 <        barostat_kinetic + barostat_potential;
97 >      barostat_kinetic + barostat_potential;
98  
99  
100      return conservedQuantity;
101  
102 < }
102 >  }
103  
104      
105 < void NPTxyz::scaleSimBox(){
105 >  void NPTxyz::scaleSimBox(){
106  
107      int i,j,k;
108      Mat3x3d scaleMat;
109 <    double eta2ij, scaleFactor;
110 <    double bigScale, smallScale, offDiagMax;
109 >    RealType eta2ij, scaleFactor;
110 >    RealType bigScale, smallScale, offDiagMax;
111      Mat3x3d hm;
112      Mat3x3d hmnew;
113  
114  
115  
116 <  // Scale the box after all the positions have been moved:
116 >    // Scale the box after all the positions have been moved:
117  
118 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
119 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
118 >    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
119 >    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
120  
121      bigScale = 1.0;
122      smallScale = 1.0;
123      offDiagMax = 0.0;
124  
125      for(i=0; i<3; i++){
126 <        for(j=0; j<3; j++){
127 <            scaleMat(i, j) = 0.0;
128 <            if(i==j) {
129 <                scaleMat(i, j) = 1.0;
130 <            }
131 <        }
126 >      for(j=0; j<3; j++){
127 >        scaleMat(i, j) = 0.0;
128 >        if(i==j) {
129 >          scaleMat(i, j) = 1.0;
130 >        }
131 >      }
132      }
133  
134      for(i=0;i<3;i++){
135  
136 <    // calculate the scaleFactors
136 >      // calculate the scaleFactors
137  
138 <        scaleFactor = exp(dt*eta(i, i));
138 >      scaleFactor = exp(dt*eta(i, i));
139  
140 <        scaleMat(i, i) = scaleFactor;
140 >      scaleMat(i, i) = scaleFactor;
141  
142 <        if (scaleMat(i, i) > bigScale) {
143 <            bigScale = scaleMat(i, i);
144 <        }
142 >      if (scaleMat(i, i) > bigScale) {
143 >        bigScale = scaleMat(i, i);
144 >      }
145          
146 <        if (scaleMat(i, i) < smallScale) {
147 <            smallScale = scaleMat(i, i);
148 <        }
146 >      if (scaleMat(i, i) < smallScale) {
147 >        smallScale = scaleMat(i, i);
148 >      }
149      }
150  
151      if ((bigScale > 1.1) || (smallScale < 0.9)) {
152 <        sprintf( painCave.errMsg,
153 <            "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
154 <            " Check your tauBarostat, as it is probably too small!\n\n"
155 <            " scaleMat = [%lf\t%lf\t%lf]\n"
156 <            "            [%lf\t%lf\t%lf]\n"
157 <            "            [%lf\t%lf\t%lf]\n",
158 <        scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
159 <        scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
160 <        scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
161 <        painCave.isFatal = 1;
162 <        simError();
152 >      sprintf( painCave.errMsg,
153 >               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
154 >               " Check your tauBarostat, as it is probably too small!\n\n"
155 >               " scaleMat = [%lf\t%lf\t%lf]\n"
156 >               "            [%lf\t%lf\t%lf]\n"
157 >               "            [%lf\t%lf\t%lf]\n",
158 >               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
159 >               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
160 >               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
161 >      painCave.isFatal = 1;
162 >      simError();
163      } else {
164  
165 <        Mat3x3d hmat = currentSnapshot_->getHmat();
166 <        hmat = hmat *scaleMat;
167 <        currentSnapshot_->setHmat(hmat);
165 >      Mat3x3d hmat = currentSnapshot_->getHmat();
166 >      hmat = hmat *scaleMat;
167 >      currentSnapshot_->setHmat(hmat);
168      }
169 < }
169 >  }
170  
171 < void NPTxyz::loadEta() {
171 >  void NPTxyz::loadEta() {
172      eta= currentSnapshot_->getEta();
173 < }
173 >  }
174  
175   }

Diff Legend

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