--- trunk/src/brains/ForceManager.cpp 2005/01/12 22:41:40 246 +++ trunk/src/brains/ForceManager.cpp 2005/12/02 15:38:03 770 @@ -1,4 +1,4 @@ - /* +/* * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. * * The University of Notre Dame grants you ("Licensee") a @@ -39,24 +39,46 @@ * such damages. */ - /** - * @file ForceManager.cpp - * @author tlin - * @date 11/09/2004 - * @time 10:39am - * @version 1.0 - */ +/** + * @file ForceManager.cpp + * @author tlin + * @date 11/09/2004 + * @time 10:39am + * @version 1.0 + */ #include "brains/ForceManager.hpp" #include "primitives/Molecule.hpp" #include "UseTheForce/doForces_interface.h" +#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) { +/* + struct BendOrderStruct { + Bend* bend; + BendDataSet dataSet; + }; + struct TorsionOrderStruct { + Torsion* torsion; + TorsionDataSet dataSet; + }; + bool BendSortFunctor(const BendOrderStruct& b1, const BendOrderStruct& b2) { + return b1.dataSet.deltaV < b2.dataSet.deltaV; + } + + bool TorsionSortFunctor(const TorsionOrderStruct& t1, const TorsionOrderStruct& t2) { + return t1.dataSet.deltaV < t2.dataSet.deltaV; + } + */ + void ForceManager::calcForces(bool needPotential, bool needStress) { + if (!info_->isFortranInitialized()) { - info_->update(); + info_->update(); } preCalculation(); @@ -66,10 +88,42 @@ void ForceManager::calcForces(bool needPotential, bool calcLongRangeInteraction(needPotential, needStress); postCalculation(); - -} -void ForceManager::preCalculation() { +/* + std::vector bendOrderStruct; + for(std::map::iterator i = bendDataSets.begin(); i != bendDataSets.end(); ++i) { + BendOrderStruct tmp; + tmp.bend= const_cast(i->first); + tmp.dataSet = i->second; + bendOrderStruct.push_back(tmp); + } + + std::vector torsionOrderStruct; + for(std::map::iterator j = torsionDataSets.begin(); j != torsionDataSets.end(); ++j) { + TorsionOrderStruct tmp; + tmp.torsion = const_cast(j->first); + tmp.dataSet = j->second; + torsionOrderStruct.push_back(tmp); + } + + std::sort(bendOrderStruct.begin(), bendOrderStruct.end(), std::ptr_fun(BendSortFunctor)); + std::sort(torsionOrderStruct.begin(), torsionOrderStruct.end(), std::ptr_fun(TorsionSortFunctor)); + std::cout << "bend" << std::endl; + for (std::vector::iterator k = bendOrderStruct.begin(); k != bendOrderStruct.end(); ++k) { + Bend* bend = k->bend; + std::cout << "atom1=" <getAtomA()->getGlobalIndex() << ",atom2 = "<< bend->getAtomB()->getGlobalIndex() << ",atom3="<getAtomC()->getGlobalIndex() << " "; + std::cout << "deltaV=" << k->dataSet.deltaV << ",p_theta=" << k->dataSet.prev.angle <<",p_pot=" << k->dataSet.prev.potential<< ",c_theta=" << k->dataSet.curr.angle << ", c_pot = " << k->dataSet.curr.potential <::iterator l = torsionOrderStruct.begin(); l != torsionOrderStruct.end(); ++l) { + Torsion* torsion = l->torsion; + std::cout << "atom1=" <getAtomA()->getGlobalIndex() << ",atom2 = "<< torsion->getAtomB()->getGlobalIndex() << ",atom3="<getAtomC()->getGlobalIndex() << ",atom4="<getAtomD()->getGlobalIndex()<< " "; + std::cout << "deltaV=" << l->dataSet.deltaV << ",p_theta=" << l->dataSet.prev.angle <<",p_pot=" << l->dataSet.prev.potential<< ",c_theta=" << l->dataSet.curr.angle << ", c_pot = " << l->dataSet.curr.potential <beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { - atom->zeroForcesAndTorques(); - } + 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)) { - rb->zeroForcesAndTorques(); - } + //change the positions of atoms which belong to the rigidbodies + for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + rb->zeroForcesAndTorques(); + } } -} + } -void ForceManager::calcShortRangeInteraction() { + void ForceManager::calcShortRangeInteraction() { Molecule* mol; RigidBody* rb; Bond* bond; @@ -103,39 +157,79 @@ void ForceManager::calcShortRangeInteraction() { Molecule::BondIterator bondIter;; Molecule::BendIterator bendIter; Molecule::TorsionIterator torsionIter; + double bondPotential = 0.0; + double bendPotential = 0.0; + double torsionPotential = 0.0; //calculate short range interactions 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(); - } + //change the positions of atoms which belong to the rigidbodies + 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 (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { - torsion->calcForce(); - } + for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { - } - - double shortRangePotential = 0.0; - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - shortRangePotential += mol->getPotential(); - } + double angle; + bend->calcForce(angle); + double 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)) { + double angle; + torsion->calcForce(angle); + double 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 = 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; @@ -163,55 +257,67 @@ void ForceManager::calcLongRangeInteraction(bool needP CutoffGroup* cg; Vector3d com; std::vector rcGroup; - - if(info_->getNCutoffGroups() > 0){ - for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { + 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)) { - cg->getCOM(com); - rcGroup.push_back(com); + cg->getCOM(com); + rcGroup.push_back(com); } - }// end for (mol) + }// end for (mol) - rc = rcGroup[0].getArrayPointer(); + 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 - rc = pos; + // 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 = 0.0; + double longRangePotential[LR_POT_TYPES]; + double lrPot = 0.0; + Mat3x3d tau; short int passedCalcPot = needPotential; short int passedCalcStress = needStress; int isError = 0; + for (int i=0; istatData[Stats::LONG_RANGE_POTENTIAL] = longRangePotential; + curSnapshot->statData[Stats::LONG_RANGE_POTENTIAL] = lrPot; + curSnapshot->statData[Stats::VANDERWAALS_POTENTIAL] = longRangePotential[VDW_POT]; + curSnapshot->statData[Stats::ELECTROSTATIC_POTENTIAL] = longRangePotential[ELECTROSTATIC_POT]; + curSnapshot->statData.setTau(tau); -} + } -void ForceManager::postCalculation() { + void ForceManager::postCalculation() { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::RigidBodyIterator rbIter; @@ -219,11 +325,11 @@ void ForceManager::postCalculation() { // 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 (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + rb->calcForcesAndTorques(); + } } -} + } } //end namespace oopse