# | 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 153 | Line 152 | template<typename T> void Integrator<T>::integrate(voi | |
152 | double sampleTime = info->sampleTime; | |
153 | double statusTime = info->statusTime; | |
154 | double thermalTime = info->thermalTime; | |
155 | + | double resetTime = info->resetTime; |
156 | ||
157 | + | |
158 | double currSample; | |
159 | double currThermal; | |
160 | double currStatus; | |
161 | < | |
161 | > | double currReset; |
162 | > | |
163 | int calcPot, calcStress; | |
164 | int isError; | |
165 | ||
# | Line 171 | Line 173 | template<typename T> void Integrator<T>::integrate(voi | |
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(); | |
191 | } | |
192 | ||
182 | – | calcPot = 0; |
183 | – | calcStress = 0; |
184 | – | currSample = sampleTime; |
185 | – | currThermal = thermalTime; |
186 | – | currStatus = statusTime; |
187 | – | |
193 | calcPot = 0; | |
194 | calcStress = 0; | |
195 | currSample = sampleTime + info->getTime(); | |
196 | currThermal = thermalTime+ info->getTime(); | |
197 | currStatus = statusTime + info->getTime(); | |
198 | + | currReset = resetTime + info->getTime(); |
199 | ||
200 | dumpOut->writeDump(info->getTime()); | |
201 | statOut->writeStat(info->getTime()); | |
202 | ||
197 | – | readyCheck(); |
203 | ||
204 | + | |
205 | #ifdef IS_MPI | |
206 | strcpy(checkPointMsg, "The integrator is ready to go."); | |
207 | MPIcheckPoint(); | |
# | Line 229 | Line 235 | template<typename T> void Integrator<T>::integrate(voi | |
235 | calcStress = 0; | |
236 | currStatus += statusTime; | |
237 | } | |
238 | + | |
239 | + | if (info->resetIntegrator){ |
240 | + | if (info->getTime() >= currReset){ |
241 | + | this->resetIntegrator(); |
242 | + | currReset += resetTime; |
243 | + | } |
244 | + | } |
245 | ||
246 | #ifdef IS_MPI | |
247 | strcpy(checkPointMsg, "successfully took a time step."); | |
# | Line 249 | Line 262 | template<typename T> void Integrator<T>::integrateStep | |
262 | ||
263 | moveA(); | |
264 | ||
252 | – | if (nConstrained){ |
253 | – | constrainA(); |
254 | – | } |
265 | ||
266 | ||
267 | + | |
268 | #ifdef IS_MPI | |
269 | strcpy(checkPointMsg, "Succesful moveA\n"); | |
270 | MPIcheckPoint(); | |
# | Line 274 | Line 285 | template<typename T> void Integrator<T>::integrateStep | |
285 | ||
286 | moveB(); | |
287 | ||
277 | – | if (nConstrained){ |
278 | – | constrainB(); |
279 | – | } |
288 | ||
289 | + | |
290 | #ifdef IS_MPI | |
291 | strcpy(checkPointMsg, "Succesful moveB\n"); | |
292 | MPIcheckPoint(); | |
# | Line 289 | Line 298 | template<typename T> void Integrator<T>::moveA(void){ | |
298 | int i, j; | |
299 | DirectionalAtom* dAtom; | |
300 | double Tb[3], ji[3]; | |
292 | – | double A[3][3], I[3][3]; |
293 | – | double angle; |
301 | double vel[3], pos[3], frc[3]; | |
302 | double mass; | |
303 | ||
# | Line 326 | Line 333 | template<typename T> void Integrator<T>::moveA(void){ | |
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 |
330 | < | // full time step |
336 | > | this->rotationPropagation( dAtom, ji ); |
337 | ||
338 | < | dAtom->getA(A); |
333 | < | dAtom->getI(I); |
334 | < | |
335 | < | // rotate about the x-axis |
336 | < | angle = dt2 * ji[0] / I[0][0]; |
337 | < | this->rotate(1, 2, angle, ji, A); |
338 | < | |
339 | < | // rotate about the y-axis |
340 | < | angle = dt2 * ji[1] / I[1][1]; |
341 | < | this->rotate(2, 0, angle, ji, A); |
342 | < | |
343 | < | // rotate about the z-axis |
344 | < | angle = dt * ji[2] / I[2][2]; |
345 | < | this->rotate(0, 1, angle, ji, A); |
346 | < | |
347 | < | // rotate about the y-axis |
348 | < | angle = dt2 * ji[1] / I[1][1]; |
349 | < | this->rotate(2, 0, angle, ji, A); |
350 | < | |
351 | < | // rotate about the x-axis |
352 | < | angle = dt2 * ji[0] / I[0][0]; |
353 | < | this->rotate(1, 2, angle, ji, A); |
354 | < | |
355 | < | |
356 | < | dAtom->setJ(ji); |
357 | < | dAtom->setA(A); |
338 | > | dAtom->setJ(ji); |
339 | } | |
340 | + | } |
341 | + | |
342 | + | if (nConstrained){ |
343 | + | constrainA(); |
344 | } | |
345 | } | |
346 | ||
# | Line 398 | Line 383 | template<typename T> void Integrator<T>::moveB(void){ | |
383 | dAtom->setJ(ji); | |
384 | } | |
385 | } | |
386 | + | |
387 | + | if (nConstrained){ |
388 | + | constrainB(); |
389 | + | } |
390 | } | |
391 | ||
392 | template<typename T> void Integrator<T>::preMove(void){ | |
# | Line 556 | Line 545 | template<typename T> void Integrator<T>::constrainA(){ | |
545 | painCave.isFatal = 1; | |
546 | simError(); | |
547 | } | |
548 | + | |
549 | } | |
550 | ||
551 | template<typename T> void Integrator<T>::constrainB(void){ | |
# | Line 657 | Line 647 | template<typename T> void Integrator<T>::constrainB(vo | |
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, | |
# | Line 749 | Line 774 | template<typename T> void Integrator<T>::thermalize(){ | |
774 | template<typename T> void Integrator<T>::thermalize(){ | |
775 | tStats->velocitize(); | |
776 | } | |
777 | + | |
778 | + | template<typename T> double Integrator<T>::getConservedQuantity(void){ |
779 | + | return tStats->getTotalE(); |
780 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |