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

Comparing trunk/src/integrators/NPTf.cpp (file contents):
Revision 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 963 by tim, Wed May 17 21:51:42 2006 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 49 | Line 49 | namespace oopse {
49  
50   namespace oopse {
51  
52 < // Basic non-isotropic thermostating and barostating via the Melchionna
53 < // modification of the Hoover algorithm:
54 < //
55 < //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
56 < //       Molec. Phys., 78, 533.
57 < //
58 < //           and
59 < //
60 < //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
52 >  // Basic non-isotropic thermostating and barostating via the Melchionna
53 >  // modification of the Hoover algorithm:
54 >  //
55 >  //    Melchionna, S., Ciccotti, G., and Holian, B. L., 1993,
56 >  //       Molec. Phys., 78, 533.
57 >  //
58 >  //           and
59 >  //
60 >  //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
61  
62 < void NPTf::evolveEtaA() {
62 >  void NPTf::evolveEtaA() {
63  
64 <  int i, j;
64 >    int i, j;
65  
66      for(i = 0; i < 3; i ++){
67 <        for(j = 0; j < 3; j++){
68 <            if( i == j) {
69 <                eta(i, j) += dt2 *  instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
70 <            } else {
71 <                eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2);
72 <            }
73 <        }
67 >      for(j = 0; j < 3; j++){
68 >        if( i == j) {
69 >          eta(i, j) += dt2 *  instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
70 >        } else {
71 >          eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2);
72 >        }
73 >      }
74      }
75    
76      for(i = 0; i < 3; i++) {
77 <        for (j = 0; j < 3; j++) {
77 >      for (j = 0; j < 3; j++) {
78          oldEta(i, j) = eta(i, j);
79 <        }
79 >      }
80      }
81    
82 < }
82 >  }
83  
84 < void NPTf::evolveEtaB() {
84 >  void NPTf::evolveEtaB() {
85  
86      int i;
87      int j;
88  
89      for(i = 0; i < 3; i++) {
90 <        for (j = 0; j < 3; j++) {
91 <            prevEta(i, j) = eta(i, j);
92 <        }
90 >      for (j = 0; j < 3; j++) {
91 >        prevEta(i, j) = eta(i, j);
92 >      }
93      }
94  
95      for(i = 0; i < 3; i ++){
96 <        for(j = 0; j < 3; j++){
97 <            if( i == j) {
98 <                eta(i, j) = oldEta(i, j) + dt2 *  instaVol *
99 <                (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
100 <            } else {
101 <                eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2);
102 <            }
103 <        }
96 >      for(j = 0; j < 3; j++){
97 >        if( i == j) {
98 >          eta(i, j) = oldEta(i, j) + dt2 *  instaVol *
99 >            (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
100 >        } else {
101 >          eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2);
102 >        }
103 >      }
104      }
105  
106    
107 < }
107 >  }
108  
109 < void NPTf::calcVelScale(){
109 >  void NPTf::calcVelScale(){
110  
111 <  for (int i = 0; i < 3; i++ ) {
112 <    for (int j = 0; j < 3; j++ ) {
113 <      vScale(i, j) = eta(i, j);
111 >    for (int i = 0; i < 3; i++ ) {
112 >      for (int j = 0; j < 3; j++ ) {
113 >        vScale(i, j) = eta(i, j);
114  
115 <      if (i == j) {
116 <        vScale(i, j) += chi;
115 >        if (i == j) {
116 >          vScale(i, j) += chi;
117 >        }
118        }
119      }
120    }
120 – }
121  
122 < void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
122 >  void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
123      sc = vScale * vel;
124 < }
124 >  }
125  
126 < void NPTf::getVelScaleB(Vector3d& sc, int index ) {
127 <  sc = vScale * oldVel[index];
128 < }
126 >  void NPTf::getVelScaleB(Vector3d& sc, int index ) {
127 >    sc = vScale * oldVel[index];
128 >  }
129  
130 < void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
130 >  void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
131  
132      /**@todo */
133 <    Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
133 >    Vector3d rj = (oldPos[index] + pos)/(RealType)2.0 -COM;
134      sc = eta * rj;
135 < }
135 >  }
136  
137 < void NPTf::scaleSimBox(){
137 >  void NPTf::scaleSimBox(){
138  
139 <  int i;
140 <  int j;
141 <  int k;
142 <  Mat3x3d scaleMat;
143 <  double eta2ij;
144 <  double bigScale, smallScale, offDiagMax;
145 <  Mat3x3d hm;
146 <  Mat3x3d hmnew;
139 >    int i;
140 >    int j;
141 >    int k;
142 >    Mat3x3d scaleMat;
143 >    RealType eta2ij;
144 >    RealType bigScale, smallScale, offDiagMax;
145 >    Mat3x3d hm;
146 >    Mat3x3d hmnew;
147  
148  
149  
150 <  // Scale the box after all the positions have been moved:
150 >    // Scale the box after all the positions have been moved:
151  
152 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
153 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
152 >    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
153 >    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
154  
155 <  bigScale = 1.0;
156 <  smallScale = 1.0;
157 <  offDiagMax = 0.0;
155 >    bigScale = 1.0;
156 >    smallScale = 1.0;
157 >    offDiagMax = 0.0;
158  
159 <  for(i=0; i<3; i++){
160 <    for(j=0; j<3; j++){
159 >    for(i=0; i<3; i++){
160 >      for(j=0; j<3; j++){
161  
162 <      // Calculate the matrix Product of the eta array (we only need
163 <      // the ij element right now):
162 >        // Calculate the matrix Product of the eta array (we only need
163 >        // the ij element right now):
164  
165 <      eta2ij = 0.0;
166 <      for(k=0; k<3; k++){
167 <        eta2ij += eta(i, k) * eta(k, j);
168 <      }
165 >        eta2ij = 0.0;
166 >        for(k=0; k<3; k++){
167 >          eta2ij += eta(i, k) * eta(k, j);
168 >        }
169  
170 <      scaleMat(i, j) = 0.0;
171 <      // identity matrix (see above):
172 <      if (i == j) scaleMat(i, j) = 1.0;
173 <      // Taylor expansion for the exponential truncated at second order:
174 <      scaleMat(i, j) += dt*eta(i, j)  + 0.5*dt*dt*eta2ij;
170 >        scaleMat(i, j) = 0.0;
171 >        // identity matrix (see above):
172 >        if (i == j) scaleMat(i, j) = 1.0;
173 >        // Taylor expansion for the exponential truncated at second order:
174 >        scaleMat(i, j) += dt*eta(i, j)  + 0.5*dt*dt*eta2ij;
175        
176  
177 <      if (i != j)
178 <        if (fabs(scaleMat(i, j)) > offDiagMax)
179 <          offDiagMax = fabs(scaleMat(i, j));
177 >        if (i != j)
178 >          if (fabs(scaleMat(i, j)) > offDiagMax)
179 >            offDiagMax = fabs(scaleMat(i, j));
180 >      }
181 >
182 >      if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
183 >      if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
184      }
185  
186 <    if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
187 <    if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
188 <  }
186 >    if ((bigScale > 1.01) || (smallScale < 0.99)) {
187 >      sprintf( painCave.errMsg,
188 >               "NPTf error: Attempting a Box scaling of more than 1 percent.\n"
189 >               " Check your tauBarostat, as it is probably too small!\n\n"
190 >               " scaleMat = [%lf\t%lf\t%lf]\n"
191 >               "            [%lf\t%lf\t%lf]\n"
192 >               "            [%lf\t%lf\t%lf]\n"
193 >               "      eta = [%lf\t%lf\t%lf]\n"
194 >               "            [%lf\t%lf\t%lf]\n"
195 >               "            [%lf\t%lf\t%lf]\n",
196 >               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
197 >               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
198 >               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
199 >               eta(0, 0),eta(0, 1),eta(0, 2),
200 >               eta(1, 0),eta(1, 1),eta(1, 2),
201 >               eta(2, 0),eta(2, 1),eta(2, 2));
202 >      painCave.isFatal = 1;
203 >      simError();
204 >    } else if (offDiagMax > 0.01) {
205 >      sprintf( painCave.errMsg,
206 >               "NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
207 >               " Check your tauBarostat, as it is probably too small!\n\n"
208 >               " scaleMat = [%lf\t%lf\t%lf]\n"
209 >               "            [%lf\t%lf\t%lf]\n"
210 >               "            [%lf\t%lf\t%lf]\n"
211 >               "      eta = [%lf\t%lf\t%lf]\n"
212 >               "            [%lf\t%lf\t%lf]\n"
213 >               "            [%lf\t%lf\t%lf]\n",
214 >               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
215 >               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
216 >               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
217 >               eta(0, 0),eta(0, 1),eta(0, 2),
218 >               eta(1, 0),eta(1, 1),eta(1, 2),
219 >               eta(2, 0),eta(2, 1),eta(2, 2));
220 >      painCave.isFatal = 1;
221 >      simError();
222 >    } else {
223  
224 <  if ((bigScale > 1.01) || (smallScale < 0.99)) {
225 <    sprintf( painCave.errMsg,
226 <             "NPTf error: Attempting a Box scaling of more than 1 percent.\n"
189 <             " Check your tauBarostat, as it is probably too small!\n\n"
190 <             " scaleMat = [%lf\t%lf\t%lf]\n"
191 <             "            [%lf\t%lf\t%lf]\n"
192 <             "            [%lf\t%lf\t%lf]\n"
193 <             "      eta = [%lf\t%lf\t%lf]\n"
194 <             "            [%lf\t%lf\t%lf]\n"
195 <             "            [%lf\t%lf\t%lf]\n",
196 <             scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
197 <             scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
198 <             scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
199 <             eta(0, 0),eta(0, 1),eta(0, 2),
200 <             eta(1, 0),eta(1, 1),eta(1, 2),
201 <             eta(2, 0),eta(2, 1),eta(2, 2));
202 <    painCave.isFatal = 1;
203 <    simError();
204 <  } else if (offDiagMax > 0.01) {
205 <    sprintf( painCave.errMsg,
206 <             "NPTf error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
207 <             " Check your tauBarostat, as it is probably too small!\n\n"
208 <             " scaleMat = [%lf\t%lf\t%lf]\n"
209 <             "            [%lf\t%lf\t%lf]\n"
210 <             "            [%lf\t%lf\t%lf]\n"
211 <             "      eta = [%lf\t%lf\t%lf]\n"
212 <             "            [%lf\t%lf\t%lf]\n"
213 <             "            [%lf\t%lf\t%lf]\n",
214 <             scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
215 <             scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
216 <             scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
217 <             eta(0, 0),eta(0, 1),eta(0, 2),
218 <             eta(1, 0),eta(1, 1),eta(1, 2),
219 <             eta(2, 0),eta(2, 1),eta(2, 2));
220 <    painCave.isFatal = 1;
221 <    simError();
222 <  } else {
223 <
224 <        Mat3x3d hmat = currentSnapshot_->getHmat();
225 <        hmat = hmat *scaleMat;
226 <        currentSnapshot_->setHmat(hmat);
224 >      Mat3x3d hmat = currentSnapshot_->getHmat();
225 >      hmat = hmat *scaleMat;
226 >      currentSnapshot_->setHmat(hmat);
227          
228 +    }
229    }
229 – }
230  
231 < bool NPTf::etaConverged() {
231 >  bool NPTf::etaConverged() {
232      int i;
233 <    double diffEta, sumEta;
233 >    RealType diffEta, sumEta;
234  
235      sumEta = 0;
236      for(i = 0; i < 3; i++) {
237 <        sumEta += pow(prevEta(i, i) - eta(i, i), 2);
237 >      sumEta += pow(prevEta(i, i) - eta(i, i), 2);
238      }
239      
240      diffEta = sqrt( sumEta / 3.0 );
241  
242      return ( diffEta <= etaTolerance );
243 < }
243 >  }
244  
245 < double NPTf::calcConservedQuantity(){
245 >  RealType NPTf::calcConservedQuantity(){
246  
247      chi= currentSnapshot_->getChi();
248      integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
# Line 258 | Line 258 | double NPTf::calcConservedQuantity(){
258      // of freedom).  
259      fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
260      
261 <    double conservedQuantity;
262 <    double totalEnergy;
263 <    double thermostat_kinetic;
264 <    double thermostat_potential;
265 <    double barostat_kinetic;
266 <    double barostat_potential;
267 <    double trEta;
261 >    RealType conservedQuantity;
262 >    RealType totalEnergy;
263 >    RealType thermostat_kinetic;
264 >    RealType thermostat_potential;
265 >    RealType barostat_kinetic;
266 >    RealType barostat_potential;
267 >    RealType trEta;
268  
269      totalEnergy = thermo.getTotalE();
270  
# Line 272 | Line 272 | double NPTf::calcConservedQuantity(){
272  
273      thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
274  
275 <    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
275 >    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
276      trEta = tmp.trace();
277      
278      barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
# Line 280 | Line 280 | double NPTf::calcConservedQuantity(){
280      barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
281  
282      conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
283 <        barostat_kinetic + barostat_potential;
283 >      barostat_kinetic + barostat_potential;
284  
285      return conservedQuantity;
286  
287 < }
287 >  }
288  
289 < void NPTf::loadEta() {
289 >  void NPTf::loadEta() {
290      eta= currentSnapshot_->getEta();
291  
292      //if (!eta.isDiagonal()) {
# Line 295 | Line 295 | void NPTf::loadEta() {
295      //    painCave.isFatal = 1;
296      //    simError();
297      //}
298 < }
298 >  }
299  
300 < void NPTf::saveEta() {
300 >  void NPTf::saveEta() {
301      currentSnapshot_->setEta(eta);
302 < }
302 >  }
303  
304   }

Diff Legend

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