--- trunk/src/brains/ForceManager.cpp 2005/11/16 23:10:02 749 +++ trunk/src/brains/ForceManager.cpp 2006/05/17 21:51:42 963 @@ -57,6 +57,7 @@ namespace oopse { #include "primitives/Bend.hpp" namespace oopse { +/* struct BendOrderStruct { Bend* bend; BendDataSet dataSet; @@ -73,7 +74,7 @@ namespace oopse { 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()) { @@ -88,6 +89,7 @@ namespace oopse { postCalculation(); +/* std::vector bendOrderStruct; for(std::map::iterator i = bendDataSets.begin(); i != bendDataSets.end(); ++i) { BendOrderStruct tmp; @@ -106,19 +108,17 @@ namespace oopse { 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 << "Bend: 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 << "Torsion: 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)) { @@ -172,12 +172,12 @@ namespace oopse { bondPotential += bond->getPotential(); } - //int i =0; + for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { - //std::cout << i++ << "\t"; - double angle; + + RealType angle; bend->calcForce(angle); - double currBendPot = bend->getPotential(); + RealType currBendPot = bend->getPotential(); bendPotential += bend->getPotential(); std::map::iterator i = bendDataSets.find(bend); if (i == bendDataSets.end()) { @@ -196,9 +196,9 @@ namespace oopse { } for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { - double angle; + RealType angle; torsion->calcForce(angle); - double currTorsionPot = torsion->getPotential(); + RealType currTorsionPot = torsion->getPotential(); torsionPotential += torsion->getPotential(); std::map::iterator i = torsionDataSets.find(torsion); if (i == torsionDataSets.end()) { @@ -218,7 +218,7 @@ namespace oopse { } - double shortRangePotential = bondPotential + bendPotential + torsionPotential; + RealType shortRangePotential = bondPotential + bendPotential + torsionPotential; Snapshot* curSnapshot = info_->getSnapshotManager()->getCurrentSnapshot(); curSnapshot->statData[Stats::SHORT_RANGE_POTENTIAL] = shortRangePotential; curSnapshot->statData[Stats::BOND_POTENTIAL] = bondPotential; @@ -230,12 +230,12 @@ namespace oopse { 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(); @@ -272,8 +272,8 @@ namespace oopse { } //initialize data before passing to fortran - double longRangePotential[LR_POT_TYPES]; - double lrPot = 0.0; + RealType longRangePotential[LR_POT_TYPES]; + RealType lrPot = 0.0; Mat3x3d tau; short int passedCalcPot = needPotential;