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 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 60 | Line 60 | namespace oopse {
60   namespace oopse {
61  
62      
63 < double NPTxyz::calcConservedQuantity(){
63 >  RealType NPTxyz::calcConservedQuantity(){
64  
65      // We need NkBT a lot, so just set it here: This is the RAW number
66      // of integrableObjects, so no subtraction or addition of constraints or
# Line 72 | Line 72 | double NPTxyz::calcConservedQuantity(){
72      // of freedom).  
73      fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;        
74  
75 <    double conservedQuantity;
76 <    double totalEnergy;
77 <    double thermostat_kinetic;
78 <    double thermostat_potential;
79 <    double barostat_kinetic;
80 <    double barostat_potential;
81 <    double trEta;
75 >    RealType conservedQuantity;
76 >    RealType totalEnergy;
77 >    RealType thermostat_kinetic;
78 >    RealType thermostat_potential;
79 >    RealType barostat_kinetic;
80 >    RealType barostat_potential;
81 >    RealType trEta;
82  
83      totalEnergy = thermo.getTotalE();
84  
# Line 86 | Line 86 | double NPTxyz::calcConservedQuantity(){
86  
87      thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
88  
89 <    SquareMatrix<double, 3> tmp = eta.transpose() * eta;
89 >    SquareMatrix<RealType, 3> tmp = eta.transpose() * eta;
90      trEta = tmp.trace();
91  
92      barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
# Line 94 | Line 94 | double NPTxyz::calcConservedQuantity(){
94      barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
95  
96      conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
97 <        barostat_kinetic + barostat_potential;
97 >      barostat_kinetic + barostat_potential;
98  
99  
100      return conservedQuantity;
101  
102 < }
102 >  }
103  
104      
105 < void NPTxyz::scaleSimBox(){
105 >  void NPTxyz::scaleSimBox(){
106  
107      int i,j,k;
108      Mat3x3d scaleMat;
109 <    double eta2ij, scaleFactor;
110 <    double bigScale, smallScale, offDiagMax;
109 >    RealType eta2ij, scaleFactor;
110 >    RealType bigScale, smallScale, offDiagMax;
111      Mat3x3d hm;
112      Mat3x3d hmnew;
113  
114  
115  
116 <  // Scale the box after all the positions have been moved:
116 >    // Scale the box after all the positions have been moved:
117  
118 <  // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
119 <  //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
118 >    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
119 >    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
120  
121      bigScale = 1.0;
122      smallScale = 1.0;
123      offDiagMax = 0.0;
124  
125      for(i=0; i<3; i++){
126 <        for(j=0; j<3; j++){
127 <            scaleMat(i, j) = 0.0;
128 <            if(i==j) {
129 <                scaleMat(i, j) = 1.0;
130 <            }
131 <        }
126 >      for(j=0; j<3; j++){
127 >        scaleMat(i, j) = 0.0;
128 >        if(i==j) {
129 >          scaleMat(i, j) = 1.0;
130 >        }
131 >      }
132      }
133  
134      for(i=0;i<3;i++){
135  
136 <    // calculate the scaleFactors
136 >      // calculate the scaleFactors
137  
138 <        scaleFactor = exp(dt*eta(i, i));
138 >      scaleFactor = exp(dt*eta(i, i));
139  
140 <        scaleMat(i, i) = scaleFactor;
140 >      scaleMat(i, i) = scaleFactor;
141  
142 <        if (scaleMat(i, i) > bigScale) {
143 <            bigScale = scaleMat(i, i);
144 <        }
142 >      if (scaleMat(i, i) > bigScale) {
143 >        bigScale = scaleMat(i, i);
144 >      }
145          
146 <        if (scaleMat(i, i) < smallScale) {
147 <            smallScale = scaleMat(i, i);
148 <        }
146 >      if (scaleMat(i, i) < smallScale) {
147 >        smallScale = scaleMat(i, i);
148 >      }
149      }
150  
151      if ((bigScale > 1.1) || (smallScale < 0.9)) {
152 <        sprintf( painCave.errMsg,
153 <            "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
154 <            " Check your tauBarostat, as it is probably too small!\n\n"
155 <            " scaleMat = [%lf\t%lf\t%lf]\n"
156 <            "            [%lf\t%lf\t%lf]\n"
157 <            "            [%lf\t%lf\t%lf]\n",
158 <        scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
159 <        scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
160 <        scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
161 <        painCave.isFatal = 1;
162 <        simError();
152 >      sprintf( painCave.errMsg,
153 >               "NPTxyz error: Attempting a Box scaling of more than 10 percent.\n"
154 >               " Check your tauBarostat, as it is probably too small!\n\n"
155 >               " scaleMat = [%lf\t%lf\t%lf]\n"
156 >               "            [%lf\t%lf\t%lf]\n"
157 >               "            [%lf\t%lf\t%lf]\n",
158 >               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
159 >               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
160 >               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2));
161 >      painCave.isFatal = 1;
162 >      simError();
163      } else {
164  
165 <        Mat3x3d hmat = currentSnapshot_->getHmat();
166 <        hmat = hmat *scaleMat;
167 <        currentSnapshot_->setHmat(hmat);
165 >      Mat3x3d hmat = currentSnapshot_->getHmat();
166 >      hmat = hmat *scaleMat;
167 >      currentSnapshot_->setHmat(hmat);
168      }
169 < }
169 >  }
170  
171 < void NPTxyz::loadEta() {
171 >  void NPTxyz::loadEta() {
172      eta= currentSnapshot_->getEta();
173 < }
173 >  }
174  
175   }

Diff Legend

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