| 174 |
|
double A[3][3]; // the rotation matrix |
| 175 |
|
double press[9]; |
| 176 |
|
|
| 177 |
– |
int time; |
| 178 |
– |
|
| 177 |
|
double dt = info->dt; |
| 178 |
|
double runTime = info->run_time; |
| 179 |
|
double sampleTime = info->sampleTime; |
| 180 |
|
double statusTime = info->statusTime; |
| 181 |
|
double thermalTime = info->thermalTime; |
| 182 |
|
|
| 183 |
< |
int n_loops = (int)( runTime / dt ); |
| 184 |
< |
int sample_n = (int)( sampleTime / dt ); |
| 185 |
< |
int status_n = (int)( statusTime / dt ); |
| 186 |
< |
int vel_n = (int)( thermalTime / dt ); |
| 183 |
> |
double currSample; |
| 184 |
> |
double currThermal; |
| 185 |
> |
double currStatus; |
| 186 |
> |
double currTime; |
| 187 |
|
|
| 188 |
|
int calcPot, calcStress; |
| 189 |
|
int isError; |
| 208 |
|
dump_out->writeDump( 0.0 ); |
| 209 |
|
e_out->writeStat( 0.0 ); |
| 210 |
|
|
| 211 |
< |
calcPot = 0; |
| 211 |
> |
calcPot = 0; |
| 212 |
> |
calcStress = 0; |
| 213 |
> |
currSample = sampleTime; |
| 214 |
> |
currThermal = thermalTime; |
| 215 |
> |
currStatus = statusTime; |
| 216 |
> |
currTime = 0.0;; |
| 217 |
|
|
| 218 |
< |
for( tl=0; tl<nLoops; tl++){ |
| 218 |
> |
while( currTime < runTime ){ |
| 219 |
|
|
| 220 |
+ |
if( (currTime+dt) >= currStatus ){ |
| 221 |
+ |
calcPot = 1; |
| 222 |
+ |
calcStress = 1; |
| 223 |
+ |
} |
| 224 |
+ |
|
| 225 |
|
integrateStep( calcPot, calcStress ); |
| 226 |
|
|
| 227 |
< |
time = tl + 1; |
| 228 |
< |
|
| 227 |
> |
currTime += dt; |
| 228 |
> |
|
| 229 |
|
if( info->setTemp ){ |
| 230 |
< |
if( !(time % vel_n) ) tStats->velocitize(); |
| 231 |
< |
} |
| 232 |
< |
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 233 |
< |
if( !((time+1) % status_n) ) { |
| 226 |
< |
calcPot = 1; |
| 227 |
< |
calcStress = 1; |
| 230 |
> |
if( currTime >= currThermal ){ |
| 231 |
> |
tStats->velocitize(); |
| 232 |
> |
currThermal += thermalTime; |
| 233 |
> |
} |
| 234 |
|
} |
| 235 |
< |
if( !(time % status_n) ){ |
| 235 |
> |
|
| 236 |
> |
if( currTime >= currSample ){ |
| 237 |
> |
dump_out->writeDump( currTime ); |
| 238 |
> |
currSample += sampleTime; |
| 239 |
> |
} |
| 240 |
> |
|
| 241 |
> |
if( currTime >= currStatus ){ |
| 242 |
|
e_out->writeStat( time * dt ); |
| 243 |
|
calcPot = 0; |
| 244 |
< |
if (!strcasecmp(info->ensemble, "NPT")) calcStress = 1; |
| 245 |
< |
else calcStress = 0; |
| 246 |
< |
} |
| 235 |
< |
|
| 236 |
< |
|
| 244 |
> |
calcStress = 0; |
| 245 |
> |
currStatus += statusTime; |
| 246 |
> |
} |
| 247 |
|
} |
| 248 |
|
|
| 249 |
|
dump_out->writeFinal(); |
| 250 |
|
|
| 251 |
|
delete dump_out; |
| 252 |
|
delete e_out; |
| 253 |
+ |
} |
| 254 |
+ |
|
| 255 |
+ |
void Integrator::integrateStep( int calcPot, int calcStress ){ |
| 256 |
+ |
|
| 257 |
+ |
// Position full step, and velocity half step |
| 258 |
+ |
|
| 259 |
+ |
preMove(); |
| 260 |
+ |
moveA(); |
| 261 |
+ |
if( nConstrained ) constrainA(); |
| 262 |
+ |
|
| 263 |
+ |
// calc forces |
| 264 |
+ |
|
| 265 |
+ |
myFF->doForces(calcPot,calcStress); |
| 266 |
+ |
|
| 267 |
+ |
// finish the velocity half step |
| 268 |
+ |
|
| 269 |
+ |
moveB(); |
| 270 |
+ |
if( nConstrained ) constrainB(); |
| 271 |
+ |
|
| 272 |
|
} |
| 273 |
|
|
| 274 |
|
|
