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

Comparing trunk/OOPSE-2.0/src/integrators/NPrT.cpp (file contents):
Revision 2234 by tim, Thu May 19 15:49:53 2005 UTC vs.
Revision 2237 by tim, Fri May 20 19:27:03 2005 UTC

# Line 50 | Line 50 | namespace oopse {
50   namespace oopse {
51    NPrT::NPrT(SimInfo* info) : NPT(info) {
52      Globals* simParams = info_->getSimParams();
53 <    if (!simParams->haveTargetStress()) {
53 >    if (!simParams->haveSurfaceTension()) {
54        sprintf(painCave.errMsg,
55                "If you use the NPT integrator, you must set tauBarostat.\n");
56        painCave.severity = OOPSE_ERROR;
57        painCave.isFatal = 1;
58        simError();
59      } else {
60 <      targetStress= simParams->getTargetStress();
60 >      surfaceTension= simParams->getSurfaceTension();
61      }
62  
63    }
# Line 67 | Line 67 | namespace oopse {
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 - targetStress) / (NkBT*tb2);
71 <    eta(1,1) -= Axy * (sy - targetStress) / (NkBT*tb2);
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    }
# Line 80 | Line 80 | namespace oopse {
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) -= Axy * (sx -targetStress) / (NkBT*tb2);
84 <    eta(1,1) -= Axy * (sy -targetStress) / (NkBT*tb2);
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 115 | Line 115 | namespace oopse {
115    }
116  
117    void NPrT::scaleSimBox(){
118 –
119 –    int i;
120 –    int j;
121 –    int k;
118      Mat3x3d scaleMat;
123 –    double eta2ij;
124 –    double bigScale, smallScale, offDiagMax;
125 –    Mat3x3d hm;
126 –    Mat3x3d hmnew;
119  
120 +    scaleMat(0, 0) = exp(dt*eta(0, 0));
121 +    scaleMat(1, 1) = exp(dt*eta(1, 1));    
122 +    scaleMat(2, 2) = exp(dt*eta(2, 2));
123 +    Mat3x3d hmat = currentSnapshot_->getHmat();
124 +    hmat = hmat *scaleMat;
125 +    currentSnapshot_->setHmat(hmat);
126  
129 –
130 –    // Scale the box after all the positions have been moved:
131 –
132 –    // Use a taylor expansion for eta products:  Hmat = Hmat . exp(dt * etaMat)
133 –    //  Hmat = Hmat . ( Ident + dt * etaMat  + dt^2 * etaMat*etaMat / 2)
134 –
135 –    bigScale = 1.0;
136 –    smallScale = 1.0;
137 –    offDiagMax = 0.0;
138 –
139 –    for(i=0; i<3; i++){
140 –      for(j=0; j<3; j++){
141 –
142 –        // Calculate the matrix Product of the eta array (we only need
143 –        // the ij element right now):
144 –
145 –        eta2ij = 0.0;
146 –        for(k=0; k<3; k++){
147 –          eta2ij += eta(i, k) * eta(k, j);
148 –        }
149 –
150 –        scaleMat(i, j) = 0.0;
151 –        // identity matrix (see above):
152 –        if (i == j) scaleMat(i, j) = 1.0;
153 –        // Taylor expansion for the exponential truncated at second order:
154 –        scaleMat(i, j) += dt*eta(i, j)  + 0.5*dt*dt*eta2ij;
155 –      
156 –
157 –        if (i != j)
158 –          if (fabs(scaleMat(i, j)) > offDiagMax)
159 –            offDiagMax = fabs(scaleMat(i, j));
160 –      }
161 –
162 –      if (scaleMat(i, i) > bigScale) bigScale = scaleMat(i, i);
163 –      if (scaleMat(i, i) < smallScale) smallScale = scaleMat(i, i);
164 –    }
165 –
166 –    if ((bigScale > 1.01) || (smallScale < 0.99)) {
167 –      sprintf( painCave.errMsg,
168 –               "NPrT error: Attempting a Box scaling of more than 1 percent.\n"
169 –               " Check your tauBarostat, as it is probably too small!\n\n"
170 –               " scaleMat = [%lf\t%lf\t%lf]\n"
171 –               "            [%lf\t%lf\t%lf]\n"
172 –               "            [%lf\t%lf\t%lf]\n"
173 –               "      eta = [%lf\t%lf\t%lf]\n"
174 –               "            [%lf\t%lf\t%lf]\n"
175 –               "            [%lf\t%lf\t%lf]\n",
176 –               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
177 –               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
178 –               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
179 –               eta(0, 0),eta(0, 1),eta(0, 2),
180 –               eta(1, 0),eta(1, 1),eta(1, 2),
181 –               eta(2, 0),eta(2, 1),eta(2, 2));
182 –      painCave.isFatal = 1;
183 –      simError();
184 –    } else if (offDiagMax > 0.01) {
185 –      sprintf( painCave.errMsg,
186 –               "NPrT error: Attempting an off-diagonal Box scaling of more than 1 percent.\n"
187 –               " Check your tauBarostat, as it is probably too small!\n\n"
188 –               " scaleMat = [%lf\t%lf\t%lf]\n"
189 –               "            [%lf\t%lf\t%lf]\n"
190 –               "            [%lf\t%lf\t%lf]\n"
191 –               "      eta = [%lf\t%lf\t%lf]\n"
192 –               "            [%lf\t%lf\t%lf]\n"
193 –               "            [%lf\t%lf\t%lf]\n",
194 –               scaleMat(0, 0),scaleMat(0, 1),scaleMat(0, 2),
195 –               scaleMat(1, 0),scaleMat(1, 1),scaleMat(1, 2),
196 –               scaleMat(2, 0),scaleMat(2, 1),scaleMat(2, 2),
197 –               eta(0, 0),eta(0, 1),eta(0, 2),
198 –               eta(1, 0),eta(1, 1),eta(1, 2),
199 –               eta(2, 0),eta(2, 1),eta(2, 2));
200 –      painCave.isFatal = 1;
201 –      simError();
202 –    } else {
203 –
204 –      Mat3x3d hmat = currentSnapshot_->getHmat();
205 –      hmat = hmat *scaleMat;
206 –      currentSnapshot_->setHmat(hmat);
207 –        
208 –    }
127    }
128  
129    bool NPrT::etaConverged() {
# Line 238 | Line 156 | namespace oopse {
156      // of freedom).  
157      fkBT = info_->getNdf()*OOPSEConstant::kB *targetTemp;    
158      
241 –    double conservedQuantity;
242 –    double totalEnergy;
243 –    double thermostat_kinetic;
244 –    double thermostat_potential;
245 –    double barostat_kinetic;
246 –    double barostat_potential;
247 –    double trEta;
159  
160 <    totalEnergy = thermo.getTotalE();
160 >    double totalEnergy = thermo.getTotalE();
161  
162 <    thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
162 >    double thermostat_kinetic = fkBT * tt2 * chi * chi /(2.0 * OOPSEConstant::energyConvert);
163  
164 <    thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
164 >    double thermostat_potential = fkBT* integralOfChidt / OOPSEConstant::energyConvert;
165  
166      SquareMatrix<double, 3> tmp = eta.transpose() * eta;
167 <    trEta = tmp.trace();
167 >    double trEta = tmp.trace();
168      
169 <    barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
169 >    double barostat_kinetic = NkBT * tb2 * trEta /(2.0 * OOPSEConstant::energyConvert);
170  
171 <    barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
171 >    double barostat_potential = (targetPressure * thermo.getVolume() / OOPSEConstant::pressureConvert) /OOPSEConstant::energyConvert;
172  
173 <    conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
174 <      barostat_kinetic + barostat_potential;
173 >    Mat3x3d hmat = currentSnapshot_->getHmat();
174 >    double hz = hmat(2, 2);
175 >    double area = hmat(0,0) * hmat(1, 1);
176  
177 +    double conservedQuantity = totalEnergy + thermostat_kinetic + thermostat_potential +
178 +      barostat_kinetic + barostat_potential - surfaceTension * area;
179 +
180      return conservedQuantity;
181  
182    }

Diff Legend

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