| 11 | 
  | 
#include "simError.h" | 
| 12 | 
  | 
 | 
| 13 | 
  | 
 | 
| 14 | 
< | 
Integrator::Integrator( SimInfo *theInfo, ForceFields* the_ff ){ | 
| 14 | 
> | 
template<typename T> Integrator<T>::Integrator( SimInfo *theInfo, ForceFields* the_ff ) { | 
| 15 | 
  | 
   | 
| 16 | 
  | 
  info = theInfo; | 
| 17 | 
  | 
  myFF = the_ff; | 
| 27 | 
  | 
 | 
| 28 | 
  | 
  nAtoms = info->n_atoms; | 
| 29 | 
  | 
 | 
| 30 | 
– | 
  std::cerr << "integ nAtoms = "  << nAtoms << "\n"; | 
| 31 | 
– | 
 | 
| 30 | 
  | 
  // check for constraints | 
| 31 | 
  | 
   | 
| 32 | 
  | 
  constrainedA    = NULL; | 
| 41 | 
  | 
  checkConstraints(); | 
| 42 | 
  | 
} | 
| 43 | 
  | 
 | 
| 44 | 
< | 
Integrator::~Integrator() { | 
| 44 | 
> | 
template<typename T> Integrator<T>::~Integrator() { | 
| 45 | 
  | 
   | 
| 46 | 
  | 
  if( nConstrained ){ | 
| 47 | 
  | 
    delete[] constrainedA; | 
| 54 | 
  | 
   | 
| 55 | 
  | 
} | 
| 56 | 
  | 
 | 
| 57 | 
< | 
void Integrator::checkConstraints( void ){ | 
| 57 | 
> | 
template<typename T> void Integrator<T>::checkConstraints( void ){ | 
| 58 | 
  | 
 | 
| 59 | 
  | 
 | 
| 60 | 
  | 
  isConstrained = 0; | 
| 73 | 
  | 
       | 
| 74 | 
  | 
      constrained = theArray[j]->is_constrained(); | 
| 75 | 
  | 
 | 
| 78 | 
– | 
      std::cerr << "Is the folowing bond constrained \n"; | 
| 79 | 
– | 
      theArray[j]->printMe(); | 
| 80 | 
– | 
       | 
| 76 | 
  | 
      if(constrained){ | 
| 82 | 
– | 
         | 
| 83 | 
– | 
        std::cerr << "Yes\n"; | 
| 77 | 
  | 
 | 
| 78 | 
  | 
        dummy_plug = theArray[j]->get_constraint(); | 
| 79 | 
  | 
        temp_con[nConstrained].set_a( dummy_plug->get_a() ); | 
| 83 | 
  | 
        nConstrained++; | 
| 84 | 
  | 
        constrained = 0; | 
| 85 | 
  | 
      }  | 
| 93 | 
– | 
      else std::cerr << "No.\n"; | 
| 86 | 
  | 
    } | 
| 87 | 
  | 
 | 
| 88 | 
  | 
    theArray = (SRI**) molecules[i].getMyBends(); | 
| 155 | 
  | 
} | 
| 156 | 
  | 
 | 
| 157 | 
  | 
 | 
| 158 | 
< | 
void Integrator::integrate( void ){ | 
| 158 | 
> | 
template<typename T> void Integrator<T>::integrate( void ){ | 
| 159 | 
  | 
 | 
| 160 | 
  | 
  int i, j;                         // loop counters | 
| 161 | 
  | 
 | 
| 167 | 
  | 
  double currSample; | 
| 168 | 
  | 
  double currThermal; | 
| 169 | 
  | 
  double currStatus; | 
| 178 | 
– | 
  double currTime; | 
| 170 | 
  | 
 | 
| 171 | 
  | 
  int calcPot, calcStress; | 
| 172 | 
  | 
  int isError; | 
| 190 | 
  | 
    tStats->velocitize(); | 
| 191 | 
  | 
  } | 
| 192 | 
  | 
   | 
| 202 | 
– | 
  dumpOut->writeDump( 0.0 ); | 
| 203 | 
– | 
  statOut->writeStat( 0.0 ); | 
| 204 | 
– | 
   | 
| 193 | 
  | 
  calcPot     = 0; | 
| 194 | 
  | 
  calcStress  = 0; | 
| 195 | 
  | 
  currSample  = sampleTime; | 
| 196 | 
  | 
  currThermal = thermalTime; | 
| 197 | 
  | 
  currStatus  = statusTime; | 
| 210 | 
– | 
  currTime    = 0.0;; | 
| 198 | 
  | 
 | 
| 199 | 
+ | 
  dumpOut->writeDump( info->getTime() ); | 
| 200 | 
+ | 
  statOut->writeStat( info->getTime() ); | 
| 201 | 
  | 
 | 
| 202 | 
  | 
  readyCheck(); | 
| 203 | 
  | 
 | 
| 207 | 
  | 
  MPIcheckPoint(); | 
| 208 | 
  | 
#endif // is_mpi | 
| 209 | 
  | 
 | 
| 210 | 
< | 
  while( currTime < runTime ){ | 
| 210 | 
> | 
  while( info->getTime() < runTime ){ | 
| 211 | 
  | 
 | 
| 212 | 
< | 
    if( (currTime+dt) >= currStatus ){ | 
| 212 | 
> | 
    if( (info->getTime()+dt) >= currStatus ){ | 
| 213 | 
  | 
      calcPot = 1; | 
| 214 | 
  | 
      calcStress = 1; | 
| 215 | 
  | 
    } | 
| 216 | 
  | 
 | 
| 228 | 
– | 
    std::cerr << currTime << "\n"; | 
| 229 | 
– | 
 | 
| 217 | 
  | 
    integrateStep( calcPot, calcStress ); | 
| 218 | 
  | 
       | 
| 219 | 
< | 
    currTime += dt; | 
| 219 | 
> | 
    info->incrTime(dt); | 
| 220 | 
  | 
 | 
| 221 | 
  | 
    if( info->setTemp ){ | 
| 222 | 
< | 
      if( currTime >= currThermal ){ | 
| 222 | 
> | 
      if( info->getTime() >= currThermal ){ | 
| 223 | 
  | 
        tStats->velocitize(); | 
| 224 | 
  | 
        currThermal += thermalTime; | 
| 225 | 
  | 
      } | 
| 226 | 
  | 
    } | 
| 227 | 
  | 
 | 
| 228 | 
< | 
    if( currTime >= currSample ){ | 
| 229 | 
< | 
      dumpOut->writeDump( currTime ); | 
| 228 | 
> | 
    if( info->getTime() >= currSample ){ | 
| 229 | 
> | 
      dumpOut->writeDump( info->getTime() ); | 
| 230 | 
  | 
      currSample += sampleTime; | 
| 231 | 
  | 
    } | 
| 232 | 
  | 
 | 
| 233 | 
< | 
    if( currTime >= currStatus ){  | 
| 234 | 
< | 
      statOut->writeStat( currTime );  | 
| 233 | 
> | 
    if( info->getTime() >= currStatus ){  | 
| 234 | 
> | 
      statOut->writeStat( info->getTime() );  | 
| 235 | 
  | 
      calcPot = 0;  | 
| 236 | 
  | 
      calcStress = 0; | 
| 237 | 
  | 
      currStatus += statusTime; | 
| 245 | 
  | 
 | 
| 246 | 
  | 
  } | 
| 247 | 
  | 
 | 
| 248 | 
< | 
  dumpOut->writeFinal(currTime); | 
| 248 | 
> | 
  dumpOut->writeFinal(info->getTime()); | 
| 249 | 
  | 
 | 
| 250 | 
  | 
  delete dumpOut; | 
| 251 | 
  | 
  delete statOut; | 
| 252 | 
  | 
} | 
| 253 | 
  | 
 | 
| 254 | 
< | 
void Integrator::integrateStep( int calcPot, int calcStress ){ | 
| 254 | 
> | 
template<typename T> void Integrator<T>::integrateStep( int calcPot, int calcStress ){ | 
| 255 | 
  | 
 | 
| 256 | 
  | 
 | 
| 257 | 
  | 
       | 
| 261 | 
  | 
  moveA(); | 
| 262 | 
  | 
  if( nConstrained ) constrainA(); | 
| 263 | 
  | 
 | 
| 264 | 
+ | 
   | 
| 265 | 
+ | 
#ifdef IS_MPI | 
| 266 | 
+ | 
  strcpy( checkPointMsg, "Succesful moveA\n" ); | 
| 267 | 
+ | 
  MPIcheckPoint(); | 
| 268 | 
+ | 
#endif // is_mpi | 
| 269 | 
+ | 
   | 
| 270 | 
+ | 
 | 
| 271 | 
  | 
  // calc forces | 
| 272 | 
  | 
 | 
| 273 | 
  | 
  myFF->doForces(calcPot,calcStress); | 
| 274 | 
  | 
 | 
| 275 | 
+ | 
#ifdef IS_MPI | 
| 276 | 
+ | 
  strcpy( checkPointMsg, "Succesful doForces\n" ); | 
| 277 | 
+ | 
  MPIcheckPoint(); | 
| 278 | 
+ | 
#endif // is_mpi | 
| 279 | 
+ | 
   | 
| 280 | 
+ | 
 | 
| 281 | 
  | 
  // finish the velocity  half step | 
| 282 | 
  | 
   | 
| 283 | 
  | 
  moveB(); | 
| 284 | 
  | 
  if( nConstrained ) constrainB(); | 
| 285 | 
< | 
    | 
| 285 | 
> | 
   | 
| 286 | 
> | 
#ifdef IS_MPI | 
| 287 | 
> | 
  strcpy( checkPointMsg, "Succesful moveB\n" ); | 
| 288 | 
> | 
  MPIcheckPoint(); | 
| 289 | 
> | 
#endif // is_mpi | 
| 290 | 
> | 
   | 
| 291 | 
> | 
  | 
| 292 | 
  | 
} | 
| 293 | 
  | 
 | 
| 294 | 
  | 
 | 
| 295 | 
< | 
void Integrator::moveA( void ){ | 
| 295 | 
> | 
template<typename T> void Integrator<T>::moveA( void ){ | 
| 296 | 
  | 
   | 
| 297 | 
  | 
  int i, j; | 
| 298 | 
  | 
  DirectionalAtom* dAtom; | 
| 371 | 
  | 
} | 
| 372 | 
  | 
 | 
| 373 | 
  | 
 | 
| 374 | 
< | 
void Integrator::moveB( void ){ | 
| 374 | 
> | 
template<typename T> void Integrator<T>::moveB( void ){ | 
| 375 | 
  | 
  int i, j; | 
| 376 | 
  | 
  DirectionalAtom* dAtom; | 
| 377 | 
  | 
  double Tb[3], ji[3]; | 
| 413 | 
  | 
  } | 
| 414 | 
  | 
} | 
| 415 | 
  | 
 | 
| 416 | 
< | 
void Integrator::preMove( void ){ | 
| 416 | 
> | 
template<typename T> void Integrator<T>::preMove( void ){ | 
| 417 | 
  | 
  int i, j; | 
| 418 | 
  | 
  double pos[3]; | 
| 419 | 
  | 
 | 
| 431 | 
  | 
  }   | 
| 432 | 
  | 
} | 
| 433 | 
  | 
 | 
| 434 | 
< | 
void Integrator::constrainA(){ | 
| 434 | 
> | 
template<typename T> void Integrator<T>::constrainA(){ | 
| 435 | 
  | 
 | 
| 436 | 
  | 
  int i,j,k; | 
| 437 | 
  | 
  int done; | 
| 582 | 
  | 
 | 
| 583 | 
  | 
} | 
| 584 | 
  | 
 | 
| 585 | 
< | 
void Integrator::constrainB( void ){ | 
| 585 | 
> | 
template<typename T> void Integrator<T>::constrainB( void ){ | 
| 586 | 
  | 
   | 
| 587 | 
  | 
  int i,j,k; | 
| 588 | 
  | 
  int done; | 
| 691 | 
  | 
 | 
| 692 | 
  | 
} | 
| 693 | 
  | 
 | 
| 694 | 
< | 
void Integrator::rotate( int axes1, int axes2, double angle, double ji[3],  | 
| 694 | 
> | 
template<typename T> void Integrator<T>::rotate( int axes1, int axes2, double angle, double ji[3],  | 
| 695 | 
  | 
                         double A[3][3] ){ | 
| 696 | 
  | 
 | 
| 697 | 
  | 
  int i,j,k; |