# | Line 1 | Line 1 | |
---|---|---|
1 | #include <iostream> | |
2 | #include <stdlib.h> | |
3 | #include <math.h> | |
4 | < | |
4 | > | #include "Rattle.hpp" |
5 | #ifdef IS_MPI | |
6 | #include "mpiSimulation.hpp" | |
7 | #include <unistd.h> | |
# | Line 31 | Line 31 | template<typename T> Integrator<T>::Integrator(SimInfo | |
31 | } | |
32 | ||
33 | nAtoms = info->n_atoms; | |
34 | + | integrableObjects = info->integrableObjects; |
35 | ||
36 | + | rattle = new RattleFramework(info); |
37 | + | |
38 | + | if(rattle == NULL){ |
39 | + | sprintf(painCave.errMsg, |
40 | + | "Integrator::Intergrator() Error: Memory allocation error for RattleFramework" ); |
41 | + | painCave.isFatal = 1; |
42 | + | simError(); |
43 | + | } |
44 | + | |
45 | + | /* |
46 | // check for constraints | |
47 | ||
48 | constrainedA = NULL; | |
# | Line 44 | Line 55 | template<typename T> Integrator<T>::Integrator(SimInfo | |
55 | nConstrained = 0; | |
56 | ||
57 | checkConstraints(); | |
58 | + | */ |
59 | } | |
60 | ||
61 | template<typename T> Integrator<T>::~Integrator(){ | |
62 | + | if (rattle != NULL) |
63 | + | delete rattle; |
64 | + | /* |
65 | if (nConstrained){ | |
66 | delete[] constrainedA; | |
67 | delete[] constrainedB; | |
# | Line 55 | Line 70 | template<typename T> Integrator<T>::~Integrator(){ | |
70 | delete[] moved; | |
71 | delete[] oldPos; | |
72 | } | |
73 | + | */ |
74 | } | |
75 | ||
76 | + | /* |
77 | template<typename T> void Integrator<T>::checkConstraints(void){ | |
78 | isConstrained = 0; | |
79 | ||
# | Line 68 | Line 85 | template<typename T> void Integrator<T>::checkConstrai | |
85 | ||
86 | SRI** theArray; | |
87 | for (int i = 0; i < nMols; i++){ | |
88 | < | theArray = (SRI * *) molecules[i].getMyBonds(); |
88 | > | |
89 | > | theArray = (SRI * *) molecules[i].getMyBonds(); |
90 | for (int j = 0; j < molecules[i].getNBonds(); j++){ | |
91 | constrained = theArray[j]->is_constrained(); | |
92 | ||
# | Line 114 | Line 132 | template<typename T> void Integrator<T>::checkConstrai | |
132 | } | |
133 | } | |
134 | ||
135 | + | |
136 | if (nConstrained > 0){ | |
137 | isConstrained = 1; | |
138 | ||
# | Line 147 | Line 166 | template<typename T> void Integrator<T>::checkConstrai | |
166 | ||
167 | delete[] temp_con; | |
168 | } | |
169 | + | */ |
170 | ||
151 | – | |
171 | template<typename T> void Integrator<T>::integrate(void){ | |
172 | ||
173 | double runTime = info->run_time; | |
# | Line 157 | Line 176 | template<typename T> void Integrator<T>::integrate(voi | |
176 | double thermalTime = info->thermalTime; | |
177 | double resetTime = info->resetTime; | |
178 | ||
179 | < | |
179 | > | double difference; |
180 | double currSample; | |
181 | double currThermal; | |
182 | double currStatus; | |
# | Line 176 | Line 195 | template<typename T> void Integrator<T>::integrate(voi | |
195 | ||
196 | readyCheck(); | |
197 | ||
198 | + | // remove center of mass drift velocity (in case we passed in a configuration |
199 | + | // that was drifting |
200 | + | tStats->removeCOMdrift(); |
201 | + | |
202 | + | // initialize the retraints if necessary |
203 | + | if (info->useSolidThermInt && !info->useLiquidThermInt) { |
204 | + | myFF->initRestraints(); |
205 | + | } |
206 | + | |
207 | // initialize the forces before the first step | |
208 | ||
209 | calcForce(1, 1); | |
210 | ||
211 | < | //temp test |
212 | < | tStats->getPotential(); |
211 | > | //execute constraint algorithm to make sure at the very beginning the system is constrained |
212 | > | rattle->doPreConstraint(); |
213 | > | rattle->doRattleA(); |
214 | > | calcForce(1, 1); |
215 | > | rattle->doRattleB(); |
216 | ||
186 | – | if (nConstrained){ |
187 | – | preMove(); |
188 | – | constrainA(); |
189 | – | calcForce(1, 1); |
190 | – | constrainB(); |
191 | – | } |
192 | – | |
217 | if (info->setTemp){ | |
218 | thermalize(); | |
219 | } | |
# | Line 210 | Line 234 | template<typename T> void Integrator<T>::integrate(voi | |
234 | MPIcheckPoint(); | |
235 | #endif // is_mpi | |
236 | ||
237 | < | while (info->getTime() < runTime){ |
238 | < | if ((info->getTime() + dt) >= currStatus){ |
237 | > | while (info->getTime() < runTime && !stopIntegrator()){ |
238 | > | difference = info->getTime() + dt - currStatus; |
239 | > | if (difference > 0 || fabs(difference) < 1e-4 ){ |
240 | calcPot = 1; | |
241 | calcStress = 1; | |
242 | } | |
# | Line 266 | Line 291 | template<typename T> void Integrator<T>::integrate(voi | |
291 | #endif // is_mpi | |
292 | } | |
293 | ||
294 | + | // dump out a file containing the omega values for the final configuration |
295 | + | if (info->useSolidThermInt && !info->useLiquidThermInt) |
296 | + | myFF->dumpzAngle(); |
297 | + | |
298 | + | |
299 | delete dumpOut; | |
300 | delete statOut; | |
301 | } | |
# | Line 278 | Line 308 | template<typename T> void Integrator<T>::integrateStep | |
308 | startProfile(pro3); | |
309 | #endif //profile | |
310 | ||
311 | < | preMove(); |
311 | > | //save old state (position, velocity etc) |
312 | > | rattle->doPreConstraint(); |
313 | ||
314 | #ifdef PROFILE | |
315 | endProfile(pro3); | |
# | Line 300 | Line 331 | template<typename T> void Integrator<T>::integrateStep | |
331 | MPIcheckPoint(); | |
332 | #endif // is_mpi | |
333 | ||
303 | – | |
334 | // calc forces | |
305 | – | |
335 | calcForce(calcPot, calcStress); | |
336 | ||
337 | #ifdef IS_MPI | |
# | Line 332 | Line 361 | template<typename T> void Integrator<T>::moveA(void){ | |
361 | ||
362 | ||
363 | template<typename T> void Integrator<T>::moveA(void){ | |
364 | < | int i, j; |
364 | > | size_t i, j; |
365 | DirectionalAtom* dAtom; | |
366 | double Tb[3], ji[3]; | |
367 | double vel[3], pos[3], frc[3]; | |
368 | double mass; | |
369 | + | double omega; |
370 | + | |
371 | + | for (i = 0; i < integrableObjects.size() ; i++){ |
372 | + | integrableObjects[i]->getVel(vel); |
373 | + | integrableObjects[i]->getPos(pos); |
374 | + | integrableObjects[i]->getFrc(frc); |
375 | + | |
376 | + | mass = integrableObjects[i]->getMass(); |
377 | ||
341 | – | for (i = 0; i < nAtoms; i++){ |
342 | – | atoms[i]->getVel(vel); |
343 | – | atoms[i]->getPos(pos); |
344 | – | atoms[i]->getFrc(frc); |
345 | – | |
346 | – | mass = atoms[i]->getMass(); |
347 | – | |
378 | for (j = 0; j < 3; j++){ | |
379 | // velocity half step | |
380 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
# | Line 352 | Line 382 | template<typename T> void Integrator<T>::moveA(void){ | |
382 | pos[j] += dt * vel[j]; | |
383 | } | |
384 | ||
385 | < | atoms[i]->setVel(vel); |
386 | < | atoms[i]->setPos(pos); |
385 | > | integrableObjects[i]->setVel(vel); |
386 | > | integrableObjects[i]->setPos(pos); |
387 | ||
388 | < | if (atoms[i]->isDirectional()){ |
359 | < | dAtom = (DirectionalAtom *) atoms[i]; |
388 | > | if (integrableObjects[i]->isDirectional()){ |
389 | ||
390 | // get and convert the torque to body frame | |
391 | ||
392 | < | dAtom->getTrq(Tb); |
393 | < | dAtom->lab2Body(Tb); |
392 | > | integrableObjects[i]->getTrq(Tb); |
393 | > | integrableObjects[i]->lab2Body(Tb); |
394 | ||
395 | // get the angular momentum, and propagate a half step | |
396 | ||
397 | < | dAtom->getJ(ji); |
397 | > | integrableObjects[i]->getJ(ji); |
398 | ||
399 | for (j = 0; j < 3; j++) | |
400 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
401 | ||
402 | < | this->rotationPropagation( dAtom, ji ); |
402 | > | this->rotationPropagation( integrableObjects[i], ji ); |
403 | ||
404 | < | dAtom->setJ(ji); |
404 | > | integrableObjects[i]->setJ(ji); |
405 | } | |
406 | } | |
407 | ||
408 | < | if (nConstrained){ |
380 | < | constrainA(); |
381 | < | } |
408 | > | rattle->doRattleA(); |
409 | } | |
410 | ||
411 | ||
412 | template<typename T> void Integrator<T>::moveB(void){ | |
413 | int i, j; | |
387 | – | DirectionalAtom* dAtom; |
414 | double Tb[3], ji[3]; | |
415 | double vel[3], frc[3]; | |
416 | double mass; | |
417 | ||
418 | < | for (i = 0; i < nAtoms; i++){ |
419 | < | atoms[i]->getVel(vel); |
420 | < | atoms[i]->getFrc(frc); |
418 | > | for (i = 0; i < integrableObjects.size(); i++){ |
419 | > | integrableObjects[i]->getVel(vel); |
420 | > | integrableObjects[i]->getFrc(frc); |
421 | ||
422 | < | mass = atoms[i]->getMass(); |
422 | > | mass = integrableObjects[i]->getMass(); |
423 | ||
424 | // velocity half step | |
425 | for (j = 0; j < 3; j++) | |
426 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
427 | ||
428 | < | atoms[i]->setVel(vel); |
428 | > | integrableObjects[i]->setVel(vel); |
429 | ||
430 | < | if (atoms[i]->isDirectional()){ |
405 | < | dAtom = (DirectionalAtom *) atoms[i]; |
430 | > | if (integrableObjects[i]->isDirectional()){ |
431 | ||
432 | // get and convert the torque to body frame | |
433 | ||
434 | < | dAtom->getTrq(Tb); |
435 | < | dAtom->lab2Body(Tb); |
434 | > | integrableObjects[i]->getTrq(Tb); |
435 | > | integrableObjects[i]->lab2Body(Tb); |
436 | ||
437 | // get the angular momentum, and propagate a half step | |
438 | ||
439 | < | dAtom->getJ(ji); |
439 | > | integrableObjects[i]->getJ(ji); |
440 | ||
441 | for (j = 0; j < 3; j++) | |
442 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
443 | ||
444 | ||
445 | < | dAtom->setJ(ji); |
445 | > | integrableObjects[i]->setJ(ji); |
446 | } | |
447 | } | |
448 | ||
449 | < | if (nConstrained){ |
425 | < | constrainB(); |
426 | < | } |
449 | > | rattle->doRattleB(); |
450 | } | |
451 | ||
452 | + | /* |
453 | template<typename T> void Integrator<T>::preMove(void){ | |
454 | int i, j; | |
455 | double pos[3]; | |
# | Line 684 | Line 708 | template<typename T> void Integrator<T>::constrainB(vo | |
708 | simError(); | |
709 | } | |
710 | } | |
711 | < | |
711 | > | */ |
712 | template<typename T> void Integrator<T>::rotationPropagation | |
713 | < | ( DirectionalAtom* dAtom, double ji[3] ){ |
713 | > | ( StuntDouble* sd, double ji[3] ){ |
714 | ||
715 | double angle; | |
716 | double A[3][3], I[3][3]; | |
717 | + | int i, j, k; |
718 | ||
719 | // use the angular velocities to propagate the rotation matrix a | |
720 | // full time step | |
721 | ||
722 | < | dAtom->getA(A); |
723 | < | dAtom->getI(I); |
699 | < | |
700 | < | // rotate about the x-axis |
701 | < | angle = dt2 * ji[0] / I[0][0]; |
702 | < | this->rotate( 1, 2, angle, ji, A ); |
722 | > | sd->getA(A); |
723 | > | sd->getI(I); |
724 | ||
725 | < | // rotate about the y-axis |
726 | < | angle = dt2 * ji[1] / I[1][1]; |
727 | < | this->rotate( 2, 0, angle, ji, A ); |
728 | < | |
729 | < | // rotate about the z-axis |
730 | < | angle = dt * ji[2] / I[2][2]; |
731 | < | this->rotate( 0, 1, angle, ji, A); |
725 | > | if (sd->isLinear()) { |
726 | > | i = sd->linearAxis(); |
727 | > | j = (i+1)%3; |
728 | > | k = (i+2)%3; |
729 | > | |
730 | > | angle = dt2 * ji[j] / I[j][j]; |
731 | > | this->rotate( k, i, angle, ji, A ); |
732 | ||
733 | < | // rotate about the y-axis |
734 | < | angle = dt2 * ji[1] / I[1][1]; |
714 | < | this->rotate( 2, 0, angle, ji, A ); |
733 | > | angle = dt * ji[k] / I[k][k]; |
734 | > | this->rotate( i, j, angle, ji, A); |
735 | ||
736 | < | // rotate about the x-axis |
737 | < | angle = dt2 * ji[0] / I[0][0]; |
718 | < | this->rotate( 1, 2, angle, ji, A ); |
736 | > | angle = dt2 * ji[j] / I[j][j]; |
737 | > | this->rotate( k, i, angle, ji, A ); |
738 | ||
739 | < | dAtom->setA( A ); |
739 | > | } else { |
740 | > | // rotate about the x-axis |
741 | > | angle = dt2 * ji[0] / I[0][0]; |
742 | > | this->rotate( 1, 2, angle, ji, A ); |
743 | > | |
744 | > | // rotate about the y-axis |
745 | > | angle = dt2 * ji[1] / I[1][1]; |
746 | > | this->rotate( 2, 0, angle, ji, A ); |
747 | > | |
748 | > | // rotate about the z-axis |
749 | > | angle = dt * ji[2] / I[2][2]; |
750 | > | sd->addZangle(angle); |
751 | > | this->rotate( 0, 1, angle, ji, A); |
752 | > | |
753 | > | // rotate about the y-axis |
754 | > | angle = dt2 * ji[1] / I[1][1]; |
755 | > | this->rotate( 2, 0, angle, ji, A ); |
756 | > | |
757 | > | // rotate about the x-axis |
758 | > | angle = dt2 * ji[0] / I[0][0]; |
759 | > | this->rotate( 1, 2, angle, ji, A ); |
760 | > | |
761 | > | } |
762 | > | sd->setA( A ); |
763 | } | |
764 | ||
765 | template<typename T> void Integrator<T>::rotate(int axes1, int axes2, |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |