45 |
|
#include "utils/OOPSEConstant.hpp" |
46 |
|
#include "hydrodynamics/Sphere.hpp" |
47 |
|
#include "hydrodynamics/Ellipsoid.hpp" |
48 |
< |
#include "openbabel/mol.hpp" |
48 |
> |
#include "utils/ElementsTable.hpp" |
49 |
|
|
50 |
– |
using namespace OpenBabel; |
50 |
|
namespace oopse { |
51 |
|
|
52 |
< |
LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){ |
52 |
> |
LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info), forceTolerance_(1e-6), maxIterNum_(4) { |
53 |
|
simParams = info->getSimParams(); |
54 |
|
veloMunge = new Velocitizer(info); |
55 |
|
|
117 |
|
hydroPropMap = parseFrictionFile(simParams->getHydroPropFile()); |
118 |
|
} else { |
119 |
|
sprintf( painCave.errMsg, |
120 |
< |
"HydroPropFile must be set to a file name if Langevin\n" |
121 |
< |
"\tDynamics is specified for rigidBodies which contain more\n" |
122 |
< |
"\tthan one atom. To create a HydroPropFile, run \"Hydro\".\n"); |
120 |
> |
"HydroPropFile must be set to a file name if Langevin Dynamics\n" |
121 |
> |
"\tis specified for rigidBodies which contain more than one atom\n" |
122 |
> |
"\tTo create a HydroPropFile, run the \"Hydro\" program.\n"); |
123 |
|
painCave.severity = OOPSE_ERROR; |
124 |
|
painCave.isFatal = 1; |
125 |
|
simError(); |
198 |
|
} |
199 |
|
} |
200 |
|
} else { |
201 |
< |
int obanum = etab.GetAtomicNum((atom->getType()).c_str()); |
202 |
< |
if (obanum != 0) { |
203 |
< |
currShape = new Sphere(atom->getPos(), etab.GetVdwRad(obanum)); |
201 |
> |
int aNum = etab.GetAtomicNum((atom->getType()).c_str()); |
202 |
> |
if (aNum != 0) { |
203 |
> |
currShape = new Sphere(atom->getPos(), etab.GetVdwRad(aNum)); |
204 |
|
} else { |
205 |
|
sprintf( painCave.errMsg, |
206 |
|
"Could not find atom type in default element.txt\n"); |
211 |
|
} |
212 |
|
} |
213 |
|
} |
214 |
+ |
|
215 |
+ |
if (!simParams->haveTargetTemp()) { |
216 |
+ |
sprintf(painCave.errMsg, "You can't use LangevinDynamics without a targetTemp!\n"); |
217 |
+ |
painCave.isFatal = 1; |
218 |
+ |
painCave.severity = OOPSE_ERROR; |
219 |
+ |
simError(); |
220 |
+ |
} |
221 |
+ |
|
222 |
+ |
if (!simParams->haveViscosity()) { |
223 |
+ |
sprintf(painCave.errMsg, "You can't use LangevinDynamics without a viscosity!\n"); |
224 |
+ |
painCave.isFatal = 1; |
225 |
+ |
painCave.severity = OOPSE_ERROR; |
226 |
+ |
simError(); |
227 |
+ |
} |
228 |
+ |
|
229 |
+ |
|
230 |
|
HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp()); |
231 |
|
std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType()); |
232 |
|
if (iter != hydroPropMap.end()) |
265 |
|
Molecule* mol; |
266 |
|
StuntDouble* integrableObject; |
267 |
|
RealType mass; |
253 |
– |
Vector3d vel; |
268 |
|
Vector3d pos; |
269 |
|
Vector3d frc; |
270 |
|
Mat3x3d A; |
276 |
|
bool freezeMolecule; |
277 |
|
int fdf; |
278 |
|
|
265 |
– |
|
266 |
– |
|
279 |
|
fdf = 0; |
280 |
|
|
281 |
|
for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { |
307 |
|
fdf += integrableObject->freeze(); |
308 |
|
|
309 |
|
if (doLangevinForces) { |
298 |
– |
vel =integrableObject->getVel(); |
310 |
|
mass = integrableObject->getMass(); |
311 |
|
if (integrableObject->isDirectional()){ |
301 |
– |
//calculate angular velocity in lab frame |
302 |
– |
Mat3x3d I = integrableObject->getI(); |
303 |
– |
Vector3d angMom = integrableObject->getJ(); |
304 |
– |
Vector3d omega; |
305 |
– |
|
306 |
– |
if (integrableObject->isLinear()) { |
307 |
– |
int linearAxis = integrableObject->linearAxis(); |
308 |
– |
int l = (linearAxis +1 )%3; |
309 |
– |
int m = (linearAxis +2 )%3; |
310 |
– |
omega[l] = angMom[l] /I(l, l); |
311 |
– |
omega[m] = angMom[m] /I(m, m); |
312 |
– |
|
313 |
– |
} else { |
314 |
– |
omega[0] = angMom[0] /I(0, 0); |
315 |
– |
omega[1] = angMom[1] /I(1, 1); |
316 |
– |
omega[2] = angMom[2] /I(2, 2); |
317 |
– |
} |
312 |
|
|
313 |
< |
//std::cerr << "I = " << I(0,0) << "\t" << I(1,1) << "\t" << I(2,2) << "\n\n"; |
313 |
> |
// preliminaries for directional objects: |
314 |
|
|
321 |
– |
//apply friction force and torque at center of resistance |
315 |
|
A = integrableObject->getA(); |
316 |
|
Atrans = A.transpose(); |
317 |
< |
//std::cerr << "A = " << integrableObject->getA() << "\n"; |
325 |
< |
//std::cerr << "Atrans = " << A.transpose() << "\n\n"; |
326 |
< |
Vector3d rcr = Atrans * hydroProps_[index]->getCOR(); |
327 |
< |
//std::cerr << "cor = " << hydroProps_[index]->getCOR() << "\n\n\n\n"; |
328 |
< |
//std::cerr << "rcr = " << rcr << "\n\n"; |
329 |
< |
Vector3d vcdLab = vel + cross(omega, rcr); |
330 |
< |
|
331 |
< |
//std::cerr << "velL = " << vel << "\n\n"; |
332 |
< |
//std::cerr << "vcdL = " << vcdLab << "\n\n"; |
333 |
< |
Vector3d vcdBody = A* vcdLab; |
334 |
< |
//std::cerr << "vcdB = " << vcdBody << "\n\n"; |
335 |
< |
Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega); |
317 |
> |
Vector3d rcrLab = Atrans * hydroProps_[index]->getCOR(); |
318 |
|
|
337 |
– |
//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n"; |
338 |
– |
//std::cerr << "ffB = " << frictionForceBody << "\n\n"; |
339 |
– |
Vector3d frictionForceLab = Atrans*frictionForceBody; |
340 |
– |
//std::cerr << "ffL = " << frictionForceLab << "\n\n"; |
341 |
– |
//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n"; |
342 |
– |
integrableObject->addFrc(frictionForceLab); |
343 |
– |
//std::cerr << "frc = " << integrableObject->getFrc() << "\n\n"; |
344 |
– |
//std::cerr << "ome = " << omega << "\n\n"; |
345 |
– |
Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega); |
346 |
– |
//std::cerr << "ftB = " << frictionTorqueBody << "\n\n"; |
347 |
– |
Vector3d frictionTorqueLab = Atrans*frictionTorqueBody; |
348 |
– |
//std::cerr << "ftL = " << frictionTorqueLab << "\n\n"; |
349 |
– |
//std::cerr << "ftL2 = " << frictionTorqueLab+cross(rcr,frictionForceLab) << "\n\n"; |
350 |
– |
//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n"; |
351 |
– |
integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab)); |
352 |
– |
//std::cerr << "trq = " << integrableObject->getTrq() << "\n\n"; |
353 |
– |
|
319 |
|
//apply random force and torque at center of resistance |
320 |
+ |
|
321 |
|
Vector3d randomForceBody; |
322 |
|
Vector3d randomTorqueBody; |
323 |
|
genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_); |
324 |
< |
//std::cerr << "rfB = " << randomForceBody << "\n\n"; |
325 |
< |
//std::cerr << "rtB = " << randomTorqueBody << "\n\n"; |
360 |
< |
Vector3d randomForceLab = Atrans*randomForceBody; |
361 |
< |
Vector3d randomTorqueLab = Atrans* randomTorqueBody; |
324 |
> |
Vector3d randomForceLab = Atrans * randomForceBody; |
325 |
> |
Vector3d randomTorqueLab = Atrans * randomTorqueBody; |
326 |
|
integrableObject->addFrc(randomForceLab); |
327 |
< |
//std::cerr << "rfL = " << randomForceLab << "\n\n"; |
328 |
< |
//std::cerr << "rtL = " << randomTorqueLab << "\n\n"; |
329 |
< |
//std::cerr << "rtL2 = " << randomTorqueLab + cross(rcr, randomForceLab) << "\n\n"; |
330 |
< |
integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab )); |
327 |
> |
integrableObject->addTrq(randomTorqueLab + cross(rcrLab, randomForceLab )); |
328 |
> |
|
329 |
> |
Mat3x3d I = integrableObject->getI(); |
330 |
> |
Vector3d omegaBody; |
331 |
> |
|
332 |
> |
// What remains contains velocity explicitly, but the velocity required |
333 |
> |
// is at the full step: v(t + h), while we have initially the velocity |
334 |
> |
// at the half step: v(t + h/2). We need to iterate to converge the |
335 |
> |
// friction force and friction torque vectors. |
336 |
> |
|
337 |
> |
// this is the velocity at the half-step: |
338 |
> |
|
339 |
> |
Vector3d vel =integrableObject->getVel(); |
340 |
> |
Vector3d angMom = integrableObject->getJ(); |
341 |
> |
|
342 |
> |
//estimate velocity at full-step using everything but friction forces: |
343 |
> |
|
344 |
> |
frc = integrableObject->getFrc(); |
345 |
> |
Vector3d velStep = vel + (dt2_ /mass * OOPSEConstant::energyConvert) * frc; |
346 |
> |
|
347 |
> |
Tb = integrableObject->lab2Body(integrableObject->getTrq()); |
348 |
> |
Vector3d angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * Tb; |
349 |
> |
|
350 |
> |
Vector3d omegaLab; |
351 |
> |
Vector3d vcdLab; |
352 |
> |
Vector3d vcdBody; |
353 |
> |
Vector3d frictionForceBody; |
354 |
> |
Vector3d frictionForceLab(0.0); |
355 |
> |
Vector3d oldFFL; // used to test for convergence |
356 |
> |
Vector3d frictionTorqueBody(0.0); |
357 |
> |
Vector3d oldFTB; // used to test for convergence |
358 |
> |
Vector3d frictionTorqueLab; |
359 |
> |
RealType fdot; |
360 |
> |
RealType tdot; |
361 |
> |
|
362 |
> |
//iteration starts here: |
363 |
> |
|
364 |
> |
for (int k = 0; k < maxIterNum_; k++) { |
365 |
> |
|
366 |
> |
if (integrableObject->isLinear()) { |
367 |
> |
int linearAxis = integrableObject->linearAxis(); |
368 |
> |
int l = (linearAxis +1 )%3; |
369 |
> |
int m = (linearAxis +2 )%3; |
370 |
> |
omegaBody[l] = angMomStep[l] /I(l, l); |
371 |
> |
omegaBody[m] = angMomStep[m] /I(m, m); |
372 |
> |
|
373 |
> |
} else { |
374 |
> |
omegaBody[0] = angMomStep[0] /I(0, 0); |
375 |
> |
omegaBody[1] = angMomStep[1] /I(1, 1); |
376 |
> |
omegaBody[2] = angMomStep[2] /I(2, 2); |
377 |
> |
} |
378 |
> |
|
379 |
> |
omegaLab = Atrans * omegaBody; |
380 |
> |
|
381 |
> |
// apply friction force and torque at center of resistance |
382 |
> |
|
383 |
> |
vcdLab = velStep + cross(omegaLab, rcrLab); |
384 |
> |
vcdBody = A * vcdLab; |
385 |
> |
frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omegaBody); |
386 |
> |
oldFFL = frictionForceLab; |
387 |
> |
frictionForceLab = Atrans * frictionForceBody; |
388 |
> |
oldFTB = frictionTorqueBody; |
389 |
> |
frictionTorqueBody = -(hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omegaBody); |
390 |
> |
frictionTorqueLab = Atrans * frictionTorqueBody; |
391 |
> |
|
392 |
> |
// re-estimate velocities at full-step using friction forces: |
393 |
> |
|
394 |
> |
velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForceLab); |
395 |
> |
angMomStep = angMom + (dt2_ * OOPSEConstant::energyConvert) * (Tb + frictionTorqueBody); |
396 |
> |
|
397 |
> |
// check for convergence (if the vectors have converged, fdot and tdot will both be 1.0): |
398 |
> |
|
399 |
> |
fdot = dot(frictionForceLab, oldFFL) / frictionForceLab.lengthSquare(); |
400 |
> |
tdot = dot(frictionTorqueBody, oldFTB) / frictionTorqueBody.lengthSquare(); |
401 |
> |
|
402 |
> |
if (fabs(1.0 - fdot) <= forceTolerance_ && fabs(1.0 - tdot) <= forceTolerance_) |
403 |
> |
break; // iteration ends here |
404 |
> |
} |
405 |
> |
|
406 |
> |
integrableObject->addFrc(frictionForceLab); |
407 |
> |
integrableObject->addTrq(frictionTorqueLab + cross(rcrLab, frictionForceLab)); |
408 |
> |
|
409 |
|
|
410 |
|
} else { |
411 |
|
//spherical atom |
412 |
< |
Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel); |
371 |
< |
//std::cerr << "xitt = " << hydroProps_[index]->getXitt() << "\n\n"; |
412 |
> |
|
413 |
|
Vector3d randomForce; |
414 |
|
Vector3d randomTorque; |
415 |
|
genRandomForceAndTorque(randomForce, randomTorque, index, variance_); |
416 |
+ |
integrableObject->addFrc(randomForce); |
417 |
+ |
|
418 |
+ |
// What remains contains velocity explicitly, but the velocity required |
419 |
+ |
// is at the full step: v(t + h), while we have initially the velocity |
420 |
+ |
// at the half step: v(t + h/2). We need to iterate to converge the |
421 |
+ |
// friction force vector. |
422 |
+ |
|
423 |
+ |
// this is the velocity at the half-step: |
424 |
|
|
425 |
< |
integrableObject->addFrc(frictionForce+randomForce); |
425 |
> |
Vector3d vel =integrableObject->getVel(); |
426 |
> |
|
427 |
> |
//estimate velocity at full-step using everything but friction forces: |
428 |
> |
|
429 |
> |
frc = integrableObject->getFrc(); |
430 |
> |
Vector3d velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * frc; |
431 |
> |
|
432 |
> |
Vector3d frictionForce(0.0); |
433 |
> |
Vector3d oldFF; // used to test for convergence |
434 |
> |
RealType fdot; |
435 |
> |
|
436 |
> |
//iteration starts here: |
437 |
> |
|
438 |
> |
for (int k = 0; k < maxIterNum_; k++) { |
439 |
> |
|
440 |
> |
oldFF = frictionForce; |
441 |
> |
frictionForce = -hydroProps_[index]->getXitt() * velStep; |
442 |
> |
|
443 |
> |
// re-estimate velocities at full-step using friction forces: |
444 |
> |
|
445 |
> |
velStep = vel + (dt2_ / mass * OOPSEConstant::energyConvert) * (frc + frictionForce); |
446 |
> |
|
447 |
> |
// check for convergence (if the vector has converged, fdot will be 1.0): |
448 |
> |
|
449 |
> |
fdot = dot(frictionForce, oldFF) / frictionForce.lengthSquare(); |
450 |
> |
|
451 |
> |
if (fabs(1.0 - fdot) <= forceTolerance_) |
452 |
> |
break; // iteration ends here |
453 |
> |
} |
454 |
> |
|
455 |
> |
integrableObject->addFrc(frictionForce); |
456 |
> |
|
457 |
|
} |
458 |
|
} |
459 |
|
|
480 |
|
Z[0] = randNumGen_.randNorm(0, variance); |
481 |
|
Z[1] = randNumGen_.randNorm(0, variance); |
482 |
|
Z[2] = randNumGen_.randNorm(0, variance); |
403 |
– |
//Z[3] = randNumGen_.randNorm(0, variance)*(2.0*M_PI); |
404 |
– |
//Z[4] = randNumGen_.randNorm(0, variance)*(2.0*M_PI); |
405 |
– |
//Z[5] = randNumGen_.randNorm(0, variance)*(2.0*M_PI); |
483 |
|
Z[3] = randNumGen_.randNorm(0, variance); |
484 |
|
Z[4] = randNumGen_.randNorm(0, variance); |
485 |
|
Z[5] = randNumGen_.randNorm(0, variance); |
486 |
|
|
410 |
– |
|
487 |
|
generalForce = hydroProps_[index]->getS()*Z; |
488 |
|
|
489 |
|
force[0] = generalForce[0]; |