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

Comparing trunk/src/integrators/NPrT.cpp (file contents):
Revision 536 by tim, Thu May 19 04:28:26 2005 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 "brains/SimInfo.hpp"
# Line 44 | Line 44
44   #include "integrators/IntegratorCreator.hpp"
45   #include "integrators/NPrT.hpp"
46   #include "primitives/Molecule.hpp"
47 < #include "utils/OOPSEConstant.hpp"
47 > #include "utils/PhysicalConstants.hpp"
48   #include "utils/simError.h"
49  
50 < namespace oopse {
50 > namespace OpenMD {
51    NPrT::NPrT(SimInfo* info) : NPT(info) {
52      Globals* simParams = info_->getSimParams();
53 <    if (!simParams->haveTargetStress())
53 >    if (!simParams->haveSurfaceTension()) {
54        sprintf(painCave.errMsg,
55                "If you use the NPT integrator, you must set tauBarostat.\n");
56 <      painCave.severity = OOPSE_ERROR;
56 >      painCave.severity = OPENMD_ERROR;
57        painCave.isFatal = 1;
58        simError();
59      } else {
60 <      targetStress= simParams->getTargetStress();
60 >      surfaceTension= simParams->getSurfaceTension()* PhysicalConstants::surfaceTensionConvert * PhysicalConstants::energyConvert;
61      }
62  
63    }
64    void NPrT::evolveEtaA() {
65 <    double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
66 <    double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert);
67 <    eta(0,0) -= Axy * (sx - targetStress) / (NkBT*tb2);
68 <    eta(1,1) -= Axy * (sy - targetStress) / (NkBT*tb2);
69 <    eta(2,2) += dt2 *  instaVol * (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
65 >    Mat3x3d hmat = currentSnapshot_->getHmat();
66 >    RealType hz = hmat(2, 2);
67 >    RealType Axy = hmat(0,0) * hmat(1, 1);
68 >    RealType sx = -hz * (press(0, 0) - targetPressure/PhysicalConstants::pressureConvert);
69 >    RealType sy = -hz * (press(1, 1) - targetPressure/PhysicalConstants::pressureConvert);
70 >    eta(0,0) -= dt2* Axy * (sx - surfaceTension) / (NkBT*tb2);
71 >    eta(1,1) -= dt2* Axy * (sy - surfaceTension) / (NkBT*tb2);
72 >    eta(2,2) += dt2 *  instaVol * (press(2, 2) - targetPressure/PhysicalConstants::pressureConvert) / (NkBT*tb2);
73      oldEta = eta;  
74    }
75  
76    void NPrT::evolveEtaB() {
77 <
77 >    Mat3x3d hmat = currentSnapshot_->getHmat();
78 >    RealType hz = hmat(2, 2);
79 >    RealType Axy = hmat(0,0) * hmat(1, 1);
80      prevEta = eta;
81 <    double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
82 <    double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert);
83 <    eta(0,0) -= Axy * (sx -targetStress) / (NkBT*tb2);
84 <    eta(1,1) -= Axy * (sy -targetStress) / (NkBT*tb2);
81 >    RealType sx = -hz * (press(0, 0) - targetPressure/PhysicalConstants::pressureConvert);
82 >    RealType sy = -hz * (press(1, 1) - targetPressure/PhysicalConstants::pressureConvert);
83 >    eta(0,0) = oldEta(0, 0) - dt2 * Axy * (sx -surfaceTension) / (NkBT*tb2);
84 >    eta(1,1) = oldEta(1, 1) - dt2 * Axy * (sy -surfaceTension) / (NkBT*tb2);
85      eta(2,2) = oldEta(2, 2) + dt2 *  instaVol *
86 <            (press(2, 2) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
86 >            (press(2, 2) - targetPressure/PhysicalConstants::pressureConvert) / (NkBT*tb2);
87    }
88  
89    void NPrT::calcVelScale(){
# Line 105 | Line 110 | namespace oopse {
110    void NPrT::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
111  
112      /**@todo */
113 <    Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
113 >    Vector3d rj = (oldPos[index] + pos)/(RealType)2.0 -COM;
114      sc = eta * rj;
115    }
116  
117    void NPrT::scaleSimBox(){
113 –
114 –    int i;
115 –    int j;
116 –    int k;
118      Mat3x3d scaleMat;
118 –    double eta2ij;
119 –    double bigScale, smallScale, offDiagMax;
120 –    Mat3x3d hm;
121 –    Mat3x3d hmnew;
119  
120 +    scaleMat(0, 0) = exp(dt*eta(0, 0));
121 +    scaleMat(1, 1) = exp(dt*eta(1, 1));    
122 +    scaleMat(2, 2) = exp(dt*eta(2, 2));
123 +    Mat3x3d hmat = currentSnapshot_->getHmat();
124 +    hmat = hmat *scaleMat;
125 +    currentSnapshot_->setHmat(hmat);
126  
124 –
125 –    // Scale the box after all the positions have been moved:
126 –
127 –    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
128 –    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
129 –
130 –    bigScale = 1.0;
131 –    smallScale = 1.0;
132 –    offDiagMax = 0.0;
133 –
134 –    for(i=0; i<3; i++){
135 –      for(j=0; j<3; j++){
136 –
137 –        // Calculate the matrix Product of the eta array (we only need
138 –        // the ij element right now):
139 –
140 –        eta2ij = 0.0;
141 –        for(k=0; k<3; k++){
142 –          eta2ij += eta(i, k) * eta(k, j);
143 –        }
144 –
145 –        scaleMat(i, j) = 0.0;
146 –        // identity matrix (see above):
147 –        if (i == j) scaleMat(i, j) = 1.0;
148 –        // Taylor expansion for the exponential truncated at second order:
149 –        scaleMat(i, j) += dt*eta(i, j)  + 0.5*dt*dt*eta2ij;
150 –      
151 –
152 –        if (i != j)
153 –          if (fabs(scaleMat(i, j)) > offDiagMax)
154 –            offDiagMax = fabs(scaleMat(i, j));
155 –      }
156 –
157 –      if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
158 –      if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
159 –    }
160 –
161 –    if ((bigScale > 1.01) || (smallScale < 0.99)) {
162 –      sprintf( painCave.errMsg,
163 –               "NPrT error: Attempting a Box scaling of more than 1 percent.\n"
164 –               " Check your tauBarostat, as it is probably too small!\n\n"
165 –               " scaleMat = [%lf\t%lf\t%lf]\n"
166 –               "            [%lf\t%lf\t%lf]\n"
167 –               "            [%lf\t%lf\t%lf]\n"
168 –               "      eta = [%lf\t%lf\t%lf]\n"
169 –               "            [%lf\t%lf\t%lf]\n"
170 –               "            [%lf\t%lf\t%lf]\n",
171 –               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
172 –               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
173 –               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
174 –               eta(0, 0),eta(0, 1),eta(0, 2),
175 –               eta(1, 0),eta(1, 1),eta(1, 2),
176 –               eta(2, 0),eta(2, 1),eta(2, 2));
177 –      painCave.isFatal = 1;
178 –      simError();
179 –    } else if (offDiagMax > 0.01) {
180 –      sprintf( painCave.errMsg,
181 –               "NPrT error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
182 –               " Check your tauBarostat, as it is probably too small!\n\n"
183 –               " scaleMat = [%lf\t%lf\t%lf]\n"
184 –               "            [%lf\t%lf\t%lf]\n"
185 –               "            [%lf\t%lf\t%lf]\n"
186 –               "      eta = [%lf\t%lf\t%lf]\n"
187 –               "            [%lf\t%lf\t%lf]\n"
188 –               "            [%lf\t%lf\t%lf]\n",
189 –               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
190 –               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
191 –               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
192 –               eta(0, 0),eta(0, 1),eta(0, 2),
193 –               eta(1, 0),eta(1, 1),eta(1, 2),
194 –               eta(2, 0),eta(2, 1),eta(2, 2));
195 –      painCave.isFatal = 1;
196 –      simError();
197 –    } else {
198 –
199 –      Mat3x3d hmat = currentSnapshot_->getHmat();
200 –      hmat = hmat *scaleMat;
201 –      currentSnapshot_->setHmat(hmat);
202 –        
203 –    }
127    }
128  
129    bool NPrT::etaConverged() {
130      int i;
131 <    double diffEta, sumEta;
131 >    RealType diffEta, sumEta;
132  
133      sumEta = 0;
134      for(i = 0; i < 3; i++) {
# Line 217 | Line 140 | namespace oopse {
140      return ( diffEta <= etaTolerance );
141    }
142  
143 <  double NPrT::calcConservedQuantity(){
143 >  RealType NPrT::calcConservedQuantity(){
144  
145      chi= currentSnapshot_->getChi();
146      integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
# Line 226 | Line 149 | namespace oopse {
149      // We need NkBT a lot, so just set it here: This is the RAW number
150      // of integrableObjects, so no subtraction or addition of constraints or
151      // orientational degrees of freedom:
152 <    NkBT = info_->getNGlobalIntegrableObjects()*OOPSEConstant::kB *targetTemp;
152 >    NkBT = info_->getNGlobalIntegrableObjects()*PhysicalConstants::kB *targetTemp;
153  
154      // fkBT is used because the thermostat operates on more degrees of freedom
155      // than the barostat (when there are particles with orientational degrees
156      // of freedom).  
157 <    fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
157 >    fkBT = info_->getNdf()*PhysicalConstants::kB *targetTemp;    
158      
236 –    double conservedQuantity;
237 –    double totalEnergy;
238 –    double thermostat_kinetic;
239 –    double thermostat_potential;
240 –    double barostat_kinetic;
241 –    double barostat_potential;
242 –    double trEta;
159  
160 <    totalEnergy = thermo.getTotalE();
160 >    RealType totalEnergy = thermo.getTotalE();
161  
162 <    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
162 >    RealType thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * PhysicalConstants::energyConvert);
163  
164 <    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
164 >    RealType thermostat_potential = fkBT* integralOfChidt / PhysicalConstants::energyConvert;
165  
166 <    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
167 <    trEta = tmp.trace();
166 >    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
167 >    RealType trEta = tmp.trace();
168      
169 <    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
169 >    RealType barostat_kinetic = NkBT * tb2 * trEta /(2.0 * PhysicalConstants::energyConvert);
170  
171 <    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
171 >    RealType barostat_potential = (targetPressure * thermo.getVolume() / PhysicalConstants::pressureConvert) /PhysicalConstants::energyConvert;
172  
173 <    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
174 <      barostat_kinetic + barostat_potential;
173 >    Mat3x3d hmat = currentSnapshot_->getHmat();
174 >    RealType hz = hmat(2, 2);
175 >    RealType area = hmat(0,0) * hmat(1, 1);
176  
177 +    RealType conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
178 +      barostat_kinetic + barostat_potential - surfaceTension * area/ PhysicalConstants::energyConvert;
179 +
180      return conservedQuantity;
181  
182    }

Diff Legend

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