# | Line 1 | Line 1 | |
---|---|---|
1 | #include <iostream> | |
2 | < | #include <cstdlib> |
3 | < | #include <cmath> |
2 | > | #include <stdlib.h> |
3 | > | #include <math.h> |
4 | ||
5 | #ifdef IS_MPI | |
6 | #include "mpiSimulation.hpp" | |
# | Line 25 | Line 25 | template<typename T> Integrator<T>::Integrator(SimInfo | |
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 | |
# | Line 147 | Line 146 | template<typename T> void Integrator<T>::integrate(voi | |
146 | ||
147 | ||
148 | template<typename T> void Integrator<T>::integrate(void){ | |
150 | – | int i, j; // loop counters |
149 | ||
150 | double runTime = info->run_time; | |
151 | double sampleTime = info->sampleTime; | |
# | Line 162 | Line 160 | template<typename T> void Integrator<T>::integrate(voi | |
160 | double currReset; | |
161 | ||
162 | int calcPot, calcStress; | |
165 | – | int isError; |
163 | ||
164 | tStats = new Thermo(info); | |
165 | statOut = new StatWriter(info); | |
166 | dumpOut = new DumpWriter(info); | |
167 | ||
168 | atoms = info->atoms; | |
172 | – | DirectionalAtom* dAtom; |
169 | ||
170 | dt = info->dt; | |
171 | dt2 = 0.5 * dt; | |
172 | ||
173 | + | readyCheck(); |
174 | + | |
175 | // initialize the forces before the first step | |
176 | ||
177 | calcForce(1, 1); | |
178 | + | |
179 | + | if (nConstrained){ |
180 | + | preMove(); |
181 | + | constrainA(); |
182 | + | calcForce(1, 1); |
183 | + | constrainB(); |
184 | + | } |
185 | ||
186 | if (info->setTemp){ | |
187 | thermalize(); | |
# | Line 192 | Line 197 | template<typename T> void Integrator<T>::integrate(voi | |
197 | dumpOut->writeDump(info->getTime()); | |
198 | statOut->writeStat(info->getTime()); | |
199 | ||
195 | – | readyCheck(); |
200 | ||
201 | + | |
202 | #ifdef IS_MPI | |
203 | strcpy(checkPointMsg, "The integrator is ready to go."); | |
204 | MPIcheckPoint(); | |
# | Line 241 | Line 246 | template<typename T> void Integrator<T>::integrate(voi | |
246 | #endif // is_mpi | |
247 | } | |
248 | ||
249 | < | dumpOut->writeFinal(info->getTime()); |
249 | > | |
250 | > | // write the last frame |
251 | > | dumpOut->writeDump(info->getTime()); |
252 | ||
253 | delete dumpOut; | |
254 | delete statOut; | |
# | Line 290 | Line 297 | template<typename T> void Integrator<T>::moveA(void){ | |
297 | int i, j; | |
298 | DirectionalAtom* dAtom; | |
299 | double Tb[3], ji[3]; | |
293 | – | double A[3][3], I[3][3]; |
294 | – | double angle; |
300 | double vel[3], pos[3], frc[3]; | |
301 | double mass; | |
302 | ||
# | Line 327 | Line 332 | template<typename T> void Integrator<T>::moveA(void){ | |
332 | for (j = 0; j < 3; j++) | |
333 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
334 | ||
335 | < | // 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); |
335 | > | this->rotationPropagation( dAtom, ji ); |
336 | ||
352 | – | // rotate about the x-axis |
353 | – | angle = dt2 * ji[0] / I[0][0]; |
354 | – | this->rotate(1, 2, angle, ji, A); |
355 | – | |
337 | dAtom->setJ(ji); | |
357 | – | dAtom->setA(A); |
338 | } | |
339 | } | |
340 | ||
# | Line 424 | Line 404 | template<typename T> void Integrator<T>::constrainA(){ | |
404 | } | |
405 | ||
406 | template<typename T> void Integrator<T>::constrainA(){ | |
407 | < | int i, j, k; |
407 | > | int i, j; |
408 | int done; | |
409 | double posA[3], posB[3]; | |
410 | double velA[3], velB[3]; | |
# | Line 568 | Line 548 | template<typename T> void Integrator<T>::constrainB(vo | |
548 | } | |
549 | ||
550 | template<typename T> void Integrator<T>::constrainB(void){ | |
551 | < | int i, j, k; |
551 | > | int i, j; |
552 | int done; | |
553 | double posA[3], posB[3]; | |
554 | double velA[3], velB[3]; | |
# | Line 577 | Line 557 | template<typename T> void Integrator<T>::constrainB(vo | |
557 | int a, b, ax, ay, az, bx, by, bz; | |
558 | double rma, rmb; | |
559 | double dx, dy, dz; | |
560 | < | double rabsq, pabsq, rvab; |
581 | < | double diffsq; |
560 | > | double rvab; |
561 | double gab; | |
562 | int iteration; | |
563 | ||
# | Line 668 | Line 647 | template<typename T> void Integrator<T>::constrainB(vo | |
647 | } | |
648 | } | |
649 | ||
650 | + | template<typename T> void Integrator<T>::rotationPropagation |
651 | + | ( DirectionalAtom* dAtom, double ji[3] ){ |
652 | + | |
653 | + | double angle; |
654 | + | double A[3][3], I[3][3]; |
655 | + | |
656 | + | // use the angular velocities to propagate the rotation matrix a |
657 | + | // full time step |
658 | + | |
659 | + | dAtom->getA(A); |
660 | + | dAtom->getI(I); |
661 | + | |
662 | + | // rotate about the x-axis |
663 | + | angle = dt2 * ji[0] / I[0][0]; |
664 | + | this->rotate( 1, 2, angle, ji, A ); |
665 | + | |
666 | + | // rotate about the y-axis |
667 | + | angle = dt2 * ji[1] / I[1][1]; |
668 | + | this->rotate( 2, 0, angle, ji, A ); |
669 | + | |
670 | + | // rotate about the z-axis |
671 | + | angle = dt * ji[2] / I[2][2]; |
672 | + | this->rotate( 0, 1, angle, ji, A); |
673 | + | |
674 | + | // rotate about the y-axis |
675 | + | angle = dt2 * ji[1] / I[1][1]; |
676 | + | this->rotate( 2, 0, angle, ji, A ); |
677 | + | |
678 | + | // rotate about the x-axis |
679 | + | angle = dt2 * ji[0] / I[0][0]; |
680 | + | this->rotate( 1, 2, angle, ji, A ); |
681 | + | |
682 | + | dAtom->setA( A ); |
683 | + | } |
684 | + | |
685 | template<typename T> void Integrator<T>::rotate(int axes1, int axes2, | |
686 | double angle, double ji[3], | |
687 | double A[3][3]){ |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |