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 481 by gezelter, Tue Apr 8 21:35:49 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 84 | Line 87 | void ExtendedSystem::NoseHooverAndersonNPT( double dt,
87  
88   void ExtendedSystem::NoseHooverAndersonNPT( double dt,
89                                              double ke,
90 <                                            double p_int ) {
90 >                                            double p_tensor[9] ) {
91  
92    // Basic barostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
93    // Hoover, Phys.Rev.A, 1986, Vol.34 (3) 2499-2500
# Line 106 | Line 109 | void ExtendedSystem::NoseHooverAndersonNPT( double dt,
109      atoms = entry_plug->atoms;
110      
111      p_ext = targetPressure * p_units;
112 <    p_mol = p_int * p_units;
113 <    
112 >    p_mol = (p_tensor[0] + p_tensor[4] + p_tensor[8])/3.0;
113 >  
114      entry_plug->getBox(oldBox);
115      volume = oldBox[0]*oldBox[1]*oldBox[2];
116      
117      ke_temp = ke * e_convert;
118      NkBT = (double)entry_plug->ndf * kB * targetTemp;
119      
120 <    // propogate the strain rate
120 >    // propagate the strain rate
121      
122      epsilonDot +=  dt * ((p_mol - p_ext) * volume /
123                           (tauBarostat*tauBarostat * kB * targetTemp) );
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";
128      
129      newBox[0] = oldBox[0] * scale;
130      newBox[1] = oldBox[1] * scale;
# Line 140 | 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++){
151 >      
152 >      vx = atoms[i]->get_vx();
153 >      vy = atoms[i]->get_vy();
154 >      vz = atoms[i]->get_vz();
155 >      
156 >      atoms[i]->set_vx(vx * (1.0 - zetaScale - epsilonScale));
157 >      atoms[i]->set_vy(vy * (1.0 - zetaScale - epsilonScale));
158 >      atoms[i]->set_vz(vz * (1.0 - zetaScale - epsilonScale));
159 >    }
160 >    if( entry_plug->n_oriented ){
161 >      
162 >      for( i=0; i < entry_plug->n_atoms; i++ ){
163 >        
164 >        if( atoms[i]->isDirectional() ){
165 >          
166 >          dAtom = (DirectionalAtom *)atoms[i];
167 >          
168 >          jx = dAtom->getJx();
169 >          jy = dAtom->getJy();
170 >          jz = dAtom->getJz();
171 >          
172 >          dAtom->setJx( jx * (1.0 - zetaScale));
173 >          dAtom->setJy( jy * (1.0 - zetaScale));
174 >          dAtom->setJz( jz * (1.0 - zetaScale));
175 >        }
176 >      }
177 >    }
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        
# Line 251 | Line 355 | short int ExtendedSystem::NVTready() {
355      }
356    }
357  
358 <  return 0;
358 >  return 1;
359   }
360  
361   short int ExtendedSystem::NPTready() {
# Line 299 | Line 403 | short int ExtendedSystem::NPTready() {
403        simError();
404      }  
405    }
406 <  return 0;
406 >  return 1;
407   }
408    

Diff Legend

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