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

Comparing branches/new_design/OOPSE-2.0/src/integrators/NPTf.cpp (file contents):
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC vs.
Revision 1822 by tim, Thu Dec 2 02:08:29 2004 UTC

# Line 1 | Line 1
1 – #include <math.h>
2 –
3 – #include "math/MatVec3.h"
4 – #include "primitives/Atom.hpp"
5 – #include "primitives/SRI.hpp"
6 – #include "primitives/AbstractClasses.hpp"
1   #include "brains/SimInfo.hpp"
8 – #include "UseTheForce/ForceFields.hpp"
2   #include "brains/Thermo.hpp"
3 < #include "io/ReadWrite.hpp"
4 < #include "integrators/Integrator.hpp"
3 > #include "integrators/NPTf.hpp"
4 > #include "primitives/Molecule.hpp"
5 > #include "utils/OOPSEConstant.hpp"
6   #include "utils/simError.h"
7  
14 – #ifdef IS_MPI
15 – #include "brains/mpiSimulation.hpp"
16 – #endif
17 –
8   // Basic non-isotropic thermostating and barostating via the Melchionna
9   // modification of the Hoover algorithm:
10   //
# Line 26 | Line 16 | NPTf::NPTf (SimInfo* info): NPT(info){
16   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
17  
18   NPTf::NPTf (SimInfo* info): NPT(info){
29 –  GenericData* data;
30 –  DoubleVectorGenericData * etaValue;
31 –  int i,j;
32 –
33 –  for(i = 0; i < 3; i++){
34 –    for (j = 0; j < 3; j++){
35 –
36 –      eta(i, j) = 0.0;
37 –      oldEta(i, j) = 0.0;
38 –    }
39 –  }
40 –
41 –
42 –  if( theInfo->useInitXSstate ){
43 –    // retrieve eta array from simInfo if it exists
44 –    data = info->getPropertyByName(ETAVALUE_ID);
45 –    if(data){
46 –      etaValue = dynamic_cast<DoubleVectorGenericData*>(data);
47 –      
48 –      if(etaValue){
49 –        
50 –        for(i = 0; i < 3; i++){
51 –          for (j = 0; j < 3; j++){
52 –            eta(i, j) = (*etaValue)[3*i+j];
53 –            oldEta(i, j) = eta(i, j);
54 –          }
55 –        }
56 –      }
57 –    }
58 –  }
19  
20   }
21  
62 – NPTf::~NPTf() {
22  
64 –  // empty for now
65 – }
66 –
23   void NPTf::evolveEtaA() {
24  
25    int i, j;
# Line 71 | Line 27 | void NPTf::evolveEtaA() {
27      for(i = 0; i < 3; i ++){
28          for(j = 0; j < 3; j++){
29              if( i == j) {
30 <                eta(i, j) += dt2 *  instaVol * (press(i, j) - targetPressure/p_convert) / (NkBT*tb2);
30 >                eta(i, j) += dt2 *  instaVol * (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
31              } else {
32                  eta(i, j) += dt2 * instaVol * press(i, j) / (NkBT*tb2);
33              }
# Line 101 | Line 57 | void NPTf::evolveEtaB() {
57          for(j = 0; j < 3; j++){
58              if( i == j) {
59                  eta(i, j) = oldEta(i, j) + dt2 *  instaVol *
60 <                (press(i, j) - targetPressure/p_convert) / (NkBT*tb2);
60 >                (press(i, j) - targetPressure/OOPSEConstant::pressureConvert) / (NkBT*tb2);
61              } else {
62                  eta(i, j) = oldEta(i, j) + dt2 * instaVol * press(i, j) / (NkBT*tb2);
63              }
# Line 124 | Line 80 | void NPTf::getVelScaleA(Vector3d& sc, const Vector3d&
80    }
81   }
82  
83 < void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel);{
83 > void NPTf::getVelScaleA(Vector3d& sc, const Vector3d& vel){
84      sc = vScale * vel;
85   }
86  
# Line 132 | Line 88 | void NPTf::getPosScale(const Vector3d& pos, const Vect
88    sc = vScale * oldVel[index];
89   }
90  
91 < void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM,
136 <                           int index, Vector3d& sc) {
91 > void NPTf::getPosScale(const Vector3d& pos, const Vector3d& COM, int index, Vector3d& sc) {
92  
93      /**@todo */
94 <    Vector3d rj = (oldPos[index] + pos[j])/2.0 -COM;
94 >    Vector3d rj = (oldPos[index] + pos)/2.0 -COM;
95      sc = eta * rj;
96   }
97  
# Line 226 | Line 181 | void NPTf::scaleSimBox(){
181      painCave.isFatal = 1;
182      simError();
183    } else {
184 <    info->getBoxM(hm);
185 <    matMul3(hm, scaleMat, hmnew);
186 <    info->setBoxM(hmnew);
184 >
185 >        Mat3x3d hmat = currentSnapshot_->getHmat();
186 >        hmat = hmat *scaleMat;
187 >        currentSnapshot_->setHmat(hmat);
188 >        
189    }
190   }
191  
# Line 256 | Line 213 | double NPTf::calcConservedQuantity(){
213      double barostat_potential;
214      double trEta;
215  
216 <    totalEnergy = tStats->getTotalE();
216 >    totalEnergy = thermo.getTotalE();
217  
218 <    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * eConvert);
218 >    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
219  
220 <    thermostat_potential = fkBT* integralOfChidt / eConvert;
220 >    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
221  
222 <    trEta = (eta.transpose() * eta).trace();
222 >    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
223 >    trEta = tmp.trace();
224      
225 <    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * eConvert);
225 >    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
226  
227 <    barostat_potential = (targetPressure * tStats->getVolume() / p_convert) /eConvert;
227 >    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
228  
229      conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
230          barostat_kinetic + barostat_potential;

Diff Legend

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