| 9 |
|
#include "simError.h" |
| 10 |
|
|
| 11 |
|
|
| 12 |
< |
// Basic isotropic thermostating and barostating via the Melchionna |
| 12 |
> |
// Basic non-isotropic thermostating and barostating via the Melchionna |
| 13 |
|
// modification of the Hoover algorithm: |
| 14 |
|
// |
| 15 |
|
// Melchionna, S., Ciccotti, G., and Holian, B. L., 1993, |
| 39 |
|
double Tb[3]; |
| 40 |
|
double ji[3]; |
| 41 |
|
double rj[3]; |
| 42 |
+ |
double ident[3][3], eta1[3][3], eta2[3][3], hmnew[3][3]; |
| 43 |
+ |
double hm[9]; |
| 44 |
+ |
double vx, vy, vz; |
| 45 |
+ |
double scx, scy, scz; |
| 46 |
|
double instaTemp, instaPress, instaVol; |
| 47 |
|
double tt2, tb2; |
| 48 |
|
double angle; |
| 164 |
|
} |
| 165 |
|
|
| 166 |
|
// Scale the box after all the positions have been moved: |
| 163 |
– |
|
| 167 |
|
|
| 168 |
+ |
// Use a taylor expansion for eta products: Hmat = Hmat . exp(dt * etaMat) |
| 169 |
+ |
// Hmat = Hmat . ( Ident + dt * etaMat + dt^2 * etaMat*etaMat / 2) |
| 170 |
|
|
| 171 |
< |
// Use a taylor expansion for eta products |
| 171 |
> |
|
| 172 |
> |
for(i=0; i<3; i++){ |
| 173 |
> |
for(j=0; j<3; j++){ |
| 174 |
> |
ident[i][j] = 0.0; |
| 175 |
> |
eta1[i][j] = eta[3*i+j]; |
| 176 |
> |
eta2[i][j] = 0.0; |
| 177 |
> |
for(k=0; k<3; k++){ |
| 178 |
> |
eta2[i][j] += eta[3*i+k] * eta[3*k+j]; |
| 179 |
> |
} |
| 180 |
> |
} |
| 181 |
> |
ident[i][i] = 1.0; |
| 182 |
> |
} |
| 183 |
> |
|
| 184 |
|
|
| 185 |
|
info->getBoxM(hm); |
| 186 |
+ |
|
| 187 |
+ |
for(i=0; i<3; i++){ |
| 188 |
+ |
for(j=0; j<3; j++){ |
| 189 |
+ |
hmnew[i][j] = 0.0; |
| 190 |
+ |
for(k=0; k<3; k++){ |
| 191 |
+ |
// remember that hmat has transpose ordering for Fortran compat: |
| 192 |
+ |
hmnew[i][j] += hm[3*k+i] * (ident[k][j] |
| 193 |
+ |
+ dt * eta1[k][j] |
| 194 |
+ |
+ 0.5 * dt * dt * eta2[k][j]); |
| 195 |
+ |
} |
| 196 |
+ |
} |
| 197 |
+ |
} |
| 198 |
|
|
| 199 |
+ |
for (i = 0; i < 3; i++) { |
| 200 |
+ |
for (j = 0; j < 3; j++) { |
| 201 |
+ |
// remember that hmat has transpose ordering for Fortran compat: |
| 202 |
+ |
hm[3*j + 1] = hmnew[i][j]; |
| 203 |
+ |
} |
| 204 |
+ |
} |
| 205 |
|
|
| 206 |
< |
|
| 207 |
< |
|
| 173 |
< |
|
| 174 |
< |
|
| 175 |
< |
info->scaleBox(exp(dt*eta)); |
| 176 |
< |
|
| 177 |
< |
|
| 206 |
> |
info->setBoxM(hm); |
| 207 |
> |
|
| 208 |
|
} |
| 209 |
|
|
| 210 |
< |
void NPTi::moveB( void ){ |
| 210 |
> |
void NPTf::moveB( void ){ |
| 211 |
|
int i,j,k; |
| 212 |
|
int atomIndex; |
| 213 |
|
DirectionalAtom* dAtom; |
| 214 |
|
double Tb[3]; |
| 215 |
|
double ji[3]; |
| 216 |
< |
double instaTemp, instaPress, instaVol; |
| 216 |
> |
double press[9]; |
| 217 |
> |
double instaTemp, instaVol; |
| 218 |
|
double tt2, tb2; |
| 219 |
+ |
double vx, vy, vz; |
| 220 |
+ |
double scx, scy, scz; |
| 221 |
+ |
const double p_convert = 1.63882576e8; |
| 222 |
|
|
| 223 |
|
tt2 = tauThermostat * tauThermostat; |
| 224 |
|
tb2 = tauBarostat * tauBarostat; |
| 225 |
|
|
| 226 |
|
instaTemp = tStats->getTemperature(); |
| 227 |
< |
instaPress = tStats->getPressure(); |
| 194 |
< |
instaVol = tStats->getVolume(); |
| 227 |
> |
tStats->getPressureTensor(press); |
| 228 |
|
|
| 229 |
+ |
for (i=0; i < 9; i++) press[i] *= p_convert; |
| 230 |
+ |
|
| 231 |
+ |
instaVol = tStats->getVolume(); |
| 232 |
+ |
|
| 233 |
+ |
// first evolve chi a half step |
| 234 |
+ |
|
| 235 |
|
chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
| 197 |
– |
eta += dt2 * ( instaVol * (instaPress - targetPressure) / (NkBT*tb2)); |
| 236 |
|
|
| 237 |
+ |
eta[0] += dt2 * instaVol * (press[0] - targetPressure) / (NkBT*tb2); |
| 238 |
+ |
eta[1] += dt2 * instaVol * press[1] / (NkBT*tb2); |
| 239 |
+ |
eta[2] += dt2 * instaVol * press[2] / (NkBT*tb2); |
| 240 |
+ |
eta[3] += dt2 * instaVol * press[3] / (NkBT*tb2); |
| 241 |
+ |
eta[4] += dt2 * instaVol * (press[4] - targetPressure) / (NkBT*tb2); |
| 242 |
+ |
eta[5] += dt2 * instaVol * press[5] / (NkBT*tb2); |
| 243 |
+ |
eta[6] += dt2 * instaVol * press[6] / (NkBT*tb2); |
| 244 |
+ |
eta[7] += dt2 * instaVol * press[7] / (NkBT*tb2); |
| 245 |
+ |
eta[8] += dt2 * instaVol * (press[8] - targetPressure) / (NkBT*tb2); |
| 246 |
+ |
|
| 247 |
|
for( i=0; i<nAtoms; i++ ){ |
| 248 |
|
atomIndex = i * 3; |
| 249 |
< |
|
| 249 |
> |
|
| 250 |
|
// velocity half step |
| 203 |
– |
for( j=atomIndex; j<(atomIndex+3); j++ ) |
| 204 |
– |
for( j=atomIndex; j<(atomIndex+3); j++ ) |
| 205 |
– |
vel[j] += dt2 * ((frc[j]/atoms[i]->getMass())*eConvert |
| 206 |
– |
- vel[j]*(chi+eta)); |
| 251 |
|
|
| 252 |
+ |
vx = vel[atomIndex]; |
| 253 |
+ |
vy = vel[atomIndex+1]; |
| 254 |
+ |
vz = vel[atomIndex+2]; |
| 255 |
+ |
|
| 256 |
+ |
scx = (chi + eta[0])*vx + eta[1]*vy + eta[2]*vz; |
| 257 |
+ |
scy = eta[3]*vx + (chi + eta[4])*vy + eta[5]*vz; |
| 258 |
+ |
scz = eta[6]*vx + eta[7]*vy + (chi + eta[8])*vz; |
| 259 |
+ |
|
| 260 |
+ |
vx += dt2 * ((frc[atomIndex] /atoms[i]->getMass())*eConvert - scx); |
| 261 |
+ |
vy += dt2 * ((frc[atomIndex+1]/atoms[i]->getMass())*eConvert - scy); |
| 262 |
+ |
vz += dt2 * ((frc[atomIndex+2]/atoms[i]->getMass())*eConvert - scz); |
| 263 |
+ |
|
| 264 |
+ |
vel[atomIndex] = vx; |
| 265 |
+ |
vel[atomIndex+1] = vy; |
| 266 |
+ |
vel[atomIndex+2] = vz; |
| 267 |
+ |
|
| 268 |
|
if( atoms[i]->isDirectional() ){ |
| 269 |
|
|
| 270 |
|
dAtom = (DirectionalAtom *)atoms[i]; |
