# | 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 | } | |
32 | < | |
32 | > | |
33 | nAtoms = info->n_atoms; | |
34 | + | integrableObjects = info->integrableObjects; |
35 | ||
36 | // check for constraints | |
37 | ||
# | Line 64 | 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 110 | Line 116 | template<typename T> void Integrator<T>::checkConstrai | |
116 | } | |
117 | } | |
118 | ||
119 | + | |
120 | if (nConstrained > 0){ | |
121 | isConstrained = 1; | |
122 | ||
# | Line 131 | 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 146 | Line 153 | template<typename T> void Integrator<T>::integrate(voi | |
153 | ||
154 | ||
155 | template<typename T> void Integrator<T>::integrate(void){ | |
149 | – | int i, j; // loop counters |
156 | ||
157 | double runTime = info->run_time; | |
158 | double sampleTime = info->sampleTime; | |
# | Line 159 | 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; | |
164 | – | int isError; |
170 | ||
171 | tStats = new Thermo(info); | |
172 | statOut = new StatWriter(info); | |
173 | dumpOut = new DumpWriter(info); | |
174 | ||
175 | atoms = info->atoms; | |
171 | – | DirectionalAtom* dAtom; |
176 | ||
177 | dt = info->dt; | |
178 | dt2 = 0.5 * dt; | |
179 | ||
180 | + | readyCheck(); |
181 | + | |
182 | + | // remove center of mass drift velocity (in case we passed in a configuration |
183 | + | // that was drifting |
184 | + | tStats->removeCOMdrift(); |
185 | + | |
186 | // initialize the forces before the first step | |
187 | ||
188 | calcForce(1, 1); | |
189 | < | |
189 | > | |
190 | if (nConstrained){ | |
191 | preMove(); | |
192 | constrainA(); | |
193 | < | calcForce(1, 1); |
193 | > | calcForce(1, 1); |
194 | constrainB(); | |
195 | } | |
196 | ||
# | Line 198 | Line 208 | template<typename T> void Integrator<T>::integrate(voi | |
208 | dumpOut->writeDump(info->getTime()); | |
209 | statOut->writeStat(info->getTime()); | |
210 | ||
201 | – | readyCheck(); |
211 | ||
212 | #ifdef IS_MPI | |
213 | strcpy(checkPointMsg, "The integrator is ready to go."); | |
214 | MPIcheckPoint(); | |
215 | #endif // is_mpi | |
216 | ||
217 | < | while (info->getTime() < runTime){ |
217 | > | while (info->getTime() < runTime && !stopIntegrator()){ |
218 | if ((info->getTime() + dt) >= currStatus){ | |
219 | calcPot = 1; | |
220 | calcStress = 1; | |
221 | } | |
222 | ||
223 | + | #ifdef PROFILE |
224 | + | startProfile( pro1 ); |
225 | + | #endif |
226 | + | |
227 | integrateStep(calcPot, calcStress); | |
228 | ||
229 | + | #ifdef PROFILE |
230 | + | endProfile( pro1 ); |
231 | + | |
232 | + | startProfile( pro2 ); |
233 | + | #endif // profile |
234 | + | |
235 | info->incrTime(dt); | |
236 | ||
237 | if (info->setTemp){ | |
# | Line 228 | Line 247 | template<typename T> void Integrator<T>::integrate(voi | |
247 | } | |
248 | ||
249 | if (info->getTime() >= currStatus){ | |
250 | < | statOut->writeStat(info->getTime()); |
251 | < | calcPot = 0; |
250 | > | statOut->writeStat(info->getTime()); |
251 | > | calcPot = 0; |
252 | calcStress = 0; | |
253 | currStatus += statusTime; | |
254 | < | } |
254 | > | } |
255 | ||
256 | if (info->resetIntegrator){ | |
257 | if (info->getTime() >= currReset){ | |
# | Line 240 | Line 259 | template<typename T> void Integrator<T>::integrate(voi | |
259 | currReset += resetTime; | |
260 | } | |
261 | } | |
262 | + | |
263 | + | #ifdef PROFILE |
264 | + | endProfile( pro2 ); |
265 | + | #endif //profile |
266 | ||
267 | #ifdef IS_MPI | |
268 | strcpy(checkPointMsg, "successfully took a time step."); | |
# | Line 247 | Line 270 | template<typename T> void Integrator<T>::integrate(voi | |
270 | #endif // is_mpi | |
271 | } | |
272 | ||
250 | – | dumpOut->writeFinal(info->getTime()); |
251 | – | |
273 | delete dumpOut; | |
274 | delete statOut; | |
275 | } | |
# | Line 256 | Line 277 | template<typename T> void Integrator<T>::integrateStep | |
277 | template<typename T> void Integrator<T>::integrateStep(int calcPot, | |
278 | int calcStress){ | |
279 | // Position full step, and velocity half step | |
280 | + | |
281 | + | #ifdef PROFILE |
282 | + | startProfile(pro3); |
283 | + | #endif //profile |
284 | + | |
285 | preMove(); | |
286 | ||
287 | < | moveA(); |
287 | > | #ifdef PROFILE |
288 | > | endProfile(pro3); |
289 | ||
290 | + | startProfile(pro4); |
291 | + | #endif // profile |
292 | ||
293 | + | moveA(); |
294 | ||
295 | + | #ifdef PROFILE |
296 | + | endProfile(pro4); |
297 | + | |
298 | + | startProfile(pro5); |
299 | + | #endif//profile |
300 | ||
301 | + | |
302 | #ifdef IS_MPI | |
303 | strcpy(checkPointMsg, "Succesful moveA\n"); | |
304 | MPIcheckPoint(); | |
# | Line 278 | Line 314 | template<typename T> void Integrator<T>::integrateStep | |
314 | MPIcheckPoint(); | |
315 | #endif // is_mpi | |
316 | ||
317 | + | #ifdef PROFILE |
318 | + | endProfile( pro5 ); |
319 | ||
320 | + | startProfile( pro6 ); |
321 | + | #endif //profile |
322 | + | |
323 | // finish the velocity half step | |
324 | ||
325 | moveB(); | |
326 | ||
327 | + | #ifdef PROFILE |
328 | + | endProfile(pro6); |
329 | + | #endif // profile |
330 | ||
287 | – | |
331 | #ifdef IS_MPI | |
332 | strcpy(checkPointMsg, "Succesful moveB\n"); | |
333 | MPIcheckPoint(); | |
# | Line 293 | Line 336 | template<typename T> void Integrator<T>::moveA(void){ | |
336 | ||
337 | ||
338 | template<typename T> void Integrator<T>::moveA(void){ | |
339 | < | int i, j; |
339 | > | size_t i, j; |
340 | DirectionalAtom* dAtom; | |
341 | double Tb[3], ji[3]; | |
342 | double vel[3], pos[3], frc[3]; | |
343 | double mass; | |
344 | + | |
345 | + | for (i = 0; i < integrableObjects.size() ; i++){ |
346 | + | integrableObjects[i]->getVel(vel); |
347 | + | integrableObjects[i]->getPos(pos); |
348 | + | integrableObjects[i]->getFrc(frc); |
349 | ||
350 | < | for (i = 0; i < nAtoms; i++){ |
351 | < | atoms[i]->getVel(vel); |
352 | < | atoms[i]->getPos(pos); |
305 | < | atoms[i]->getFrc(frc); |
350 | > | std::cerr << "i =\t" << i << "\t" << frc[0] << "\t" << frc[1]<< "\t" << frc[2] << "\n"; |
351 | > | |
352 | > | mass = integrableObjects[i]->getMass(); |
353 | ||
307 | – | mass = atoms[i]->getMass(); |
308 | – | |
354 | for (j = 0; j < 3; j++){ | |
355 | // velocity half step | |
356 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
# | Line 313 | Line 358 | template<typename T> void Integrator<T>::moveA(void){ | |
358 | pos[j] += dt * vel[j]; | |
359 | } | |
360 | ||
361 | < | atoms[i]->setVel(vel); |
362 | < | atoms[i]->setPos(pos); |
361 | > | integrableObjects[i]->setVel(vel); |
362 | > | integrableObjects[i]->setPos(pos); |
363 | ||
364 | < | if (atoms[i]->isDirectional()){ |
320 | < | dAtom = (DirectionalAtom *) atoms[i]; |
364 | > | if (integrableObjects[i]->isDirectional()){ |
365 | ||
366 | // get and convert the torque to body frame | |
367 | ||
368 | < | dAtom->getTrq(Tb); |
369 | < | dAtom->lab2Body(Tb); |
368 | > | integrableObjects[i]->getTrq(Tb); |
369 | > | integrableObjects[i]->lab2Body(Tb); |
370 | ||
371 | // get the angular momentum, and propagate a half step | |
372 | ||
373 | < | dAtom->getJ(ji); |
373 | > | integrableObjects[i]->getJ(ji); |
374 | ||
375 | for (j = 0; j < 3; j++) | |
376 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
377 | ||
378 | < | this->rotationPropagation( dAtom, ji ); |
378 | > | this->rotationPropagation( integrableObjects[i], ji ); |
379 | ||
380 | < | dAtom->setJ(ji); |
380 | > | integrableObjects[i]->setJ(ji); |
381 | } | |
382 | } | |
383 | ||
# | Line 345 | Line 389 | template<typename T> void Integrator<T>::moveB(void){ | |
389 | ||
390 | template<typename T> void Integrator<T>::moveB(void){ | |
391 | int i, j; | |
348 | – | DirectionalAtom* dAtom; |
392 | double Tb[3], ji[3]; | |
393 | double vel[3], frc[3]; | |
394 | double mass; | |
395 | ||
396 | < | for (i = 0; i < nAtoms; i++){ |
397 | < | atoms[i]->getVel(vel); |
398 | < | atoms[i]->getFrc(frc); |
396 | > | for (i = 0; i < integrableObjects.size(); i++){ |
397 | > | integrableObjects[i]->getVel(vel); |
398 | > | integrableObjects[i]->getFrc(frc); |
399 | ||
400 | < | mass = atoms[i]->getMass(); |
400 | > | mass = integrableObjects[i]->getMass(); |
401 | ||
402 | // velocity half step | |
403 | for (j = 0; j < 3; j++) | |
404 | vel[j] += (dt2 * frc[j] / mass) * eConvert; | |
405 | ||
406 | < | atoms[i]->setVel(vel); |
406 | > | integrableObjects[i]->setVel(vel); |
407 | ||
408 | < | if (atoms[i]->isDirectional()){ |
366 | < | dAtom = (DirectionalAtom *) atoms[i]; |
408 | > | if (integrableObjects[i]->isDirectional()){ |
409 | ||
410 | < | // get and convert the torque to body frame |
410 | > | // get and convert the torque to body frame |
411 | ||
412 | < | dAtom->getTrq(Tb); |
413 | < | dAtom->lab2Body(Tb); |
412 | > | integrableObjects[i]->getTrq(Tb); |
413 | > | integrableObjects[i]->lab2Body(Tb); |
414 | ||
415 | // get the angular momentum, and propagate a half step | |
416 | ||
417 | < | dAtom->getJ(ji); |
417 | > | integrableObjects[i]->getJ(ji); |
418 | ||
419 | for (j = 0; j < 3; j++) | |
420 | ji[j] += (dt2 * Tb[j]) * eConvert; | |
421 | ||
422 | ||
423 | < | dAtom->setJ(ji); |
423 | > | integrableObjects[i]->setJ(ji); |
424 | } | |
425 | } | |
426 | ||
# | Line 403 | Line 445 | template<typename T> void Integrator<T>::constrainA(){ | |
445 | } | |
446 | ||
447 | template<typename T> void Integrator<T>::constrainA(){ | |
448 | < | int i, j, k; |
448 | > | int i, j; |
449 | int done; | |
450 | double posA[3], posB[3]; | |
451 | double velA[3], velB[3]; | |
# | Line 547 | Line 589 | template<typename T> void Integrator<T>::constrainB(vo | |
589 | } | |
590 | ||
591 | template<typename T> void Integrator<T>::constrainB(void){ | |
592 | < | int i, j, k; |
592 | > | int i, j; |
593 | int done; | |
594 | double posA[3], posB[3]; | |
595 | double velA[3], velB[3]; | |
# | Line 556 | Line 598 | template<typename T> void Integrator<T>::constrainB(vo | |
598 | int a, b, ax, ay, az, bx, by, bz; | |
599 | double rma, rmb; | |
600 | double dx, dy, dz; | |
601 | < | double rabsq, pabsq, rvab; |
560 | < | double diffsq; |
601 | > | double rvab; |
602 | double gab; | |
603 | int iteration; | |
604 | ||
# | Line 648 | Line 689 | template<typename T> void Integrator<T>::rotationPropa | |
689 | } | |
690 | ||
691 | template<typename T> void Integrator<T>::rotationPropagation | |
692 | < | ( DirectionalAtom* dAtom, double ji[3] ){ |
692 | > | ( StuntDouble* sd, double ji[3] ){ |
693 | ||
694 | double angle; | |
695 | double A[3][3], I[3][3]; | |
696 | + | int i, j, k; |
697 | ||
698 | // use the angular velocities to propagate the rotation matrix a | |
699 | // full time step | |
700 | ||
701 | < | dAtom->getA(A); |
702 | < | dAtom->getI(I); |
703 | < | |
704 | < | // rotate about the x-axis |
705 | < | angle = dt2 * ji[0] / I[0][0]; |
706 | < | this->rotate( 1, 2, angle, ji, A ); |
707 | < | |
708 | < | // rotate about the y-axis |
709 | < | angle = dt2 * ji[1] / I[1][1]; |
710 | < | this->rotate( 2, 0, angle, ji, A ); |
711 | < | |
712 | < | // rotate about the z-axis |
713 | < | angle = dt * ji[2] / I[2][2]; |
714 | < | this->rotate( 0, 1, angle, ji, A); |
715 | < | |
716 | < | // rotate about the y-axis |
717 | < | angle = dt2 * ji[1] / I[1][1]; |
718 | < | this->rotate( 2, 0, angle, ji, A ); |
719 | < | |
720 | < | // rotate about the x-axis |
721 | < | angle = dt2 * ji[0] / I[0][0]; |
722 | < | this->rotate( 1, 2, angle, ji, A ); |
723 | < | |
724 | < | dAtom->setA( A ); |
701 | > | sd->getA(A); |
702 | > | sd->getI(I); |
703 | > | |
704 | > | if (sd->isLinear()) { |
705 | > | i = sd->linearAxis(); |
706 | > | j = (i+1)%3; |
707 | > | k = (i+2)%3; |
708 | > | |
709 | > | angle = dt2 * ji[j] / I[j][j]; |
710 | > | this->rotate( k, i, angle, ji, A ); |
711 | > | |
712 | > | angle = dt * ji[k] / I[k][k]; |
713 | > | this->rotate( i, j, angle, ji, A); |
714 | > | |
715 | > | angle = dt2 * ji[j] / I[j][j]; |
716 | > | this->rotate( k, i, angle, ji, A ); |
717 | > | |
718 | > | } else { |
719 | > | // rotate about the x-axis |
720 | > | angle = dt2 * ji[0] / I[0][0]; |
721 | > | this->rotate( 1, 2, angle, ji, A ); |
722 | > | |
723 | > | // rotate about the y-axis |
724 | > | angle = dt2 * ji[1] / I[1][1]; |
725 | > | this->rotate( 2, 0, angle, ji, A ); |
726 | > | |
727 | > | // rotate about the z-axis |
728 | > | angle = dt * ji[2] / I[2][2]; |
729 | > | this->rotate( 0, 1, angle, ji, A); |
730 | > | |
731 | > | // rotate about the y-axis |
732 | > | angle = dt2 * ji[1] / I[1][1]; |
733 | > | this->rotate( 2, 0, angle, ji, A ); |
734 | > | |
735 | > | // rotate about the x-axis |
736 | > | angle = dt2 * ji[0] / I[0][0]; |
737 | > | this->rotate( 1, 2, angle, ji, A ); |
738 | > | |
739 | > | } |
740 | > | sd->setA( A ); |
741 | } | |
742 | ||
743 | template<typename T> void Integrator<T>::rotate(int axes1, int axes2, | |
# | Line 747 | Line 805 | template<typename T> void Integrator<T>::rotate(int ax | |
805 | } | |
806 | } | |
807 | ||
808 | < | // rotate the Rotation matrix acording to: |
808 | > | // rotate the Rotation matrix acording to: |
809 | // A[][] = A[][] * transpose(rot[][]) | |
810 | ||
811 | ||
# | Line 776 | Line 834 | template<typename T> double Integrator<T>::getConserve | |
834 | template<typename T> double Integrator<T>::getConservedQuantity(void){ | |
835 | return tStats->getTotalE(); | |
836 | } | |
837 | + | template<typename T> string Integrator<T>::getAdditionalParameters(void){ |
838 | + | //By default, return a null string |
839 | + | //The reason we use string instead of char* is that if we use char*, we will |
840 | + | //return a pointer point to local variable which might cause problem |
841 | + | return string(); |
842 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |