| 1 |
|
#include <iostream> |
| 2 |
|
#include <stdlib.h> |
| 3 |
|
#include <math.h> |
| 4 |
< |
|
| 4 |
> |
#include "Rattle.hpp" |
| 5 |
|
#ifdef IS_MPI |
| 6 |
|
#include "mpiSimulation.hpp" |
| 7 |
|
#include <unistd.h> |
| 32 |
|
|
| 33 |
|
nAtoms = info->n_atoms; |
| 34 |
|
integrableObjects = info->integrableObjects; |
| 35 |
< |
|
| 35 |
> |
|
| 36 |
> |
rattle = new RattleFramework(info); |
| 37 |
> |
|
| 38 |
> |
if(rattle == NULL){ |
| 39 |
> |
sprintf(painCave.errMsg, |
| 40 |
> |
"Integrator::Intergrator() Error: Memory allocation error for RattleFramework" ); |
| 41 |
> |
painCave.isFatal = 1; |
| 42 |
> |
simError(); |
| 43 |
> |
} |
| 44 |
> |
|
| 45 |
> |
/* |
| 46 |
|
// check for constraints |
| 47 |
|
|
| 48 |
|
constrainedA = NULL; |
| 55 |
|
nConstrained = 0; |
| 56 |
|
|
| 57 |
|
checkConstraints(); |
| 58 |
< |
|
| 58 |
> |
*/ |
| 59 |
|
} |
| 60 |
|
|
| 61 |
|
template<typename T> Integrator<T>::~Integrator(){ |
| 62 |
+ |
if (rattle != NULL) |
| 63 |
+ |
delete rattle; |
| 64 |
+ |
/* |
| 65 |
|
if (nConstrained){ |
| 66 |
|
delete[] constrainedA; |
| 67 |
|
delete[] constrainedB; |
| 70 |
|
delete[] moved; |
| 71 |
|
delete[] oldPos; |
| 72 |
|
} |
| 73 |
+ |
*/ |
| 74 |
|
} |
| 75 |
|
|
| 76 |
+ |
/* |
| 77 |
|
template<typename T> void Integrator<T>::checkConstraints(void){ |
| 78 |
|
isConstrained = 0; |
| 79 |
|
|
| 166 |
|
|
| 167 |
|
delete[] temp_con; |
| 168 |
|
} |
| 169 |
+ |
*/ |
| 170 |
|
|
| 155 |
– |
|
| 171 |
|
template<typename T> void Integrator<T>::integrate(void){ |
| 172 |
|
|
| 173 |
|
double runTime = info->run_time; |
| 200 |
|
tStats->removeCOMdrift(); |
| 201 |
|
|
| 202 |
|
// initialize the retraints if necessary |
| 203 |
< |
if (info->useThermInt) { |
| 203 |
> |
if (info->useSolidThermInt && !info->useLiquidThermInt) { |
| 204 |
|
myFF->initRestraints(); |
| 205 |
|
} |
| 206 |
|
|
| 207 |
|
// initialize the forces before the first step |
| 208 |
|
|
| 209 |
|
calcForce(1, 1); |
| 210 |
< |
|
| 211 |
< |
if (nConstrained){ |
| 212 |
< |
preMove(); |
| 213 |
< |
constrainA(); |
| 214 |
< |
calcForce(1, 1); |
| 215 |
< |
constrainB(); |
| 201 |
< |
} |
| 210 |
> |
|
| 211 |
> |
//execute constraint algorithm to make sure at the very beginning the system is constrained |
| 212 |
> |
rattle->doPreConstraint(); |
| 213 |
> |
rattle->doRattleA(); |
| 214 |
> |
calcForce(1, 1); |
| 215 |
> |
rattle->doRattleB(); |
| 216 |
|
|
| 217 |
|
if (info->setTemp){ |
| 218 |
|
thermalize(); |
| 269 |
|
|
| 270 |
|
if (info->getTime() >= currStatus){ |
| 271 |
|
statOut->writeStat(info->getTime()); |
| 258 |
– |
if (info->useThermInt) |
| 259 |
– |
statOut->writeRaw(info->getTime()); |
| 272 |
|
calcPot = 0; |
| 273 |
|
calcStress = 0; |
| 274 |
|
currStatus += statusTime; |
| 292 |
|
} |
| 293 |
|
|
| 294 |
|
// dump out a file containing the omega values for the final configuration |
| 295 |
< |
if (info->useThermInt) |
| 295 |
> |
if (info->useSolidThermInt && !info->useLiquidThermInt) |
| 296 |
|
myFF->dumpzAngle(); |
| 297 |
|
|
| 298 |
|
|
| 308 |
|
startProfile(pro3); |
| 309 |
|
#endif //profile |
| 310 |
|
|
| 311 |
< |
preMove(); |
| 311 |
> |
//save old state (position, velocity etc) |
| 312 |
> |
rattle->doPreConstraint(); |
| 313 |
|
|
| 314 |
|
#ifdef PROFILE |
| 315 |
|
endProfile(pro3); |
| 405 |
|
} |
| 406 |
|
} |
| 407 |
|
|
| 408 |
< |
if (nConstrained){ |
| 396 |
< |
constrainA(); |
| 397 |
< |
} |
| 408 |
> |
rattle->doRattleA(); |
| 409 |
|
} |
| 410 |
|
|
| 411 |
|
|
| 446 |
|
} |
| 447 |
|
} |
| 448 |
|
|
| 449 |
< |
if (nConstrained){ |
| 439 |
< |
constrainB(); |
| 440 |
< |
} |
| 449 |
> |
rattle->doRattleB(); |
| 450 |
|
} |
| 451 |
|
|
| 452 |
+ |
/* |
| 453 |
|
template<typename T> void Integrator<T>::preMove(void){ |
| 454 |
|
int i, j; |
| 455 |
|
double pos[3]; |
| 708 |
|
simError(); |
| 709 |
|
} |
| 710 |
|
} |
| 711 |
< |
|
| 711 |
> |
*/ |
| 712 |
|
template<typename T> void Integrator<T>::rotationPropagation |
| 713 |
|
( StuntDouble* sd, double ji[3] ){ |
| 714 |
|
|
