# | 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" | |
7 | #include <unistd.h> | |
8 | #endif //is_mpi | |
9 | ||
10 | + | #ifdef PROFILE |
11 | + | #include "mdProfile.hpp" |
12 | + | #endif // profile |
13 | + | |
14 | #include "Integrator.hpp" | |
15 | #include "simError.h" | |
16 | ||
# | Line 25 | Line 29 | template<typename T> Integrator<T>::Integrator(SimInfo | |
29 | if (info->the_integrator != NULL){ | |
30 | delete info->the_integrator; | |
31 | } | |
28 | – | info->the_integrator = this; |
32 | ||
33 | nAtoms = info->n_atoms; | |
34 | + | integrableObjects = info->integrableObjects; |
35 | ||
36 | // check for constraints | |
37 | ||
# | Line 65 | Line 69 | template<typename T> void Integrator<T>::checkConstrai | |
69 | ||
70 | SRI** theArray; | |
71 | for (int i = 0; i < nMols; i++){ | |
72 | < | theArray = (SRI * *) molecules[i].getMyBonds(); |
72 | > | |
73 | > | theArray = (SRI * *) molecules[i].getMyBonds(); |
74 | for (int j = 0; j < molecules[i].getNBonds(); j++){ | |
75 | constrained = theArray[j]->is_constrained(); | |
76 | ||
# | Line 111 | Line 116 | template<typename T> void Integrator<T>::checkConstrai | |
116 | } | |
117 | } | |
118 | ||
119 | + | |
120 | if (nConstrained > 0){ | |
121 | isConstrained = 1; | |
122 | ||
# | Line 132 | Line 138 | template<typename T> void Integrator<T>::checkConstrai | |
138 | } | |
139 | ||
140 | ||
141 | < | // save oldAtoms to check for lode balanceing later on. |
141 | > | // save oldAtoms to check for lode balancing later on. |
142 | ||
143 | oldAtoms = nAtoms; | |
144 | ||
# | Line 147 | Line 153 | template<typename T> void Integrator<T>::integrate(voi | |
153 | ||
154 | ||
155 | template<typename T> void Integrator<T>::integrate(void){ | |
150 | – | int i, j; // loop counters |
156 | ||
157 | double runTime = info->run_time; | |
158 | double sampleTime = info->sampleTime; | |
# | Line 160 | Line 165 | template<typename T> void Integrator<T>::integrate(voi | |
165 | double currThermal; | |
166 | double currStatus; | |
167 | double currReset; | |
168 | < | |
168 | > | |
169 | int calcPot, calcStress; | |
165 | – | int isError; |
170 | ||
171 | tStats = new Thermo(info); | |
172 | statOut = new StatWriter(info); | |
173 | dumpOut = new DumpWriter(info); | |
174 | ||
175 | atoms = info->atoms; | |
172 | – | DirectionalAtom* dAtom; |
176 | ||
177 | dt = info->dt; | |
178 | dt2 = 0.5 * dt; | |
179 | ||
180 | + | readyCheck(); |
181 | + | |
182 | // initialize the forces before the first step | |
183 | ||
184 | calcForce(1, 1); | |
185 | < | |
185 | > | |
186 | if (nConstrained){ | |
187 | preMove(); | |
188 | constrainA(); | |
189 | < | calcForce(1, 1); |
189 | > | calcForce(1, 1); |
190 | constrainB(); | |
191 | } | |
192 | ||
# | Line 199 | Line 204 | template<typename T> void Integrator<T>::integrate(voi | |
204 | dumpOut->writeDump(info->getTime()); | |
205 | statOut->writeStat(info->getTime()); | |
206 | ||
202 | – | readyCheck(); |
207 | ||
208 | #ifdef IS_MPI | |
209 | strcpy(checkPointMsg, "The integrator is ready to go."); | |
210 | MPIcheckPoint(); | |
211 | #endif // is_mpi | |
212 | ||
213 | < | while (info->getTime() < runTime){ |
213 | > | while (info->getTime() < runTime && !stopIntegrator()){ |
214 | if ((info->getTime() + dt) >= currStatus){ | |
215 | calcPot = 1; | |
216 | calcStress = 1; | |
217 | } | |
218 | ||
219 | + | #ifdef PROFILE |
220 | + | startProfile( pro1 ); |
221 | + | #endif |
222 | + | |
223 | integrateStep(calcPot, calcStress); | |
224 | ||
225 | + | #ifdef PROFILE |
226 | + | endProfile( pro1 ); |
227 | + | |
228 | + | startProfile( pro2 ); |
229 | + | #endif // profile |
230 | + | |
231 | info->incrTime(dt); | |
232 | ||
233 | if (info->setTemp){ | |
# | Line 229 | Line 243 | template<typename T> void Integrator<T>::integrate(voi | |
243 | } | |
244 | ||
245 | if (info->getTime() >= currStatus){ | |
246 | < | statOut->writeStat(info->getTime()); |
247 | < | calcPot = 0; |
246 | > | statOut->writeStat(info->getTime()); |
247 | > | calcPot = 0; |
248 | calcStress = 0; | |
249 | currStatus += statusTime; | |
250 | < | } |
250 | > | } |
251 | ||
252 | if (info->resetIntegrator){ | |
253 | if (info->getTime() >= currReset){ | |
# | Line 241 | Line 255 | template<typename T> void Integrator<T>::integrate(voi | |
255 | currReset += resetTime; | |
256 | } | |
257 | } | |
258 | + | |
259 | + | #ifdef PROFILE |
260 | + | endProfile( pro2 ); |
261 | + | #endif //profile |
262 | ||
263 | #ifdef IS_MPI | |
264 | strcpy(checkPointMsg, "successfully took a time step."); | |
# | Line 248 | Line 266 | template<typename T> void Integrator<T>::integrate(voi | |
266 | #endif // is_mpi | |
267 | } | |
268 | ||
251 | – | dumpOut->writeFinal(info->getTime()); |
252 | – | |
269 | delete dumpOut; | |
270 | delete statOut; | |
271 | } | |
# | Line 257 | Line 273 | template<typename T> void Integrator<T>::integrateStep | |
273 | template<typename T> void Integrator<T>::integrateStep(int calcPot, | |
274 | int calcStress){ | |
275 | // Position full step, and velocity half step | |
276 | + | |
277 | + | #ifdef PROFILE |
278 | + | startProfile(pro3); |
279 | + | #endif //profile |
280 | + | |
281 | preMove(); | |
282 | ||
283 | < | moveA(); |
283 | > | #ifdef PROFILE |
284 | > | endProfile(pro3); |
285 | ||
286 | + | startProfile(pro4); |
287 | + | #endif // profile |
288 | ||
289 | + | moveA(); |
290 | ||
291 | + | #ifdef PROFILE |
292 | + | endProfile(pro4); |
293 | + | |
294 | + | startProfile(pro5); |
295 | + | #endif//profile |
296 | ||
297 | + | |
298 | #ifdef IS_MPI | |
299 | strcpy(checkPointMsg, "Succesful moveA\n"); | |
300 | MPIcheckPoint(); | |
# | Line 279 | Line 310 | template<typename T> void Integrator<T>::integrateStep | |
310 | MPIcheckPoint(); | |
311 | #endif // is_mpi | |
312 | ||
313 | + | #ifdef PROFILE |
314 | + | endProfile( pro5 ); |
315 | ||
316 | + | startProfile( pro6 ); |
317 | + | #endif //profile |
318 | + | |
319 | // finish the velocity half step | |
320 | ||
321 | moveB(); | |
322 | ||
323 | + | #ifdef PROFILE |
324 | + | endProfile(pro6); |
325 | + | #endif // profile |
326 | ||
288 | – | |
327 | #ifdef IS_MPI | |
328 | strcpy(checkPointMsg, "Succesful moveB\n"); | |
329 | MPIcheckPoint(); | |
# | Line 294 | Line 332 | template<typename T> void Integrator<T>::moveA(void){ | |
332 | ||
333 | ||
334 | template<typename T> void Integrator<T>::moveA(void){ | |
335 | < | int i, j; |
335 | > | size_t i, j; |
336 | DirectionalAtom* dAtom; | |
337 | double Tb[3], ji[3]; | |
338 | double vel[3], pos[3], frc[3]; | |
339 | double mass; | |
340 | < | |
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(); |
340 | > | |
341 | > | for (i = 0; i < integrableObjects.size() ; i++){ |
342 | > | integrableObjects[i]->getVel(vel); |
343 | > | integrableObjects[i]->getPos(pos); |
344 | > | integrableObjects[i]->getFrc(frc); |
345 | > | |
346 | > | mass = integrableObjects[i]->getMass(); |
347 | ||
348 | for (j = 0; j < 3; j++){ | |
349 | // velocity half step | |
# | Line 314 | Line 352 | template<typename T> void Integrator<T>::moveA(void){ | |
352 | pos[j] += dt * vel[j]; | |
353 | } | |
354 | ||
355 | < | atoms[i]->setVel(vel); |
356 | < | atoms[i]->setPos(pos); |
355 | > | integrableObjects[i]->setVel(vel); |
356 | > | integrableObjects[i]->setPos(pos); |
357 | ||
358 | < | if (atoms[i]->isDirectional()){ |
321 | < | dAtom = (DirectionalAtom *) atoms[i]; |
358 | > | if (integrableObjects[i]->isDirectional()){ |
359 | ||
360 | // get and convert the torque to body frame | |
361 | ||
362 | < | dAtom->getTrq(Tb); |
363 | < | dAtom->lab2Body(Tb); |
362 | > | integrableObjects[i]->getTrq(Tb); |
363 | > | integrableObjects[i]->lab2Body(Tb); |
364 | ||
365 | // get the angular momentum, and propagate a half step | |
366 | ||
367 | < | dAtom->getJ(ji); |
367 | > | integrableObjects[i]->getJ(ji); |
368 | ||
369 | for (j = 0; j < 3; j++) | |
370 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
371 | ||
372 | < | this->rotationPropagation( dAtom, ji ); |
372 | > | this->rotationPropagation( integrableObjects[i], ji ); |
373 | ||
374 | < | dAtom->setJ(ji); |
374 | > | integrableObjects[i]->setJ(ji); |
375 | } | |
376 | } | |
377 | ||
# | Line 346 | Line 383 | template<typename T> void Integrator<T>::moveB(void){ | |
383 | ||
384 | template<typename T> void Integrator<T>::moveB(void){ | |
385 | int i, j; | |
349 | – | DirectionalAtom* dAtom; |
386 | double Tb[3], ji[3]; | |
387 | double vel[3], frc[3]; | |
388 | double mass; | |
389 | ||
390 | < | for (i = 0; i < nAtoms; i++){ |
391 | < | atoms[i]->getVel(vel); |
392 | < | atoms[i]->getFrc(frc); |
390 | > | for (i = 0; i < integrableObjects.size(); i++){ |
391 | > | integrableObjects[i]->getVel(vel); |
392 | > | integrableObjects[i]->getFrc(frc); |
393 | ||
394 | < | mass = atoms[i]->getMass(); |
394 | > | mass = integrableObjects[i]->getMass(); |
395 | ||
396 | // velocity half step | |
397 | for (j = 0; j < 3; j++) | |
398 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
399 | ||
400 | < | atoms[i]->setVel(vel); |
400 | > | integrableObjects[i]->setVel(vel); |
401 | ||
402 | < | if (atoms[i]->isDirectional()){ |
367 | < | dAtom = (DirectionalAtom *) atoms[i]; |
402 | > | if (integrableObjects[i]->isDirectional()){ |
403 | ||
404 | < | // get and convert the torque to body frame |
404 | > | // get and convert the torque to body frame |
405 | ||
406 | < | dAtom->getTrq(Tb); |
407 | < | dAtom->lab2Body(Tb); |
406 | > | integrableObjects[i]->getTrq(Tb); |
407 | > | integrableObjects[i]->lab2Body(Tb); |
408 | ||
409 | // get the angular momentum, and propagate a half step | |
410 | ||
411 | < | dAtom->getJ(ji); |
411 | > | integrableObjects[i]->getJ(ji); |
412 | ||
413 | for (j = 0; j < 3; j++) | |
414 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
415 | ||
416 | ||
417 | < | dAtom->setJ(ji); |
417 | > | integrableObjects[i]->setJ(ji); |
418 | } | |
419 | } | |
420 | ||
# | Line 404 | Line 439 | template<typename T> void Integrator<T>::constrainA(){ | |
439 | } | |
440 | ||
441 | template<typename T> void Integrator<T>::constrainA(){ | |
442 | < | int i, j, k; |
442 | > | int i, j; |
443 | int done; | |
444 | double posA[3], posB[3]; | |
445 | double velA[3], velB[3]; | |
# | Line 548 | Line 583 | template<typename T> void Integrator<T>::constrainB(vo | |
583 | } | |
584 | ||
585 | template<typename T> void Integrator<T>::constrainB(void){ | |
586 | < | int i, j, k; |
586 | > | int i, j; |
587 | int done; | |
588 | double posA[3], posB[3]; | |
589 | double velA[3], velB[3]; | |
# | Line 557 | Line 592 | template<typename T> void Integrator<T>::constrainB(vo | |
592 | int a, b, ax, ay, az, bx, by, bz; | |
593 | double rma, rmb; | |
594 | double dx, dy, dz; | |
595 | < | double rabsq, pabsq, rvab; |
561 | < | double diffsq; |
595 | > | double rvab; |
596 | double gab; | |
597 | int iteration; | |
598 | ||
# | Line 649 | Line 683 | template<typename T> void Integrator<T>::rotationPropa | |
683 | } | |
684 | ||
685 | template<typename T> void Integrator<T>::rotationPropagation | |
686 | < | ( DirectionalAtom* dAtom, double ji[3] ){ |
686 | > | ( StuntDouble* sd, double ji[3] ){ |
687 | ||
688 | double angle; | |
689 | double A[3][3], I[3][3]; | |
# | Line 657 | Line 691 | template<typename T> void Integrator<T>::rotationPropa | |
691 | // use the angular velocities to propagate the rotation matrix a | |
692 | // full time step | |
693 | ||
694 | < | dAtom->getA(A); |
695 | < | dAtom->getI(I); |
696 | < | |
697 | < | // rotate about the x-axis |
694 | > | sd->getA(A); |
695 | > | sd->getI(I); |
696 | > | |
697 | > | // rotate about the x-axis |
698 | angle = dt2 * ji[0] / I[0][0]; | |
699 | < | this->rotate( 1, 2, angle, ji, A ); |
700 | < | |
699 | > | this->rotate( 1, 2, angle, ji, A ); |
700 | > | |
701 | // rotate about the y-axis | |
702 | angle = dt2 * ji[1] / I[1][1]; | |
703 | this->rotate( 2, 0, angle, ji, A ); | |
704 | < | |
704 | > | |
705 | // rotate about the z-axis | |
706 | angle = dt * ji[2] / I[2][2]; | |
707 | this->rotate( 0, 1, angle, ji, A); | |
708 | < | |
708 | > | |
709 | // rotate about the y-axis | |
710 | angle = dt2 * ji[1] / I[1][1]; | |
711 | this->rotate( 2, 0, angle, ji, A ); | |
712 | < | |
712 | > | |
713 | // rotate about the x-axis | |
714 | angle = dt2 * ji[0] / I[0][0]; | |
715 | this->rotate( 1, 2, angle, ji, A ); | |
716 | < | |
717 | < | dAtom->setA( A ); |
716 | > | |
717 | > | sd->setA( A ); |
718 | } | |
719 | ||
720 | template<typename T> void Integrator<T>::rotate(int axes1, int axes2, | |
# | Line 748 | Line 782 | template<typename T> void Integrator<T>::rotate(int ax | |
782 | } | |
783 | } | |
784 | ||
785 | < | // rotate the Rotation matrix acording to: |
785 | > | // rotate the Rotation matrix acording to: |
786 | // A[][] = A[][] * transpose(rot[][]) | |
787 | ||
788 | ||
# | Line 777 | Line 811 | template<typename T> double Integrator<T>::getConserve | |
811 | template<typename T> double Integrator<T>::getConservedQuantity(void){ | |
812 | return tStats->getTotalE(); | |
813 | } | |
814 | + | template<typename T> string Integrator<T>::getAdditionalParameters(void){ |
815 | + | //By default, return a null string |
816 | + | //The reason we use string instead of char* is that if we use char*, we will |
817 | + | //return a pointer point to local variable which might cause problem |
818 | + | return string(); |
819 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |