ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/ExtendedSystem.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/ExtendedSystem.cpp (file contents):
Revision 483 by gezelter, Wed Apr 9 04:06:43 2003 UTC vs.
Revision 488 by gezelter, Thu Apr 10 16:35:31 2003 UTC

# Line 13 | Line 13 | ExtendedSystem::ExtendedSystem( SimInfo* the_entry_plu
13    entry_plug = the_entry_plug;
14    zeta = 0.0;
15    epsilonDot = 0.0;
16 +  epsilonDotX = 0.0;
17 +  epsilonDotY = 0.0;
18 +  epsilonDotZ = 0.0;  
19    have_tau_thermostat = 0;
20    have_tau_barostat = 0;
21    have_target_temp = 0;
# Line 47 | Line 50 | void ExtendedSystem::NoseHooverNVT( double dt, double
50      
51      zetaScale = zeta * dt;
52      
53 <    std::cerr << "zetaScale = " << zetaScale << "\n";
53 >    //std::cerr << "zetaScale = " << zetaScale << "\n";
54      
55      // perform thermostat scaling on linear velocities and angular momentum
56      for(i = 0; i < entry_plug->n_atoms; i++){
# Line 121 | Line 124 | void ExtendedSystem::NoseHooverAndersonNPT( double dt,
124      
125      // determine the change in cell volume
126      scale = pow( (1.0 + dt * 3.0 * epsilonDot), (1.0 / 3.0));
127 <    std::cerr << "pmol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n";
127 >    //std::cerr << "pmol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n";
128      
129      newBox[0] = oldBox[0] * scale;
130      newBox[1] = oldBox[1] * scale;
# Line 141 | Line 144 | void ExtendedSystem::NoseHooverAndersonNPT( double dt,
144      zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass );
145      zetaScale = zeta * dt;
146      
147 <    std::cerr << "zetaScale = " << zetaScale << " epsilonScale = " << epsilonScale <<  "\n";
147 >    //std::cerr << "zetaScale = " << zetaScale << " epsilonScale = " << epsilonScale <<  "\n";
148      
149    // apply barostating and thermostating to velocities and angular momenta
150      for(i = 0; i < entry_plug->n_atoms; i++){
# Line 175 | Line 178 | void ExtendedSystem::NoseHooverAndersonNPT( double dt,
178    }
179   }
180  
181 +
182 + void ExtendedSystem::ConstantStress( double dt,
183 +                                     double ke,
184 +                                     double p_tensor[9] ) {
185 +
186 +  double oldBox[3];
187 +  double newBox[3];
188 +  const double kB = 8.31451e-7;     // boltzmann constant in amu*Ang^2*fs^-2/K
189 +  const double p_units = 6.10192996e-9; // converts atm to amu*fs^-2*Ang^-1
190 +  const double e_convert = 4.184e-4;    // to convert ke from kcal/mol to
191 +                                        // amu*Ang^2*fs^-2/K
192 +
193 +  int i;
194 +  double p_ext, zetaScale, epsilonScale, scale, NkBT, ke_temp;
195 +  double pX_ext, pY_ext, pZ_ext;
196 +  double volume, p_mol;
197 +  double vx, vy, vz, jx, jy, jz;
198 +  DirectionalAtom* dAtom;
199 +
200 +  if (this->NPTready()) {
201 +    atoms = entry_plug->atoms;
202 +    
203 +    p_ext = targetPressure * p_units;
204 +
205 +    pX_ext = p_ext / 3.0;
206 +    pY_ext = p_ext / 3.0;
207 +    pZ_ext = p_ext / 3.0;
208 +  
209 +    entry_plug->getBox(oldBox);
210 +    volume = oldBox[0]*oldBox[1]*oldBox[2];
211 +    
212 +    ke_temp = ke * e_convert;
213 +    NkBT = (double)entry_plug->ndf * kB * targetTemp;
214 +    
215 +    // propagate the strain rate
216 +    
217 +    epsilonDotX +=  dt * ((p_tensor[0] - pX_ext) * volume /
218 +                          (tauBarostat*tauBarostat * kB * targetTemp) );
219 +    epsilonDotY +=  dt * ((p_tensor[4] - pY_ext) * volume /
220 +                          (tauBarostat*tauBarostat * kB * targetTemp) );
221 +    epsilonDotZ +=  dt * ((p_tensor[8] - pZ_ext) * volume /
222 +                          (tauBarostat*tauBarostat * kB * targetTemp) );
223 +    
224 +    // determine the change in cell volume
225 +
226 +    //scale = pow( (1.0 + dt * 3.0 * (epsilonDot), (1.0 / 3.0));
227 +    //std::cerr << "pmol = " << p_mol << " p_ext = " << p_ext << " scale = " << scale << "\n";
228 +    
229 +    newBox[0] = oldBox[0] * scale;
230 +    newBox[1] = oldBox[1] * scale;
231 +    newBox[2] = oldBox[2] * scale;
232 +    volume = newBox[0]*newBox[1]*newBox[2];
233 +    
234 +    entry_plug->setBox(newBox);
235 +    
236 +    // perform affine transform to update positions with volume fluctuations
237 +    this->AffineTransform( oldBox, newBox );
238 +    
239 +    epsilonScale = epsilonDot * dt;
240 +    
241 +    // advance the zeta term to zeta(t + dt) - zeta is 0.0d0 on config. readin
242 +    // qmass is set in the parameter file
243 +    
244 +    zeta += dt * ( (ke_temp*2.0 - NkBT) / qmass );
245 +    zetaScale = zeta * dt;
246 +    
247 +    //std::cerr << "zetaScale = " << zetaScale << " epsilonScale = " << epsilonScale <<  "\n";
248 +    
249 +  // apply barostating and thermostating to velocities and angular momenta
250 +    for(i = 0; i < entry_plug->n_atoms; i++){
251 +      
252 +      vx = atoms[i]->get_vx();
253 +      vy = atoms[i]->get_vy();
254 +      vz = atoms[i]->get_vz();
255 +      
256 +      atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale));
257 +      atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale));
258 +      atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale));
259 +    }
260 +    if( entry_plug->n_oriented ){
261 +      
262 +      for( i=0; i < entry_plug->n_atoms; i++ ){
263 +        
264 +        if( atoms[i]->isDirectional() ){
265 +          
266 +          dAtom = (DirectionalAtom *)atoms[i];
267 +          
268 +          jx = dAtom->getJx();
269 +          jy = dAtom->getJy();
270 +          jz = dAtom->getJz();
271 +          
272 +          dAtom->setJx( jx * (1.0 - zetaScale));
273 +          dAtom->setJy( jy * (1.0 - zetaScale));
274 +          dAtom->setJz( jz * (1.0 - zetaScale));
275 +        }
276 +      }
277 +    }
278 +  }
279 + }
280 +
281   void ExtendedSystem::AffineTransform( double oldBox[3], double newBox[3] ){
282  
283    int i;

Diff Legend

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