# | 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 | ||
# | Line 132 | Line 135 | template<typename T> void Integrator<T>::checkConstrai | |
135 | } | |
136 | ||
137 | ||
138 | < | // save oldAtoms to check for lode balanceing later on. |
138 | > | // save oldAtoms to check for lode balancing later on. |
139 | ||
140 | oldAtoms = nAtoms; | |
141 | ||
# | Line 147 | Line 150 | template<typename T> void Integrator<T>::integrate(voi | |
150 | ||
151 | ||
152 | template<typename T> void Integrator<T>::integrate(void){ | |
150 | – | int i, j; // loop counters |
153 | ||
154 | double runTime = info->run_time; | |
155 | double sampleTime = info->sampleTime; | |
# | Line 160 | Line 162 | template<typename T> void Integrator<T>::integrate(voi | |
162 | double currThermal; | |
163 | double currStatus; | |
164 | double currReset; | |
165 | < | |
165 | > | |
166 | int calcPot, calcStress; | |
165 | – | int isError; |
167 | ||
168 | tStats = new Thermo(info); | |
169 | statOut = new StatWriter(info); | |
170 | dumpOut = new DumpWriter(info); | |
171 | ||
172 | atoms = info->atoms; | |
172 | – | DirectionalAtom* dAtom; |
173 | ||
174 | dt = info->dt; | |
175 | dt2 = 0.5 * dt; | |
176 | ||
177 | + | readyCheck(); |
178 | + | |
179 | // initialize the forces before the first step | |
180 | ||
181 | calcForce(1, 1); | |
182 | ||
183 | + | //temp test |
184 | + | tStats->getPotential(); |
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."); | |
# | Line 212 | Line 216 | template<typename T> void Integrator<T>::integrate(voi | |
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 | < | |
323 | > | #ifdef PROFILE |
324 | > | endProfile(pro6); |
325 | > | #endif // profile |
326 | ||
327 | #ifdef IS_MPI | |
328 | strcpy(checkPointMsg, "Succesful moveB\n"); | |
# | Line 366 | Line 404 | template<typename T> void Integrator<T>::moveB(void){ | |
404 | if (atoms[i]->isDirectional()){ | |
405 | dAtom = (DirectionalAtom *) atoms[i]; | |
406 | ||
407 | < | // get and convert the torque to body frame |
407 | > | // get and convert the torque to body frame |
408 | ||
409 | dAtom->getTrq(Tb); | |
410 | dAtom->lab2Body(Tb); | |
# | Line 404 | Line 442 | template<typename T> void Integrator<T>::constrainA(){ | |
442 | } | |
443 | ||
444 | template<typename T> void Integrator<T>::constrainA(){ | |
445 | < | int i, j, k; |
445 | > | int i, j; |
446 | int done; | |
447 | double posA[3], posB[3]; | |
448 | double velA[3], velB[3]; | |
# | Line 548 | Line 586 | template<typename T> void Integrator<T>::constrainB(vo | |
586 | } | |
587 | ||
588 | template<typename T> void Integrator<T>::constrainB(void){ | |
589 | < | int i, j, k; |
589 | > | int i, j; |
590 | int done; | |
591 | double posA[3], posB[3]; | |
592 | double velA[3], velB[3]; | |
# | Line 557 | Line 595 | template<typename T> void Integrator<T>::constrainB(vo | |
595 | int a, b, ax, ay, az, bx, by, bz; | |
596 | double rma, rmb; | |
597 | double dx, dy, dz; | |
598 | < | double rabsq, pabsq, rvab; |
561 | < | double diffsq; |
598 | > | double rvab; |
599 | double gab; | |
600 | int iteration; | |
601 | ||
# | Line 659 | Line 696 | template<typename T> void Integrator<T>::rotationPropa | |
696 | ||
697 | dAtom->getA(A); | |
698 | dAtom->getI(I); | |
699 | < | |
700 | < | // rotate about the x-axis |
699 | > | |
700 | > | // rotate about the x-axis |
701 | angle = dt2 * ji[0] / I[0][0]; | |
702 | < | this->rotate( 1, 2, angle, ji, A ); |
703 | < | |
702 | > | this->rotate( 1, 2, angle, ji, A ); |
703 | > | |
704 | // rotate about the y-axis | |
705 | angle = dt2 * ji[1] / I[1][1]; | |
706 | this->rotate( 2, 0, angle, ji, A ); | |
707 | < | |
707 | > | |
708 | // rotate about the z-axis | |
709 | angle = dt * ji[2] / I[2][2]; | |
710 | this->rotate( 0, 1, angle, ji, A); | |
711 | < | |
711 | > | |
712 | // rotate about the y-axis | |
713 | angle = dt2 * ji[1] / I[1][1]; | |
714 | this->rotate( 2, 0, angle, ji, A ); | |
715 | < | |
715 | > | |
716 | // rotate about the x-axis | |
717 | angle = dt2 * ji[0] / I[0][0]; | |
718 | this->rotate( 1, 2, angle, ji, A ); | |
719 | < | |
720 | < | dAtom->setA( A ); |
719 | > | |
720 | > | dAtom->setA( A ); |
721 | } | |
722 | ||
723 | template<typename T> void Integrator<T>::rotate(int axes1, int axes2, | |
# | Line 748 | Line 785 | template<typename T> void Integrator<T>::rotate(int ax | |
785 | } | |
786 | } | |
787 | ||
788 | < | // rotate the Rotation matrix acording to: |
788 | > | // rotate the Rotation matrix acording to: |
789 | // A[][] = A[][] * transpose(rot[][]) | |
790 | ||
791 | ||
# | Line 777 | Line 814 | template<typename T> double Integrator<T>::getConserve | |
814 | template<typename T> double Integrator<T>::getConservedQuantity(void){ | |
815 | return tStats->getTotalE(); | |
816 | } | |
817 | + | template<typename T> string Integrator<T>::getAdditionalParameters(void){ |
818 | + | //By default, return a null string |
819 | + | //The reason we use string instead of char* is that if we use char*, we will |
820 | + | //return a pointer point to local variable which might cause problem |
821 | + | return string(); |
822 | + | } |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |