| 25 |
|
if (info->the_integrator != NULL){ |
| 26 |
|
delete info->the_integrator; |
| 27 |
|
} |
| 28 |
< |
info->the_integrator = this; |
| 29 |
< |
|
| 28 |
> |
|
| 29 |
|
nAtoms = info->n_atoms; |
| 30 |
|
|
| 31 |
|
// check for constraints |
| 173 |
|
dt = info->dt; |
| 174 |
|
dt2 = 0.5 * dt; |
| 175 |
|
|
| 176 |
+ |
readyCheck(); |
| 177 |
+ |
|
| 178 |
|
// initialize the forces before the first step |
| 179 |
|
|
| 180 |
|
calcForce(1, 1); |
| 181 |
+ |
|
| 182 |
+ |
if (nConstrained){ |
| 183 |
+ |
preMove(); |
| 184 |
+ |
constrainA(); |
| 185 |
+ |
calcForce(1, 1); |
| 186 |
+ |
constrainB(); |
| 187 |
+ |
} |
| 188 |
|
|
| 189 |
|
if (info->setTemp){ |
| 190 |
|
thermalize(); |
| 200 |
|
dumpOut->writeDump(info->getTime()); |
| 201 |
|
statOut->writeStat(info->getTime()); |
| 202 |
|
|
| 195 |
– |
readyCheck(); |
| 203 |
|
|
| 204 |
+ |
|
| 205 |
|
#ifdef IS_MPI |
| 206 |
|
strcpy(checkPointMsg, "The integrator is ready to go."); |
| 207 |
|
MPIcheckPoint(); |
| 298 |
|
int i, j; |
| 299 |
|
DirectionalAtom* dAtom; |
| 300 |
|
double Tb[3], ji[3]; |
| 293 |
– |
double A[3][3], I[3][3]; |
| 294 |
– |
double angle; |
| 301 |
|
double vel[3], pos[3], frc[3]; |
| 302 |
|
double mass; |
| 303 |
|
|
| 333 |
|
for (j = 0; j < 3; j++) |
| 334 |
|
ji[j] += (dt2 * Tb[j]) * eConvert; |
| 335 |
|
|
| 336 |
< |
// use the angular velocities to propagate the rotation matrix a |
| 331 |
< |
// full time step |
| 332 |
< |
|
| 333 |
< |
dAtom->getA(A); |
| 334 |
< |
dAtom->getI(I); |
| 335 |
< |
|
| 336 |
< |
// rotate about the x-axis |
| 337 |
< |
angle = dt2 * ji[0] / I[0][0]; |
| 338 |
< |
this->rotate(1, 2, angle, ji, A); |
| 339 |
< |
|
| 340 |
< |
// rotate about the y-axis |
| 341 |
< |
angle = dt2 * ji[1] / I[1][1]; |
| 342 |
< |
this->rotate(2, 0, angle, ji, A); |
| 343 |
< |
|
| 344 |
< |
// rotate about the z-axis |
| 345 |
< |
angle = dt * ji[2] / I[2][2]; |
| 346 |
< |
this->rotate(0, 1, angle, ji, A); |
| 347 |
< |
|
| 348 |
< |
// rotate about the y-axis |
| 349 |
< |
angle = dt2 * ji[1] / I[1][1]; |
| 350 |
< |
this->rotate(2, 0, angle, ji, A); |
| 351 |
< |
|
| 352 |
< |
// rotate about the x-axis |
| 353 |
< |
angle = dt2 * ji[0] / I[0][0]; |
| 354 |
< |
this->rotate(1, 2, angle, ji, A); |
| 336 |
> |
this->rotationPropagation( dAtom, ji ); |
| 337 |
|
|
| 338 |
|
dAtom->setJ(ji); |
| 357 |
– |
dAtom->setA(A); |
| 339 |
|
} |
| 340 |
|
} |
| 341 |
|
|
| 647 |
|
painCave.isFatal = 1; |
| 648 |
|
simError(); |
| 649 |
|
} |
| 650 |
+ |
} |
| 651 |
+ |
|
| 652 |
+ |
template<typename T> void Integrator<T>::rotationPropagation |
| 653 |
+ |
( DirectionalAtom* dAtom, double ji[3] ){ |
| 654 |
+ |
|
| 655 |
+ |
double angle; |
| 656 |
+ |
double A[3][3], I[3][3]; |
| 657 |
+ |
|
| 658 |
+ |
// use the angular velocities to propagate the rotation matrix a |
| 659 |
+ |
// full time step |
| 660 |
+ |
|
| 661 |
+ |
dAtom->getA(A); |
| 662 |
+ |
dAtom->getI(I); |
| 663 |
+ |
|
| 664 |
+ |
// rotate about the x-axis |
| 665 |
+ |
angle = dt2 * ji[0] / I[0][0]; |
| 666 |
+ |
this->rotate( 1, 2, angle, ji, A ); |
| 667 |
+ |
|
| 668 |
+ |
// rotate about the y-axis |
| 669 |
+ |
angle = dt2 * ji[1] / I[1][1]; |
| 670 |
+ |
this->rotate( 2, 0, angle, ji, A ); |
| 671 |
+ |
|
| 672 |
+ |
// rotate about the z-axis |
| 673 |
+ |
angle = dt * ji[2] / I[2][2]; |
| 674 |
+ |
this->rotate( 0, 1, angle, ji, A); |
| 675 |
+ |
|
| 676 |
+ |
// rotate about the y-axis |
| 677 |
+ |
angle = dt2 * ji[1] / I[1][1]; |
| 678 |
+ |
this->rotate( 2, 0, angle, ji, A ); |
| 679 |
+ |
|
| 680 |
+ |
// rotate about the x-axis |
| 681 |
+ |
angle = dt2 * ji[0] / I[0][0]; |
| 682 |
+ |
this->rotate( 1, 2, angle, ji, A ); |
| 683 |
+ |
|
| 684 |
+ |
dAtom->setA( A ); |
| 685 |
|
} |
| 686 |
|
|
| 687 |
|
template<typename T> void Integrator<T>::rotate(int axes1, int axes2, |
