| 2 |
|
#include <stdlib.h> |
| 3 |
|
#include <math.h> |
| 4 |
|
#ifdef IS_MPI |
| 5 |
< |
#include "mpiSimulation.hpp" |
| 5 |
> |
#include "brains/mpiSimulation.hpp" |
| 6 |
|
#include <unistd.h> |
| 7 |
|
#endif //is_mpi |
| 8 |
|
|
| 9 |
|
#ifdef PROFILE |
| 10 |
< |
#include "mdProfile.hpp" |
| 10 |
> |
#include "profiling/mdProfile.hpp" |
| 11 |
|
#endif // profile |
| 12 |
|
|
| 13 |
< |
#include "Integrator.hpp" |
| 14 |
< |
#include "simError.h" |
| 13 |
> |
#include "integrators/Integrator.hpp" |
| 14 |
> |
#include "utils/simError.h" |
| 15 |
|
|
| 16 |
|
|
| 17 |
|
template<typename T> Integrator<T>::Integrator(SimInfo* theInfo, |
| 171 |
|
double currReset; |
| 172 |
|
|
| 173 |
|
int calcPot, calcStress; |
| 174 |
+ |
int i; |
| 175 |
+ |
int localIndex; |
| 176 |
|
|
| 177 |
+ |
#ifdef IS_MPI |
| 178 |
+ |
int which_node; |
| 179 |
+ |
#endif // is_mpi |
| 180 |
+ |
|
| 181 |
+ |
vector<StuntDouble*> particles; |
| 182 |
+ |
string inAngle; |
| 183 |
+ |
|
| 184 |
|
tStats = new Thermo(info); |
| 185 |
|
statOut = new StatWriter(info); |
| 186 |
|
dumpOut = new DumpWriter(info); |
| 187 |
|
|
| 188 |
+ |
if (info->useSolidThermInt && !info->useLiquidThermInt) { |
| 189 |
+ |
restOut = new RestraintWriter(info); |
| 190 |
+ |
initRestraints = new RestraintReader(info); |
| 191 |
+ |
} |
| 192 |
+ |
|
| 193 |
|
atoms = info->atoms; |
| 194 |
|
|
| 195 |
|
dt = info->dt; |
| 203 |
|
|
| 204 |
|
// initialize the retraints if necessary |
| 205 |
|
if (info->useSolidThermInt && !info->useLiquidThermInt) { |
| 206 |
< |
myFF->initRestraints(); |
| 206 |
> |
initRestraints->zeroZangle(); |
| 207 |
> |
inAngle = info->zAngleName + "_in"; |
| 208 |
> |
initRestraints->readZangle(inAngle.c_str()); |
| 209 |
> |
initRestraints->readIdealCrystal(); |
| 210 |
|
} |
| 211 |
|
|
| 212 |
|
// initialize the forces before the first step |
| 196 |
– |
|
| 213 |
|
calcForce(1, 1); |
| 214 |
|
|
| 215 |
|
//execute constraint algorithm to make sure at the very beginning the system is constrained |
| 233 |
|
|
| 234 |
|
dumpOut->writeDump(info->getTime()); |
| 235 |
|
statOut->writeStat(info->getTime()); |
| 236 |
+ |
restOut->writeZangle(info->getTime()); |
| 237 |
|
|
| 221 |
– |
|
| 238 |
|
#ifdef IS_MPI |
| 239 |
|
strcpy(checkPointMsg, "The integrator is ready to go."); |
| 240 |
|
MPIcheckPoint(); |
| 270 |
|
|
| 271 |
|
if (info->getTime() >= currSample){ |
| 272 |
|
dumpOut->writeDump(info->getTime()); |
| 273 |
+ |
// write a zAng file to coincide with each dump or eor file |
| 274 |
+ |
if (info->useSolidThermInt && !info->useLiquidThermInt) |
| 275 |
+ |
restOut->writeZangle(info->getTime()); |
| 276 |
|
currSample += sampleTime; |
| 277 |
|
} |
| 278 |
|
|
| 302 |
|
|
| 303 |
|
dumpOut->writeFinal(info->getTime()); |
| 304 |
|
|
| 305 |
< |
// dump out a file containing the omega values for the final configuration |
| 306 |
< |
if (info->useSolidThermInt && !info->useLiquidThermInt) |
| 307 |
< |
myFF->dumpzAngle(); |
| 308 |
< |
|
| 305 |
> |
// write the file containing the omega values of the final configuration |
| 306 |
> |
if (info->useSolidThermInt && !info->useLiquidThermInt){ |
| 307 |
> |
restOut->writeZangle(info->getTime()); |
| 308 |
> |
restOut->writeZangle(info->getTime(), inAngle.c_str()); |
| 309 |
> |
} |
| 310 |
|
|
| 311 |
|
delete dumpOut; |
| 312 |
|
delete statOut; |
| 383 |
|
integrableObjects[i]->getVel(vel); |
| 384 |
|
integrableObjects[i]->getPos(pos); |
| 385 |
|
integrableObjects[i]->getFrc(frc); |
| 386 |
+ |
// std::cerr << "f = " << frc[0] << "\t" << frc[1] << "\t" << frc[2] << "\n"; |
| 387 |
|
|
| 388 |
|
mass = integrableObjects[i]->getMass(); |
| 389 |
|
|
| 397 |
|
integrableObjects[i]->setVel(vel); |
| 398 |
|
integrableObjects[i]->setPos(pos); |
| 399 |
|
|
| 400 |
+ |
|
| 401 |
|
if (integrableObjects[i]->isDirectional()){ |
| 402 |
|
|
| 403 |
|
// get and convert the torque to body frame |
| 404 |
|
|
| 405 |
|
integrableObjects[i]->getTrq(Tb); |
| 406 |
+ |
|
| 407 |
+ |
// std::cerr << "t = " << Tb[0] << "\t" << Tb[1] << "\t" << Tb[2] << "\n"; |
| 408 |
|
integrableObjects[i]->lab2Body(Tb); |
| 409 |
|
|
| 410 |
|
// get the angular momentum, and propagate a half step |
| 740 |
|
sd->getI(I); |
| 741 |
|
|
| 742 |
|
if (sd->isLinear()) { |
| 743 |
+ |
|
| 744 |
|
i = sd->linearAxis(); |
| 745 |
|
j = (i+1)%3; |
| 746 |
|
k = (i+2)%3; |
| 747 |
< |
|
| 747 |
> |
|
| 748 |
|
angle = dt2 * ji[j] / I[j][j]; |
| 749 |
|
this->rotate( k, i, angle, ji, A ); |
| 750 |
|
|
