81 |
|
double COM[3]; |
82 |
|
|
83 |
|
instaTemp = tStats->getTemperature(); |
84 |
+ |
|
85 |
|
tStats->getPressureTensor( press ); |
86 |
|
instaPress = p_convert * (press[0][0] + press[1][1] + press[2][2]) / 3.0; |
87 |
+ |
|
88 |
|
instaVol = tStats->getVolume(); |
89 |
|
|
90 |
|
tStats->getCOM(COM); |
179 |
|
int i, j, k; |
180 |
|
DirectionalAtom* dAtom; |
181 |
|
double Tb[3], ji[3], sc[3]; |
182 |
< |
double vel[3], frc[3]; |
182 |
> |
double vel[3], frc[3], qVel[3]; |
183 |
|
double mass; |
184 |
|
|
185 |
|
// Set things up for the iteration: |
209 |
|
|
210 |
|
for (k=0; k < 4; k++) { |
211 |
|
|
212 |
+ |
atoms[0]->getVel(vel); |
213 |
+ |
|
214 |
|
instaTemp = tStats->getTemperature(); |
215 |
+ |
|
216 |
|
instaPress = tStats->getPressure(); |
217 |
|
|
218 |
|
// evolve chi another half step using the temperature at t + dt/2 |
229 |
|
|
230 |
|
getVelScaleB( sc, i ); |
231 |
|
|
232 |
+ |
for(j=0;j<3;j++) |
233 |
+ |
qVel[j] = vel[j]; |
234 |
+ |
|
235 |
|
// velocity half step |
236 |
|
for (j=0; j < 3; j++) |
237 |
|
vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]); |
273 |
|
} |
274 |
|
|
275 |
|
void NPT::evolveChiA() { |
276 |
+ |
|
277 |
|
chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2; |
278 |
|
oldChi = chi; |
279 |
|
} |