--- trunk/src/brains/ForceManager.cpp 2005/01/12 22:41:40 246 +++ trunk/src/brains/ForceManager.cpp 2005/10/17 23:13:44 681 @@ -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,26 @@ * 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" namespace oopse { -void ForceManager::calcForces(bool needPotential, bool needStress) { + void ForceManager::calcForces(bool needPotential, bool needStress) { if (!info_->isFortranInitialized()) { - info_->update(); + info_->update(); } preCalculation(); @@ -67,9 +69,9 @@ void ForceManager::calcForces(bool needPotential, bool postCalculation(); -} + } -void ForceManager::preCalculation() { + void ForceManager::preCalculation() { SimInfo::MoleculeIterator mi; Molecule* mol; Molecule::AtomIterator ai; @@ -80,19 +82,19 @@ void ForceManager::preCalculation() { // 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(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; @@ -107,35 +109,56 @@ void ForceManager::calcShortRangeInteraction() { //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(); + } - for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { - bend->calcForce(); - } + 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 (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { + torsion->calcForce(); + } } - double shortRangePotential = 0.0; + + double bondPotential = 0.0; + double bendPotential = 0.0; + double torsionPotential = 0.0; + for (mol = info_->beginMolecule(mi); mol != NULL; mol = info_->nextMolecule(mi)) { - shortRangePotential += mol->getPotential(); - } + for (bond = mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + bondPotential += bond->getPotential(); + } + + for (bend = mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + bendPotential += bend->getPotential(); + } + + for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { + torsionPotential += torsion->getPotential(); + } + + } + + 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 +186,70 @@ 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::LONG_RANGE_POTENTIAL] = longRangePotential; + 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 +257,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