| 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 | 
  | 
 | 
| 29 | 
  | 
  if (info->the_integrator != NULL){ | 
| 30 | 
  | 
    delete info->the_integrator; | 
| 31 | 
  | 
  } | 
| 32 | 
< | 
   | 
| 32 | 
> | 
 | 
| 33 | 
  | 
  nAtoms = info->n_atoms; | 
| 34 | 
  | 
 | 
| 35 | 
  | 
  // check for constraints | 
| 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 | 
  | 
 | 
| 162 | 
  | 
  double currThermal; | 
| 163 | 
  | 
  double currStatus; | 
| 164 | 
  | 
  double currReset; | 
| 165 | 
< | 
   | 
| 165 | 
> | 
 | 
| 166 | 
  | 
  int calcPot, calcStress; | 
| 167 | 
  | 
 | 
| 168 | 
  | 
  tStats = new Thermo(info); | 
| 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 | 
  | 
   | 
| 205 | 
  | 
  statOut->writeStat(info->getTime()); | 
| 206 | 
  | 
 | 
| 207 | 
  | 
 | 
| 201 | 
– | 
 | 
| 208 | 
  | 
#ifdef IS_MPI | 
| 209 | 
  | 
  strcpy(checkPointMsg, "The integrator is ready to go."); | 
| 210 | 
  | 
  MPIcheckPoint(); | 
| 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 | 
  | 
 | 
| 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){ | 
| 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."); | 
| 266 | 
  | 
#endif // is_mpi | 
| 267 | 
  | 
  } | 
| 268 | 
  | 
 | 
| 249 | 
– | 
 | 
| 250 | 
– | 
  // write the last frame | 
| 251 | 
– | 
  dumpOut->writeDump(info->getTime()); | 
| 252 | 
– | 
 | 
| 269 | 
  | 
  delete dumpOut; | 
| 270 | 
  | 
  delete statOut; | 
| 271 | 
  | 
} | 
| 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 | 
| 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(); | 
| 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); | 
| 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, | 
| 785 | 
  | 
    } | 
| 786 | 
  | 
  } | 
| 787 | 
  | 
 | 
| 788 | 
< | 
  // rotate the Rotation matrix acording to:  | 
| 788 | 
> | 
  // rotate the Rotation matrix acording to: | 
| 789 | 
  | 
  //            A[][] = A[][] * transpose(rot[][]) | 
| 790 | 
  | 
 | 
| 791 | 
  | 
 | 
| 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 | 
+ | 
} |