| 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; |
| 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; |
