# | Line 32 | Line 32 | template<typename T> Integrator<T>::Integrator(SimInfo | |
---|---|---|
32 | ||
33 | nAtoms = info->n_atoms; | |
34 | integrableObjects = info->integrableObjects; | |
35 | < | |
35 | > | |
36 | // check for constraints | |
37 | ||
38 | constrainedA = NULL; | |
# | Line 45 | Line 45 | template<typename T> Integrator<T>::Integrator(SimInfo | |
45 | nConstrained = 0; | |
46 | ||
47 | checkConstraints(); | |
48 | + | |
49 | } | |
50 | ||
51 | template<typename T> Integrator<T>::~Integrator(){ | |
# | Line 160 | Line 161 | template<typename T> void Integrator<T>::integrate(voi | |
161 | double thermalTime = info->thermalTime; | |
162 | double resetTime = info->resetTime; | |
163 | ||
164 | < | |
164 | > | double difference; |
165 | double currSample; | |
166 | double currThermal; | |
167 | double currStatus; | |
# | Line 179 | Line 180 | template<typename T> void Integrator<T>::integrate(voi | |
180 | ||
181 | readyCheck(); | |
182 | ||
183 | + | // remove center of mass drift velocity (in case we passed in a configuration |
184 | + | // that was drifting |
185 | + | tStats->removeCOMdrift(); |
186 | + | |
187 | + | // initialize the retraints if necessary |
188 | + | if (info->useSolidThermInt && !info->useLiquidThermInt) { |
189 | + | myFF->initRestraints(); |
190 | + | } |
191 | + | |
192 | // initialize the forces before the first step | |
193 | ||
194 | calcForce(1, 1); | |
185 | – | |
186 | – | //temp test |
187 | – | tStats->getPotential(); |
195 | ||
196 | if (nConstrained){ | |
197 | preMove(); | |
# | Line 213 | Line 220 | template<typename T> void Integrator<T>::integrate(voi | |
220 | MPIcheckPoint(); | |
221 | #endif // is_mpi | |
222 | ||
223 | < | while (info->getTime() < runTime){ |
224 | < | if ((info->getTime() + dt) >= currStatus){ |
223 | > | while (info->getTime() < runTime && !stopIntegrator()){ |
224 | > | difference = info->getTime() + dt - currStatus; |
225 | > | if (difference > 0 || fabs(difference) < 1e-4 ){ |
226 | calcPot = 1; | |
227 | calcStress = 1; | |
228 | } | |
# | Line 247 | Line 255 | template<typename T> void Integrator<T>::integrate(voi | |
255 | ||
256 | if (info->getTime() >= currStatus){ | |
257 | statOut->writeStat(info->getTime()); | |
258 | + | if (info->useSolidThermInt || info->useLiquidThermInt) |
259 | + | statOut->writeRaw(info->getTime()); |
260 | calcPot = 0; | |
261 | calcStress = 0; | |
262 | currStatus += statusTime; | |
# | Line 269 | Line 279 | template<typename T> void Integrator<T>::integrate(voi | |
279 | #endif // is_mpi | |
280 | } | |
281 | ||
282 | + | // dump out a file containing the omega values for the final configuration |
283 | + | if (info->useSolidThermInt && !info->useLiquidThermInt) |
284 | + | myFF->dumpzAngle(); |
285 | + | |
286 | + | |
287 | delete dumpOut; | |
288 | delete statOut; | |
289 | } | |
# | Line 303 | Line 318 | template<typename T> void Integrator<T>::integrateStep | |
318 | MPIcheckPoint(); | |
319 | #endif // is_mpi | |
320 | ||
306 | – | |
321 | // calc forces | |
308 | – | |
322 | calcForce(calcPot, calcStress); | |
323 | ||
324 | #ifdef IS_MPI | |
# | Line 340 | Line 353 | template<typename T> void Integrator<T>::moveA(void){ | |
353 | double Tb[3], ji[3]; | |
354 | double vel[3], pos[3], frc[3]; | |
355 | double mass; | |
356 | + | double omega; |
357 | ||
358 | for (i = 0; i < integrableObjects.size() ; i++){ | |
359 | integrableObjects[i]->getVel(vel); | |
# | Line 690 | Line 704 | template<typename T> void Integrator<T>::rotationPropa | |
704 | ||
705 | double angle; | |
706 | double A[3][3], I[3][3]; | |
707 | + | int i, j, k; |
708 | ||
709 | // use the angular velocities to propagate the rotation matrix a | |
710 | // full time step | |
# | Line 697 | Line 712 | template<typename T> void Integrator<T>::rotationPropa | |
712 | sd->getA(A); | |
713 | sd->getI(I); | |
714 | ||
715 | < | // rotate about the x-axis |
716 | < | angle = dt2 * ji[0] / I[0][0]; |
717 | < | this->rotate( 1, 2, angle, ji, A ); |
718 | < | |
719 | < | // rotate about the y-axis |
720 | < | angle = dt2 * ji[1] / I[1][1]; |
721 | < | this->rotate( 2, 0, angle, ji, A ); |
707 | < | |
708 | < | // rotate about the z-axis |
709 | < | angle = dt * ji[2] / I[2][2]; |
710 | < | this->rotate( 0, 1, angle, ji, A); |
715 | > | if (sd->isLinear()) { |
716 | > | i = sd->linearAxis(); |
717 | > | j = (i+1)%3; |
718 | > | k = (i+2)%3; |
719 | > | |
720 | > | angle = dt2 * ji[j] / I[j][j]; |
721 | > | this->rotate( k, i, angle, ji, A ); |
722 | ||
723 | < | // rotate about the y-axis |
724 | < | angle = dt2 * ji[1] / I[1][1]; |
714 | < | this->rotate( 2, 0, angle, ji, A ); |
723 | > | angle = dt * ji[k] / I[k][k]; |
724 | > | this->rotate( i, j, angle, ji, A); |
725 | ||
726 | < | // rotate about the x-axis |
727 | < | angle = dt2 * ji[0] / I[0][0]; |
718 | < | this->rotate( 1, 2, angle, ji, A ); |
726 | > | angle = dt2 * ji[j] / I[j][j]; |
727 | > | this->rotate( k, i, angle, ji, A ); |
728 | ||
729 | + | } else { |
730 | + | // rotate about the x-axis |
731 | + | angle = dt2 * ji[0] / I[0][0]; |
732 | + | this->rotate( 1, 2, angle, ji, A ); |
733 | + | |
734 | + | // rotate about the y-axis |
735 | + | angle = dt2 * ji[1] / I[1][1]; |
736 | + | this->rotate( 2, 0, angle, ji, A ); |
737 | + | |
738 | + | // rotate about the z-axis |
739 | + | angle = dt * ji[2] / I[2][2]; |
740 | + | sd->addZangle(angle); |
741 | + | this->rotate( 0, 1, angle, ji, A); |
742 | + | |
743 | + | // rotate about the y-axis |
744 | + | angle = dt2 * ji[1] / I[1][1]; |
745 | + | this->rotate( 2, 0, angle, ji, A ); |
746 | + | |
747 | + | // rotate about the x-axis |
748 | + | angle = dt2 * ji[0] / I[0][0]; |
749 | + | this->rotate( 1, 2, angle, ji, A ); |
750 | + | |
751 | + | } |
752 | sd->setA( A ); | |
753 | } | |
754 |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |