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