| 39 |
|
* such damages. |
| 40 |
|
*/ |
| 41 |
|
#include <fstream> |
| 42 |
+ |
#include <iostream> |
| 43 |
|
#include "integrators/LDForceManager.hpp" |
| 44 |
|
#include "math/CholeskyDecomposition.hpp" |
| 45 |
|
#include "utils/OOPSEConstant.hpp" |
| 51 |
|
namespace oopse { |
| 52 |
|
|
| 53 |
|
LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){ |
| 54 |
< |
Globals* simParams = info->getSimParams(); |
| 55 |
< |
|
| 54 |
> |
simParams = info->getSimParams(); |
| 55 |
> |
veloMunge = new Velocitizer(info); |
| 56 |
> |
|
| 57 |
|
sphericalBoundaryConditions_ = false; |
| 58 |
|
if (simParams->getUseSphericalBoundaryConditions()) { |
| 59 |
|
sphericalBoundaryConditions_ = true; |
| 247 |
|
return props; |
| 248 |
|
} |
| 249 |
|
|
| 250 |
< |
void LDForceManager::postCalculation() { |
| 250 |
> |
void LDForceManager::postCalculation(bool needStress){ |
| 251 |
|
SimInfo::MoleculeIterator i; |
| 252 |
|
Molecule::IntegrableObjectIterator j; |
| 253 |
|
Molecule* mol; |
| 259 |
|
Mat3x3d Atrans; |
| 260 |
|
Vector3d Tb; |
| 261 |
|
Vector3d ji; |
| 260 |
– |
RealType mass; |
| 262 |
|
unsigned int index = 0; |
| 263 |
|
bool doLangevinForces; |
| 264 |
|
bool freezeMolecule; |
| 265 |
|
int fdf; |
| 266 |
< |
|
| 266 |
> |
|
| 267 |
|
fdf = 0; |
| 268 |
+ |
|
| 269 |
|
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
| 270 |
|
|
| 271 |
|
doLangevinForces = true; |
| 275 |
|
|
| 276 |
|
Vector3d molPos = mol->getCom(); |
| 277 |
|
RealType molRad = molPos.length(); |
| 278 |
< |
|
| 278 |
> |
|
| 279 |
|
doLangevinForces = false; |
| 280 |
|
|
| 281 |
|
if (molRad > langevinBufferRadius_) { |
| 294 |
|
if (freezeMolecule) |
| 295 |
|
fdf += integrableObject->freeze(); |
| 296 |
|
|
| 297 |
< |
if (doLangevinForces) { |
| 297 |
> |
if (doLangevinForces) { |
| 298 |
|
vel =integrableObject->getVel(); |
| 299 |
|
if (integrableObject->isDirectional()){ |
| 300 |
|
//calculate angular velocity in lab frame |
| 352 |
|
|
| 353 |
|
} |
| 354 |
|
} |
| 355 |
+ |
|
| 356 |
|
info_->setFdf(fdf); |
| 357 |
< |
|
| 358 |
< |
ForceManager::postCalculation(); |
| 357 |
> |
veloMunge->removeComDrift(); |
| 358 |
> |
// Remove angular drift if we are not using periodic boundary conditions. |
| 359 |
> |
if(!simParams->getUsePeriodicBoundaryConditions()) |
| 360 |
> |
veloMunge->removeAngularDrift(); |
| 361 |
> |
|
| 362 |
> |
ForceManager::postCalculation(needStress); |
| 363 |
|
} |
| 364 |
|
|
| 365 |
|
void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) { |
| 367 |
|
|
| 368 |
|
Vector<RealType, 6> Z; |
| 369 |
|
Vector<RealType, 6> generalForce; |
| 363 |
– |
|
| 370 |
|
|
| 371 |
|
Z[0] = randNumGen_.randNorm(0, variance); |
| 372 |
|
Z[1] = randNumGen_.randNorm(0, variance); |
