| 57 |
|
#include "primitives/Bend.hpp" |
| 58 |
|
#include "primitives/Torsion.hpp" |
| 59 |
|
#include "primitives/Inversion.hpp" |
| 60 |
+ |
#include "parallel/ForceDecomposition.hpp" |
| 61 |
+ |
//#include "parallel/SerialDecomposition.hpp" |
| 62 |
|
|
| 63 |
|
namespace OpenMD { |
| 64 |
|
|
| 65 |
|
ForceManager::ForceManager(SimInfo * info) : info_(info), |
| 66 |
|
NBforcesInitialized_(false) { |
| 67 |
+ |
#ifdef IS_MPI |
| 68 |
+ |
decomp_ = new ForceDecomposition(info_); |
| 69 |
+ |
#else |
| 70 |
+ |
// decomp_ = new SerialDecomposition(info); |
| 71 |
+ |
#endif |
| 72 |
|
} |
| 73 |
|
|
| 74 |
|
void ForceManager::calcForces() { |
| 77 |
|
if (!info_->isFortranInitialized()) { |
| 78 |
|
info_->update(); |
| 79 |
|
nbiMan_->setSimInfo(info_); |
| 80 |
< |
nbiMan_->initialize(); |
| 80 |
> |
nbiMan_->initialize(); |
| 81 |
> |
decomp_->distributeInitialData(); |
| 82 |
|
info_->setupFortran(); |
| 83 |
|
} |
| 84 |
|
|
| 85 |
< |
preCalculation(); |
| 78 |
< |
|
| 85 |
> |
preCalculation(); |
| 86 |
|
calcShortRangeInteraction(); |
| 80 |
– |
|
| 87 |
|
calcLongRangeInteraction(); |
| 82 |
– |
|
| 88 |
|
postCalculation(); |
| 89 |
|
|
| 90 |
|
} |
| 96 |
|
Atom* atom; |
| 97 |
|
Molecule::RigidBodyIterator rbIter; |
| 98 |
|
RigidBody* rb; |
| 99 |
+ |
Molecule::CutoffGroupIterator ci; |
| 100 |
+ |
CutoffGroup* cg; |
| 101 |
|
|
| 102 |
|
// forces are zeroed here, before any are accumulated. |
| 103 |
|
// NOTE: do not rezero the forces in Fortran. |
| 113 |
|
rb = mol->nextRigidBody(rbIter)) { |
| 114 |
|
rb->zeroForcesAndTorques(); |
| 115 |
|
} |
| 116 |
< |
|
| 116 |
> |
|
| 117 |
> |
if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){ |
| 118 |
> |
std::cerr << "should not see me \n"; |
| 119 |
> |
for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
| 120 |
> |
cg = mol->nextCutoffGroup(ci)) { |
| 121 |
> |
//calculate the center of mass of cutoff group |
| 122 |
> |
cg->updateCOM(); |
| 123 |
> |
} |
| 124 |
> |
} |
| 125 |
|
} |
| 126 |
< |
|
| 126 |
> |
|
| 127 |
|
// Zero out the stress tensor |
| 128 |
|
tau *= 0.0; |
| 129 |
|
|
| 250 |
|
void ForceManager::calcLongRangeInteraction() { |
| 251 |
|
Snapshot* curSnapshot; |
| 252 |
|
DataStorage* config; |
| 253 |
+ |
DataStorage* cgConfig; |
| 254 |
|
RealType* frc; |
| 255 |
|
RealType* pos; |
| 256 |
|
RealType* trq; |
| 264 |
|
|
| 265 |
|
//get array pointers |
| 266 |
|
config = &(curSnapshot->atomData); |
| 267 |
+ |
cgConfig = &(curSnapshot->cgData); |
| 268 |
|
frc = config->getArrayPointer(DataStorage::dslForce); |
| 269 |
|
pos = config->getArrayPointer(DataStorage::dslPosition); |
| 270 |
|
trq = config->getArrayPointer(DataStorage::dslTorque); |
| 272 |
|
electroFrame = config->getArrayPointer(DataStorage::dslElectroFrame); |
| 273 |
|
particlePot = config->getArrayPointer(DataStorage::dslParticlePot); |
| 274 |
|
|
| 275 |
< |
//calculate the center of mass of cutoff group |
| 276 |
< |
SimInfo::MoleculeIterator mi; |
| 277 |
< |
Molecule* mol; |
| 261 |
< |
Molecule::CutoffGroupIterator ci; |
| 262 |
< |
CutoffGroup* cg; |
| 263 |
< |
Vector3d com; |
| 264 |
< |
std::vector<Vector3d> rcGroup; |
| 265 |
< |
|
| 266 |
< |
if(info_->getNCutoffGroups() > 0){ |
| 267 |
< |
|
| 268 |
< |
for (mol = info_->beginMolecule(mi); mol != NULL; |
| 269 |
< |
mol = info_->nextMolecule(mi)) { |
| 270 |
< |
for(cg = mol->beginCutoffGroup(ci); cg != NULL; |
| 271 |
< |
cg = mol->nextCutoffGroup(ci)) { |
| 272 |
< |
cg->getCOM(com); |
| 273 |
< |
rcGroup.push_back(com); |
| 274 |
< |
} |
| 275 |
< |
}// end for (mol) |
| 276 |
< |
|
| 277 |
< |
rc = rcGroup[0].getArrayPointer(); |
| 275 |
> |
if(info_->getNGlobalCutoffGroups() != info_->getNGlobalAtoms()){ |
| 276 |
> |
std::cerr << "should not see me \n"; |
| 277 |
> |
rc = cgConfig->getArrayPointer(DataStorage::dslPosition); |
| 278 |
|
} else { |
| 279 |
|
// center of mass of the group is the same as position of the atom |
| 280 |
|
// if cutoff group does not exist |
| 290 |
|
longRangePotential[i]=0.0; //Initialize array |
| 291 |
|
} |
| 292 |
|
|
| 293 |
< |
doForceLoop(pos, |
| 294 |
< |
rc, |
| 295 |
< |
A, |
| 296 |
< |
electroFrame, |
| 297 |
< |
frc, |
| 298 |
< |
trq, |
| 299 |
< |
tau.getArrayPointer(), |
| 300 |
< |
longRangePotential, |
| 301 |
< |
particlePot, |
| 302 |
< |
&isError ); |
| 303 |
< |
|
| 293 |
> |
decomp_->distributeData(); |
| 294 |
> |
|
| 295 |
> |
int nLoops = 1; |
| 296 |
> |
for (int iLoop = 0; iLoop < nLoops; iLoop++) { |
| 297 |
> |
doForceLoop(pos, |
| 298 |
> |
rc, |
| 299 |
> |
A, |
| 300 |
> |
electroFrame, |
| 301 |
> |
frc, |
| 302 |
> |
trq, |
| 303 |
> |
tau.getArrayPointer(), |
| 304 |
> |
longRangePotential, |
| 305 |
> |
particlePot, |
| 306 |
> |
&isError ); |
| 307 |
> |
|
| 308 |
> |
if (nLoops > 1) { |
| 309 |
> |
decomp_->collectIntermediateData(); |
| 310 |
> |
decomp_->distributeIntermediateData(); |
| 311 |
> |
} |
| 312 |
> |
} |
| 313 |
> |
|
| 314 |
> |
decomp_->collectData(); |
| 315 |
> |
|
| 316 |
|
if( isError ){ |
| 317 |
|
sprintf( painCave.errMsg, |
| 318 |
|
"Error returned from the fortran force calculation.\n" ); |
