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 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1497 by gezelter, Thu Sep 9 13:02:24 2010 UTC

# Line 1 | Line 1
1 < #include <math.h>
2 < #include "math/MatVec3.h"
3 < #include "primitives/Atom.hpp"
4 < #include "primitives/SRI.hpp"
5 < #include "primitives/AbstractClasses.hpp"
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]  Vardeman & Gezelter, in progress (2009).                        
40 > */
41 >
42   #include "brains/SimInfo.hpp"
7 #include "UseTheForce/ForceFields.hpp"
43   #include "brains/Thermo.hpp"
44 < #include "io/ReadWrite.hpp"
45 < #include "integrators/Integrator.hpp"
44 > #include "integrators/IntegratorCreator.hpp"
45 > #include "integrators/NPTxyz.hpp"
46 > #include "primitives/Molecule.hpp"
47 > #include "utils/PhysicalConstants.hpp"
48   #include "utils/simError.h"
49  
13 #ifdef IS_MPI
14 #include "brains/mpiSimulation.hpp"
15 #endif
16
50   // Basic non-isotropic thermostating and barostating via the Melchionna
51   // modification of the Hoover algorithm:
52   //
# Line 24 | Line 57
57   //
58   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
59  
60 < template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theInfo, ForceFields* the_ff):
28 <  T( theInfo, the_ff )
29 < {
30 <  GenericData* data;
31 <  DoubleArrayData * etaValue;
32 <  vector<double> etaArray;
33 <  int i,j;
60 > namespace OpenMD {
61  
62 <  for(i = 0; i < 3; i++){
63 <    for (j = 0; j < 3; j++){
62 >    
63 >  RealType NPTxyz::calcConservedQuantity(){
64  
65 <      eta[i][j] = 0.0;
66 <      oldEta[i][j] = 0.0;
67 <    }
68 <  }
65 >    chi= currentSnapshot_->getChi();
66 >    integralOfChidt = currentSnapshot_->getIntegralOfChiDt();
67 >    loadEta();
68 >    
69 >    // We need NkBT a lot, so just set it here: This is the RAW number
70 >    // of integrableObjects, so no subtraction or addition of constraints or
71 >    // orientational degrees of freedom:
72 >    NkBT = info_->getNGlobalIntegrableObjects()*PhysicalConstants::kB *targetTemp;
73  
74 +    // fkBT is used because the thermostat operates on more degrees of freedom
75 +    // than the barostat (when there are particles with orientational degrees
76 +    // of freedom).  
77 +    fkBT = info_->getNdf()*PhysicalConstants::kB *targetTemp;        
78  
79 <  if( theInfo->useInitXSstate ){
79 >    RealType conservedQuantity;
80 >    RealType totalEnergy;
81 >    RealType thermostat_kinetic;
82 >    RealType thermostat_potential;
83 >    RealType barostat_kinetic;
84 >    RealType barostat_potential;
85 >    RealType trEta;
86  
87 <    // retrieve eta array from simInfo if it exists
47 <    data = info->getProperty(ETAVALUE_ID);
48 <    if(data){
49 <      etaValue = dynamic_cast<DoubleArrayData*>(data);
50 <      
51 <      if(etaValue){
52 <        etaArray = etaValue->getData();
53 <        
54 <        for(i = 0; i < 3; i++){
55 <          for (j = 0; j < 3; j++){
56 <            eta[i][j] = etaArray[3*i+j];
57 <            oldEta[i][j] = eta[i][j];
58 <          }
59 <        }
60 <      }
61 <    }
62 <  }
63 < }
87 >    totalEnergy = thermo.getTotalE();
88  
89 < template<typename T> NPTxyz<T>::~NPTxyz() {
89 >    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * PhysicalConstants::energyConvert);
90  
91 <  // empty for now
68 < }
91 >    thermostat_potential = fkBT* integralOfChidt / PhysicalConstants::energyConvert;
92  
93 < template<typename T> void NPTxyz<T>::resetIntegrator() {
93 >    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
94 >    trEta = tmp.trace();
95  
96 <  int i, j;
96 >    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * PhysicalConstants::energyConvert);
97  
98 <  for(i = 0; i < 3; i++)
75 <    for (j = 0; j < 3; j++)
76 <      eta[i][j] = 0.0;
98 >    barostat_potential = (targetPressure * thermo.getVolume() / PhysicalConstants::pressureConvert) /PhysicalConstants::energyConvert;
99  
100 <  T::resetIntegrator();
101 < }
100 >    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
101 >      barostat_kinetic + barostat_potential;
102  
81 template<typename T> void NPTxyz<T>::evolveEtaA() {
103  
104 <  int i, j;
104 >    return conservedQuantity;
105  
85  for(i = 0; i < 3; i ++){
86    for(j = 0; j < 3; j++){
87      if( i == j)
88        eta[i][j] += dt2 *  instaVol *
89          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
90      else
91        eta[i][j] = 0.0;
92    }
106    }
107  
108 <  for(i = 0; i < 3; i++)
109 <    for (j = 0; j < 3; j++)
97 <      oldEta[i][j] = eta[i][j];
98 < }
108 >    
109 >  void NPTxyz::scaleSimBox(){
110  
111 < template<typename T> void NPTxyz<T>::evolveEtaB() {
111 >    int i,j,k;
112 >    Mat3x3d scaleMat;
113 >    RealType eta2ij, scaleFactor;
114 >    RealType bigScale, smallScale, offDiagMax;
115 >    Mat3x3d hm;
116 >    Mat3x3d hmnew;
117  
102  int i,j;
118  
104  for(i = 0; i < 3; i++)
105    for (j = 0; j < 3; j++)
106      prevEta[i][j] = eta[i][j];
119  
120 <  for(i = 0; i < 3; i ++){
109 <    for(j = 0; j < 3; j++){
110 <      if( i == j) {
111 <        eta[i][j] = oldEta[i][j] + dt2 *  instaVol *
112 <          (press[i][j] - targetPressure/p_convert) / (NkBT*tb2);
113 <      } else {
114 <        eta[i][j] = 0.0;
115 <      }
116 <    }
117 <  }
118 < }
120 >    // Scale the box after all the positions have been moved:
121  
122 < template<typename T> void NPTxyz<T>::calcVelScale(void) {
123 <  int i,j;
122 >    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
123 >    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
124  
125 <  for (i = 0; i < 3; i++ ) {
126 <    for (j = 0; j < 3; j++ ) {
127 <      vScale[i][j] = eta[i][j];
125 >    bigScale = 1.0;
126 >    smallScale = 1.0;
127 >    offDiagMax = 0.0;
128  
129 <      if (i == j) {
130 <        vScale[i][j] += chi;
129 >    for(i=0; i<3; i++){
130 >      for(j=0; j<3; j++){
131 >        scaleMat(i, j) = 0.0;
132 >        if(i==j) {
133 >          scaleMat(i, j) = 1.0;
134 >        }
135        }
136      }
131  }
132 }
137  
138 < template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
135 <  matVecMul3( vScale, vel, sc );
136 < }
138 >    for(i=0;i<3;i++){
139  
140 < template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
139 <  int j;
140 <  double myVel[3];
140 >      // calculate the scaleFactors
141  
142 <  for (j = 0; j < 3; j++)
143 <    myVel[j] = oldVel[3*index + j];
142 >      scaleFactor = exp(dt*eta(i, i));
143  
144 <  matVecMul3( vScale, myVel, sc );
146 < }
144 >      scaleMat(i, i) = scaleFactor;
145  
146 < template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
147 <                                               int index, double sc[3]){
148 <  int j;
149 <  double rj[3];
146 >      if (scaleMat(i, i) > bigScale) {
147 >        bigScale = scaleMat(i, i);
148 >      }
149 >        
150 >      if (scaleMat(i, i) < smallScale) {
151 >        smallScale = scaleMat(i, i);
152 >      }
153 >    }
154  
155 <  for(j=0; j<3; j++)
156 <    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
155 >    if ((bigScale > 1.1) || (smallScale < 0.9)) {
156 >      sprintf( painCave.errMsg,
157 >               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
158 >               " Check your tauBarostat, as it is probably too small!\n\n"
159 >               " scaleMat = [%lf\t%lf\t%lf]\n"
160 >               "            [%lf\t%lf\t%lf]\n"
161 >               "            [%lf\t%lf\t%lf]\n",
162 >               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
163 >               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
164 >               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
165 >      painCave.isFatal = 1;
166 >      simError();
167 >    } else {
168  
169 <  matVecMul3( eta, rj, sc );
170 < }
171 <
159 < template<typename T> void NPTxyz<T>::scaleSimBox( void ){
160 <
161 <  int i,j,k;
162 <  double scaleMat[3][3];
163 <  double eta2ij, scaleFactor;
164 <  double bigScale, smallScale, offDiagMax;
165 <  double hm[3][3], hmnew[3][3];
166 <
167 <
168 <
169 <  // Scale the box after all the positions have been moved:
170 <
171 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
172 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
173 <
174 <  bigScale = 1.0;
175 <  smallScale = 1.0;
176 <  offDiagMax = 0.0;
177 <
178 <  for(i=0; i<3; i++){
179 <    for(j=0; j<3; j++){
180 <      scaleMat[i][j] = 0.0;
181 <      if(i==j) scaleMat[i][j] = 1.0;
169 >      Mat3x3d hmat = currentSnapshot_->getHmat();
170 >      hmat = hmat *scaleMat;
171 >      currentSnapshot_->setHmat(hmat);
172      }
173    }
174  
175 <  for(i=0;i<3;i++){
176 <
187 <    // calculate the scaleFactors
188 <
189 <    scaleFactor = exp(dt*eta[i][i]);
190 <
191 <    scaleMat[i][i] = scaleFactor;
192 <
193 <    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
194 <    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
175 >  void NPTxyz::loadEta() {
176 >    eta= currentSnapshot_->getEta();
177    }
178  
197 //   for(i=0; i<3; i++){
198 //     for(j=0; j<3; j++){
199
200 //       // Calculate the matrix Product of the eta array (we only need
201 //       // the ij element right now):
202
203 //       eta2ij = 0.0;
204 //       for(k=0; k<3; k++){
205 //         eta2ij += eta[i][k] * eta[k][j];
206 //       }
207
208 //       scaleMat[i][j] = 0.0;
209 //       // identity matrix (see above):
210 //       if (i == j) scaleMat[i][j] = 1.0;
211 //       // Taylor expansion for the exponential truncated at second order:
212 //       scaleMat[i][j] += dt*eta[i][j]  + 0.5*dt*dt*eta2ij;
213
214 //       if (i != j)
215 //         if (fabs(scaleMat[i][j]) > offDiagMax)
216 //           offDiagMax = fabs(scaleMat[i][j]);
217 //     }
218
219 //     if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
220 //     if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
221 //   }
222
223  if ((bigScale > 1.1) || (smallScale < 0.9)) {
224    sprintf( painCave.errMsg,
225             "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
226             " Check your tauBarostat, as it is probably too small!\n\n"
227             " scaleMat = [%lf\t%lf\t%lf]\n"
228             "            [%lf\t%lf\t%lf]\n"
229             "            [%lf\t%lf\t%lf]\n",
230             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
231             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
232             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
233    painCave.isFatal = 1;
234    simError();
235  } else {
236    info->getBoxM(hm);
237    matMul3(hm, scaleMat, hmnew);
238    info->setBoxM(hmnew);
239  }
179   }
241
242 template<typename T> bool NPTxyz<T>::etaConverged() {
243  int i;
244  double diffEta, sumEta;
245
246  sumEta = 0;
247  for(i = 0; i < 3; i++)
248    sumEta += pow(prevEta[i][i] - eta[i][i], 2);
249
250  diffEta = sqrt( sumEta / 3.0 );
251
252  return ( diffEta <= etaTolerance );
253 }
254
255 template<typename T> double NPTxyz<T>::getConservedQuantity(void){
256
257  double conservedQuantity;
258  double totalEnergy;
259  double thermostat_kinetic;
260  double thermostat_potential;
261  double barostat_kinetic;
262  double barostat_potential;
263  double trEta;
264  double a[3][3], b[3][3];
265
266  totalEnergy = tStats->getTotalE();
267
268  thermostat_kinetic = fkBT * tt2 * chi * chi /
269    (2.0 * eConvert);
270
271  thermostat_potential = fkBT* integralOfChidt / eConvert;
272
273  transposeMat3(eta, a);
274  matMul3(a, eta, b);
275  trEta = matTrace3(b);
276
277  barostat_kinetic = NkBT * tb2 * trEta /
278    (2.0 * eConvert);
279
280  barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /
281    eConvert;
282
283  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
284    barostat_kinetic + barostat_potential;
285
286 //   cout.width(8);
287 //   cout.precision(8);
288
289 //   cerr << info->getTime() << "\t" << Energy << "\t" << thermostat_kinetic <<
290 //       "\t" << thermostat_potential << "\t" << barostat_kinetic <<
291 //       "\t" << barostat_potential << "\t" << conservedQuantity << endl;
292
293  return conservedQuantity;
294
295 }
296
297 template<typename T> string NPTxyz<T>::getAdditionalParameters(void){
298  string parameters;
299  const int BUFFERSIZE = 2000; // size of the read buffer
300  char buffer[BUFFERSIZE];
301
302  sprintf(buffer,"\t%G\t%G;", chi, integralOfChidt);
303  parameters += buffer;
304
305  for(int i = 0; i < 3; i++){
306    sprintf(buffer,"\t%G\t%G\t%G;", eta[i][0], eta[i][1], eta[i][2]);
307    parameters += buffer;
308  }
309
310  return parameters;
311
312 }

Comparing trunk/src/integrators/NPTxyz.cpp (property svn:keywords):
Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1497 by gezelter, Thu Sep 9 13:02:24 2010 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines