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; |
90 |
|
} |
91 |
|
|
92 |
|
// Build the hydroProp map: |
93 |
< |
std::map<std::string, HydroProp> hydroPropMap; |
93 |
> |
std::map<std::string, HydroProp*> hydroPropMap; |
94 |
|
|
95 |
|
Molecule* mol; |
96 |
|
StuntDouble* integrableObject; |
132 |
|
integrableObject != NULL; |
133 |
|
integrableObject = mol->nextIntegrableObject(j)) { |
134 |
|
|
135 |
< |
std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType()); |
135 |
> |
std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType()); |
136 |
|
if (iter != hydroPropMap.end()) { |
137 |
|
hydroProps_.push_back(iter->second); |
138 |
|
} else { |
145 |
|
} |
146 |
|
} |
147 |
|
} else { |
148 |
< |
|
149 |
< |
std::map<std::string, HydroProp> hydroPropMap; |
148 |
> |
|
149 |
> |
std::map<std::string, HydroProp*> hydroPropMap; |
150 |
|
for (mol = info->beginMolecule(i); mol != NULL; |
151 |
|
mol = info->nextMolecule(i)) { |
152 |
|
for (integrableObject = mol->beginIntegrableObject(j); |
166 |
|
if (gayBerneData != NULL) { |
167 |
|
GayBerneParam gayBerneParam = gayBerneData->getData(); |
168 |
|
currShape = new Ellipsoid(V3Zero, |
169 |
< |
gayBerneParam.GB_sigma/2.0, |
170 |
< |
gayBerneParam.GB_l2b_ratio*gayBerneParam.GB_sigma/2.0, |
169 |
> |
gayBerneParam.GB_d / 2.0, |
170 |
> |
gayBerneParam.GB_l / 2.0, |
171 |
|
Mat3x3d::identity()); |
172 |
|
} else { |
173 |
|
sprintf( painCave.errMsg, |
215 |
|
} |
216 |
|
} |
217 |
|
} |
218 |
< |
HydroProps currHydroProp = currShape->getHydroProps(simParams->getViscosity(),simParams->getTargetTemp()); |
219 |
< |
std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType()); |
218 |
> |
HydroProp* currHydroProp = currShape->getHydroProp(simParams->getViscosity(),simParams->getTargetTemp()); |
219 |
> |
std::map<std::string, HydroProp*>::iterator iter = hydroPropMap.find(integrableObject->getType()); |
220 |
|
if (iter != hydroPropMap.end()) |
221 |
|
hydroProps_.push_back(iter->second); |
222 |
|
else { |
223 |
< |
HydroProp myProp; |
224 |
< |
myProp.cor = V3Zero; |
225 |
< |
for (int i1 = 0; i1 < 3; i1++) { |
224 |
< |
for (int j1 = 0; j1 < 3; j1++) { |
225 |
< |
myProp.Xirtt(i1,j1) = currHydroProp.Xi(i1,j1); |
226 |
< |
myProp.Xirrt(i1,j1) = currHydroProp.Xi(i1,j1+3); |
227 |
< |
myProp.Xirtr(i1,j1) = currHydroProp.Xi(i1+3,j1); |
228 |
< |
myProp.Xirrr(i1,j1) = currHydroProp.Xi(i1+3,j1+3); |
229 |
< |
} |
230 |
< |
} |
231 |
< |
CholeskyDecomposition(currHydroProp.Xi, myProp.S); |
232 |
< |
hydroPropMap.insert(std::map<std::string, HydroProp>::value_type(integrableObject->getType(), myProp)); |
233 |
< |
hydroProps_.push_back(myProp); |
223 |
> |
currHydroProp->complete(); |
224 |
> |
hydroPropMap.insert(std::map<std::string, HydroProp*>::value_type(integrableObject->getType(), currHydroProp)); |
225 |
> |
hydroProps_.push_back(currHydroProp); |
226 |
|
} |
227 |
|
} |
228 |
|
} |
229 |
|
} |
230 |
|
variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt(); |
231 |
< |
} |
240 |
< |
|
241 |
< |
|
231 |
> |
} |
232 |
|
|
233 |
< |
|
234 |
< |
|
245 |
< |
std::map<std::string, HydroProp> LDForceManager::parseFrictionFile(const std::string& filename) { |
246 |
< |
std::map<std::string, HydroProp> props; |
233 |
> |
std::map<std::string, HydroProp*> LDForceManager::parseFrictionFile(const std::string& filename) { |
234 |
> |
std::map<std::string, HydroProp*> props; |
235 |
|
std::ifstream ifs(filename.c_str()); |
236 |
|
if (ifs.is_open()) { |
237 |
|
|
240 |
|
const unsigned int BufferSize = 65535; |
241 |
|
char buffer[BufferSize]; |
242 |
|
while (ifs.getline(buffer, BufferSize)) { |
243 |
< |
StringTokenizer tokenizer(buffer); |
244 |
< |
HydroProp currProp; |
257 |
< |
if (tokenizer.countTokens() >= 40) { |
258 |
< |
std::string atomName = tokenizer.nextToken(); |
259 |
< |
currProp.cor[0] = tokenizer.nextTokenAsDouble(); |
260 |
< |
currProp.cor[1] = tokenizer.nextTokenAsDouble(); |
261 |
< |
currProp.cor[2] = tokenizer.nextTokenAsDouble(); |
262 |
< |
|
263 |
< |
currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble(); |
264 |
< |
currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble(); |
265 |
< |
currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble(); |
266 |
< |
currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble(); |
267 |
< |
currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble(); |
268 |
< |
currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble(); |
269 |
< |
currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble(); |
270 |
< |
currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble(); |
271 |
< |
currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble(); |
272 |
< |
|
273 |
< |
currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble(); |
274 |
< |
currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble(); |
275 |
< |
currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble(); |
276 |
< |
currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble(); |
277 |
< |
currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble(); |
278 |
< |
currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble(); |
279 |
< |
currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble(); |
280 |
< |
currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble(); |
281 |
< |
currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble(); |
282 |
< |
|
283 |
< |
currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble(); |
284 |
< |
currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble(); |
285 |
< |
currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble(); |
286 |
< |
currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble(); |
287 |
< |
currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble(); |
288 |
< |
currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble(); |
289 |
< |
currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble(); |
290 |
< |
currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble(); |
291 |
< |
currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble(); |
292 |
< |
|
293 |
< |
currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble(); |
294 |
< |
currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble(); |
295 |
< |
currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble(); |
296 |
< |
currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble(); |
297 |
< |
currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble(); |
298 |
< |
currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble(); |
299 |
< |
currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble(); |
300 |
< |
currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble(); |
301 |
< |
currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble(); |
302 |
< |
|
303 |
< |
SquareMatrix<RealType, 6> Xir; |
304 |
< |
Xir.setSubMatrix(0, 0, currProp.Xirtt); |
305 |
< |
Xir.setSubMatrix(0, 3, currProp.Xirrt); |
306 |
< |
Xir.setSubMatrix(3, 0, currProp.Xirtr); |
307 |
< |
Xir.setSubMatrix(3, 3, currProp.Xirrr); |
308 |
< |
CholeskyDecomposition(Xir, currProp.S); |
309 |
< |
|
310 |
< |
props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp)); |
311 |
< |
} |
243 |
> |
HydroProp* currProp = new HydroProp(buffer); |
244 |
> |
props.insert(std::map<std::string, HydroProp*>::value_type(currProp->getName(), currProp)); |
245 |
|
} |
246 |
< |
|
246 |
> |
|
247 |
|
return props; |
248 |
|
} |
249 |
< |
|
249 |
> |
|
250 |
|
void LDForceManager::postCalculation() { |
251 |
|
SimInfo::MoleculeIterator i; |
252 |
|
Molecule::IntegrableObjectIterator j; |
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 |
321 |
|
//apply friction force and torque at center of resistance |
322 |
|
A = integrableObject->getA(); |
323 |
|
Atrans = A.transpose(); |
324 |
< |
Vector3d rcr = Atrans * hydroProps_[index].cor; |
324 |
> |
Vector3d rcr = Atrans * hydroProps_[index]->getCOR(); |
325 |
|
Vector3d vcdLab = vel + cross(omega, rcr); |
326 |
|
Vector3d vcdBody = A* vcdLab; |
327 |
< |
Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega); |
327 |
> |
Vector3d frictionForceBody = -(hydroProps_[index]->getXitt() * vcdBody + hydroProps_[index]->getXirt() * omega); |
328 |
|
Vector3d frictionForceLab = Atrans*frictionForceBody; |
329 |
|
integrableObject->addFrc(frictionForceLab); |
330 |
< |
Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega); |
330 |
> |
Vector3d frictionTorqueBody = - (hydroProps_[index]->getXitr() * vcdBody + hydroProps_[index]->getXirr() * omega); |
331 |
|
Vector3d frictionTorqueLab = Atrans*frictionTorqueBody; |
332 |
|
integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab)); |
333 |
|
|
342 |
|
|
343 |
|
} else { |
344 |
|
//spherical atom |
345 |
< |
Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel); |
345 |
> |
Vector3d frictionForce = -(hydroProps_[index]->getXitt() * vel); |
346 |
|
Vector3d randomForce; |
347 |
|
Vector3d randomTorque; |
348 |
|
genRandomForceAndTorque(randomForce, randomTorque, index, variance_); |
355 |
|
|
356 |
|
} |
357 |
|
} |
358 |
+ |
|
359 |
|
info_->setFdf(fdf); |
360 |
< |
|
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(); |
366 |
|
} |
367 |
|
|
380 |
|
Z[5] = randNumGen_.randNorm(0, variance); |
381 |
|
|
382 |
|
|
383 |
< |
generalForce = hydroProps_[index].S*Z; |
383 |
> |
generalForce = hydroProps_[index]->getS()*Z; |
384 |
|
|
385 |
|
force[0] = generalForce[0]; |
386 |
|
force[1] = generalForce[1]; |