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; |
264 |
|
bool doLangevinForces; |
265 |
|
bool freezeMolecule; |
266 |
|
int fdf; |
267 |
< |
|
267 |
> |
int nIntegrated; |
268 |
> |
int nFrozen; |
269 |
> |
|
270 |
|
fdf = 0; |
271 |
+ |
|
272 |
|
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
273 |
|
|
274 |
|
doLangevinForces = true; |
278 |
|
|
279 |
|
Vector3d molPos = mol->getCom(); |
280 |
|
RealType molRad = molPos.length(); |
281 |
< |
|
281 |
> |
|
282 |
|
doLangevinForces = false; |
283 |
|
|
284 |
|
if (molRad > langevinBufferRadius_) { |
297 |
|
if (freezeMolecule) |
298 |
|
fdf += integrableObject->freeze(); |
299 |
|
|
300 |
< |
if (doLangevinForces) { |
300 |
> |
if (doLangevinForces) { |
301 |
|
vel =integrableObject->getVel(); |
302 |
|
if (integrableObject->isDirectional()){ |
303 |
|
//calculate angular velocity in lab frame |
355 |
|
|
356 |
|
} |
357 |
|
} |
358 |
+ |
|
359 |
|
info_->setFdf(fdf); |
360 |
< |
|
361 |
< |
ForceManager::postCalculation(); |
360 |
> |
veloMunge->removeComDrift(); |
361 |
> |
// Remove angular drift if we are not using periodic boundary conditions. |
362 |
> |
if(!simParams->getUsePeriodicBoundaryConditions()) |
363 |
> |
veloMunge->removeAngularDrift(); |
364 |
> |
|
365 |
> |
ForceManager::postCalculation(needStress); |
366 |
|
} |
367 |
|
|
368 |
|
void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, RealType variance) { |
370 |
|
|
371 |
|
Vector<RealType, 6> Z; |
372 |
|
Vector<RealType, 6> generalForce; |
363 |
– |
|
373 |
|
|
374 |
|
Z[0] = randNumGen_.randNorm(0, variance); |
375 |
|
Z[1] = randNumGen_.randNorm(0, variance); |