ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-2.0/src/integrators/NPTxyz.cpp
(Generate patch)

Comparing:
trunk/OOPSE-2.0/src/integrators/NPTxyz.cpp (file contents), Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
branches/new_design/OOPSE-2.0/src/integrators/NPTxyz.cpp (file contents), Revision 1837 by tim, Thu Dec 2 22:15:31 2004 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   #include "brains/SimInfo.hpp"
7 – #include "UseTheForce/ForceFields.hpp"
2   #include "brains/Thermo.hpp"
3 < #include "io/ReadWrite.hpp"
4 < #include "integrators/Integrator.hpp"
3 > #include "integrators/IntegratorCreator.hpp"
4 > #include "integrators/NPTxyz.hpp"
5 > #include "primitives/Molecule.hpp"
6 > #include "utils/OOPSEConstant.hpp"
7   #include "utils/simError.h"
8  
13 – #ifdef IS_MPI
14 – #include "brains/mpiSimulation.hpp"
15 – #endif
16 –
9   // Basic non-isotropic thermostating and barostating via the Melchionna
10   // modification of the Hoover algorithm:
11   //
# Line 24 | Line 16 | template<typename T> NPTxyz<T>::NPTxyz ( SimInfo *theI
16   //
17   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
18  
19 < 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;
19 > namespace oopse {
20  
21 <  for(i = 0; i < 3; i++){
22 <    for (j = 0; j < 3; j++){
21 > static IntegratorBuilder<NPTxyz>* NPTxyzCreator = new IntegratorBuilder<NPTxyz>("NPTxyz");
22 >    
23 > double NPTxyz::calcConservedQuantity(){
24  
25 <      eta[i][j] = 0.0;
26 <      oldEta[i][j] = 0.0;
27 <    }
28 <  }
25 >  double conservedQuantity;
26 >  double totalEnergy;
27 >  double thermostat_kinetic;
28 >  double thermostat_potential;
29 >  double barostat_kinetic;
30 >  double barostat_potential;
31 >  double trEta;
32  
33 +  totalEnergy = thermo.getTotalE();
34  
35 <  if( theInfo->useInitXSstate ){
35 >  thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
36  
37 <    // 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 < }
37 >  thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
38  
39 < template<typename T> NPTxyz<T>::~NPTxyz() {
39 >    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
40 >    trEta = tmp.trace();
41  
42 <  // empty for now
68 < }
42 >  barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
43  
44 < template<typename T> void NPTxyz<T>::resetIntegrator() {
44 >  barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
45  
46 <  int i, j;
46 >  conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
47 >    barostat_kinetic + barostat_potential;
48  
74 –  for(i = 0; i < 3; i++)
75 –    for (j = 0; j < 3; j++)
76 –      eta[i][j] = 0.0;
49  
50 <  T::resetIntegrator();
79 < }
50 >  return conservedQuantity;
51  
81 – template<typename T> void NPTxyz<T>::evolveEtaA() {
82 –
83 –  int i, j;
84 –
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 –    }
93 –  }
94 –
95 –  for(i = 0; i < 3; i++)
96 –    for (j = 0; j < 3; j++)
97 –      oldEta[i][j] = eta[i][j];
52   }
53  
54 < template<typename T> void NPTxyz<T>::evolveEtaB() {
54 >    
55 > void NPTxyz::scaleSimBox(){
56  
102 –  int i,j;
103 –
104 –  for(i = 0; i < 3; i++)
105 –    for (j = 0; j < 3; j++)
106 –      prevEta[i][j] = eta[i][j];
107 –
108 –  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 – }
119 –
120 – template<typename T> void NPTxyz<T>::calcVelScale(void) {
121 –  int i,j;
122 –
123 –  for (i = 0; i < 3; i++ ) {
124 –    for (j = 0; j < 3; j++ ) {
125 –      vScale[i][j] = eta[i][j];
126 –
127 –      if (i == j) {
128 –        vScale[i][j] += chi;
129 –      }
130 –    }
131 –  }
132 – }
133 –
134 – template<typename T> void NPTxyz<T>::getVelScaleA(double sc[3], double vel[3]) {
135 –  matVecMul3( vScale, vel, sc );
136 – }
137 –
138 – template<typename T> void NPTxyz<T>::getVelScaleB(double sc[3], int index ){
139 –  int j;
140 –  double myVel[3];
141 –
142 –  for (j = 0; j < 3; j++)
143 –    myVel[j] = oldVel[3*index + j];
144 –
145 –  matVecMul3( vScale, myVel, sc );
146 – }
147 –
148 – template<typename T> void NPTxyz<T>::getPosScale(double pos[3], double COM[3],
149 –                                               int index, double sc[3]){
150 –  int j;
151 –  double rj[3];
152 –
153 –  for(j=0; j<3; j++)
154 –    rj[j] = ( oldPos[index*3+j] + pos[j]) / 2.0 - COM[j];
155 –
156 –  matVecMul3( eta, rj, sc );
157 – }
158 –
159 – template<typename T> void NPTxyz<T>::scaleSimBox( void ){
160 –
57    int i,j,k;
58 <  double scaleMat[3][3];
58 >  Mat3x3d scaleMat;
59    double eta2ij, scaleFactor;
60    double bigScale, smallScale, offDiagMax;
61 <  double hm[3][3], hmnew[3][3];
61 >  Mat3x3d hm;
62 >  Mat3x3d hmnew;
63  
64  
65  
# Line 177 | Line 74 | template<typename T> void NPTxyz<T>::scaleSimBox( void
74  
75    for(i=0; i<3; i++){
76      for(j=0; j<3; j++){
77 <      scaleMat[i][j] = 0.0;
78 <      if(i==j) scaleMat[i][j] = 1.0;
77 >      scaleMat(i, j) = 0.0;
78 >      if(i==j) scaleMat(i, j) = 1.0;
79      }
80    }
81  
# Line 186 | Line 83 | template<typename T> void NPTxyz<T>::scaleSimBox( void
83  
84      // calculate the scaleFactors
85  
86 <    scaleFactor = exp(dt*eta[i][i]);
86 >    scaleFactor = exp(dt*eta(i, i));
87  
88 <    scaleMat[i][i] = scaleFactor;
88 >    scaleMat(i, i) = scaleFactor;
89  
90 <    if (scaleMat[i][i] > bigScale) bigScale = scaleMat[i][i];
91 <    if (scaleMat[i][i] < smallScale) smallScale = scaleMat[i][i];
90 >    if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
91 >    if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
92    }
93  
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 –
94    if ((bigScale > 1.1) || (smallScale < 0.9)) {
95      sprintf( painCave.errMsg,
96               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
# Line 227 | Line 98 | template<typename T> void NPTxyz<T>::scaleSimBox( void
98               " scaleMat = [%lf\t%lf\t%lf]\n"
99               "            [%lf\t%lf\t%lf]\n"
100               "            [%lf\t%lf\t%lf]\n",
101 <             scaleMat[0][0],scaleMat[0][1],scaleMat[0][2],
102 <             scaleMat[1][0],scaleMat[1][1],scaleMat[1][2],
103 <             scaleMat[2][0],scaleMat[2][1],scaleMat[2][2]);
101 >             scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
102 >             scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
103 >             scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
104      painCave.isFatal = 1;
105      simError();
106    } else {
107 <    info->getBoxM(hm);
108 <    matMul3(hm, scaleMat, hmnew);
109 <    info->setBoxM(hmnew);
107 >
108 >        Mat3x3d hmat = currentSnapshot_->getHmat();
109 >        hmat = hmat *scaleMat;
110 >        currentSnapshot_->setHmat(hmat);
111    }
112   }
113  
114 < template<typename T> bool NPTxyz<T>::etaConverged() {
115 <  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 );
114 > void NPTf::loadEta() {
115 >    eta= currentSnapshot_->getEta();
116   }
117  
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 –
118   }
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 – }

Diff Legend

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