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 539 by tim, Fri May 20 16:01:38 2005 UTC vs.
Revision 1356 by gezelter, Wed Jul 22 15:00:21 2009 UTC

# Line 57 | Line 57 | namespace oopse {
57        painCave.isFatal = 1;
58        simError();
59      } else {
60 <      surfaceTension= simParams->getSurfaceTension();
60 >      surfaceTension= simParams->getSurfaceTension()* OOPSEConstant::surfaceTensionConvert * OOPSEConstant::energyConvert;
61      }
62  
63    }
64    void NPrT::evolveEtaA() {
65      Mat3x3d hmat = currentSnapshot_->getHmat();
66 <    double hz = hmat(2, 2);
67 <    double Axy = hmat(0,0) * hmat(1, 1);
68 <    double sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
69 <    double sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::pressureConvert);
70 <    eta(0,0) -= Axy * (sx - surfaceTension) / (NkBT*tb2);
71 <    eta(1,1) -= Axy * (sy - surfaceTension) / (NkBT*tb2);
66 >    RealType hz = hmat(2, 2);
67 >    RealType Axy = hmat(0,0) * hmat(1, 1);
68 >    RealType sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
69 >    RealType sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::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/OOPSEConstant::pressureConvert) / (NkBT*tb2);
73      oldEta = eta;  
74    }
75  
76    void NPrT::evolveEtaB() {
77      Mat3x3d hmat = currentSnapshot_->getHmat();
78 <    double hz = hmat(2, 2);
79 <    double Axy = hmat(0,0) * hmat(1, 1);
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) = oldEta(0, 0) - Axy * (sx -surfaceTension) / (NkBT*tb2);
84 <    eta(1,1) = oldEta(1, 1) - Axy * (sy -surfaceTension) / (NkBT*tb2);
81 >    RealType sx = -hz * (press(0, 0) - targetPressure/OOPSEConstant::pressureConvert);
82 >    RealType sy = -hz * (press(1, 1) - targetPressure/OOPSEConstant::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);
87    }
# Line 110 | 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  
# Line 128 | Line 128 | namespace oopse {
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 140 | 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 157 | Line 157 | namespace oopse {
157      fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
158      
159  
160 <    double totalEnergy = thermo.getTotalE();
160 >    RealType totalEnergy = thermo.getTotalE();
161  
162 <    double thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
162 >    RealType thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
163  
164 <    double thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
164 >    RealType thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
165  
166 <    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
167 <    double trEta = tmp.trace();
166 >    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
167 >    RealType trEta = tmp.trace();
168      
169 <    double barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
169 >    RealType barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
170  
171 <    double barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
171 >    RealType barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
172  
173      Mat3x3d hmat = currentSnapshot_->getHmat();
174 <    double hz = hmat(2, 2);
175 <    double area = hmat(0,0) * hmat(1, 1);
174 >    RealType hz = hmat(2, 2);
175 >    RealType area = hmat(0,0) * hmat(1, 1);
176  
177 <    double conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
178 <      barostat_kinetic + barostat_potential - surfaceTension * area;
177 >    RealType conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
178 >      barostat_kinetic + barostat_potential - surfaceTension * area/ OOPSEConstant::energyConvert;
179  
180      return conservedQuantity;
181  

Diff Legend

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