--- trunk/src/brains/ForceManager.cpp 2005/10/12 21:57:16 664 +++ trunk/src/brains/ForceManager.cpp 2007/04/06 21:53:43 1126 @@ -53,24 +53,26 @@ #define __C #include "UseTheForce/DarkSide/fInteractionMap.h" #include "utils/simError.h" +#include "primitives/Bend.hpp" +#include "primitives/Bend.hpp" namespace oopse { void ForceManager::calcForces(bool needPotential, bool needStress) { - + if (!info_->isFortranInitialized()) { info_->update(); } - + preCalculation(); calcShortRangeInteraction(); calcLongRangeInteraction(needPotential, needStress); - postCalculation(); - + postCalculation(needStress); + } - + void ForceManager::preCalculation() { SimInfo::MoleculeIterator mi; Molecule* mol; @@ -81,19 +83,25 @@ namespace oopse { // forces are zeroed here, before any are accumulated. // NOTE: do not rezero the forces in Fortran. - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { atom->zeroForcesAndTorques(); } - + //change the positions of atoms which belong to the rigidbodies - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { rb->zeroForcesAndTorques(); } } + // Zero out the stress tensor + tau *= 0.0; + } - + void ForceManager::calcShortRangeInteraction() { Molecule* mol; RigidBody* rb; @@ -105,51 +113,98 @@ namespace oopse { Molecule::BondIterator bondIter;; Molecule::BendIterator bendIter; Molecule::TorsionIterator torsionIter; + RealType bondPotential = 0.0; + RealType bendPotential = 0.0; + RealType torsionPotential = 0.0; //calculate short range interactions - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { //change the positions of atoms which belong to the rigidbodies - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { - rb->updateAtoms(); + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + rb->updateAtoms(); } - for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { - bond->calcForce(); + for (bond = mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { + bond->calcForce(); + bondPotential += bond->getPotential(); } - for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { - bend->calcForce(); + for (bend = mol->beginBend(bendIter); bend != NULL; + bend = mol->nextBend(bendIter)) { + + RealType angle; + bend->calcForce(angle); + RealType currBendPot = bend->getPotential(); + bendPotential += bend->getPotential(); + std::map::iterator i = bendDataSets.find(bend); + if (i == bendDataSets.end()) { + BendDataSet dataSet; + dataSet.prev.angle = dataSet.curr.angle = angle; + dataSet.prev.potential = dataSet.curr.potential = currBendPot; + dataSet.deltaV = 0.0; + bendDataSets.insert(std::map::value_type(bend, dataSet)); + }else { + i->second.prev.angle = i->second.curr.angle; + i->second.prev.potential = i->second.curr.potential; + i->second.curr.angle = angle; + i->second.curr.potential = currBendPot; + i->second.deltaV = fabs(i->second.curr.potential - + i->second.prev.potential); + } } - - for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { - torsion->calcForce(); - } - + + for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; + torsion = mol->nextTorsion(torsionIter)) { + RealType angle; + torsion->calcForce(angle); + RealType currTorsionPot = torsion->getPotential(); + torsionPotential += torsion->getPotential(); + std::map::iterator i = torsionDataSets.find(torsion); + if (i == torsionDataSets.end()) { + TorsionDataSet dataSet; + dataSet.prev.angle = dataSet.curr.angle = angle; + dataSet.prev.potential = dataSet.curr.potential = currTorsionPot; + dataSet.deltaV = 0.0; + torsionDataSets.insert(std::map::value_type(torsion, dataSet)); + }else { + i->second.prev.angle = i->second.curr.angle; + i->second.prev.potential = i->second.curr.potential; + i->second.curr.angle = angle; + i->second.curr.potential = currTorsionPot; + i->second.deltaV = fabs(i->second.curr.potential - + i->second.prev.potential); + } + } } - double shortRangePotential = 0.0; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - shortRangePotential += mol->getPotential(); - } - + RealType shortRangePotential = bondPotential + bendPotential + + torsionPotential; Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; + curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential; + curSnapshot->statData[Stats::BEND_POTENTIAL] = bendPotential; + curSnapshot->statData[Stats::DIHEDRAL_POTENTIAL] = torsionPotential; + } - - void ForceManager::calcLongRangeInteraction(bool needPotential, bool needStress) { + + void ForceManager::calcLongRangeInteraction(bool needPotential, + bool needStress) { Snapshot* curSnapshot; DataStorage* config; - double* frc; - double* pos; - double* trq; - double* A; - double* electroFrame; - double* rc; + RealType* frc; + RealType* pos; + RealType* trq; + RealType* A; + RealType* electroFrame; + RealType* rc; //get current snapshot from SimInfo curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); - + //get array pointers config = &(curSnapshot->atomData); frc = config->getArrayPointer(DataStorage::dslForce); @@ -165,11 +220,13 @@ namespace oopse { CutoffGroup* cg; Vector3d com; std::vector rcGroup; - + if(info_->getNCutoffGroups() > 0){ - - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for(cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { + + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for(cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { cg->getCOM(com); rcGroup.push_back(com); } @@ -177,15 +234,15 @@ namespace oopse { rc = rcGroup[0].getArrayPointer(); } else { - // center of mass of the group is the same as position of the atom if cutoff group does not exist + // center of mass of the group is the same as position of the atom + // if cutoff group does not exist rc = pos; } - - //initialize data before passing to fortran - double longRangePotential[LR_POT_TYPES]; - double lrPot = 0.0; - Mat3x3d tau; + //initialize data before passing to fortran + RealType longRangePotential[LR_POT_TYPES]; + RealType lrPot = 0.0; + Vector3d totalDipole; short int passedCalcPot = needPotential; short int passedCalcStress = needStress; int isError = 0; @@ -193,9 +250,7 @@ namespace oopse { for (int i=0; igetCalcBoxDipole()){ + getAccumulatedBoxDipole(totalDipole.getArrayPointer()); + + curSnapshot->statData[Stats::BOX_DIPOLE_X] = totalDipole(0); + curSnapshot->statData[Stats::BOX_DIPOLE_Y] = totalDipole(1); + curSnapshot->statData[Stats::BOX_DIPOLE_Z] = totalDipole(2); + } + //store the tau and long range potential curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; - // curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; - curSnapshot->statData.setTau(tau); + curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; + curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; } - - void ForceManager::postCalculation() { + + void ForceManager::postCalculation(bool needStress) { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::RigidBodyIterator rbIter; RigidBody* rb; + Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); // collect the atomic forces onto rigid bodies - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { - rb->calcForcesAndTorques(); + + for (mol = info_->beginMolecule(mi); mol != NULL; + mol = info_->nextMolecule(mi)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + if (needStress) { + Mat3x3d rbTau = rb->calcForcesAndTorquesAndVirial(); + tau += rbTau; + } else{ + rb->calcForcesAndTorques(); + } } } + if (needStress) { +#ifdef IS_MPI + Mat3x3d tmpTau(tau); + MPI_Allreduce(tmpTau.getArrayPointer(), tau.getArrayPointer(), + 9, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); +#endif + curSnapshot->statData.setTau(tau); + } } } //end namespace oopse