--- trunk/src/primitives/Molecule.cpp 2004/09/24 16:27:58 3 +++ trunk/src/primitives/Molecule.cpp 2006/05/17 21:51:42 963 @@ -1,256 +1,269 @@ -#include +/* + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. + * + * The University of Notre Dame grants you ("Licensee") a + * non-exclusive, royalty free, license to use, modify and + * redistribute this software in source and binary code form, provided + * that the following conditions are met: + * + * 1. Acknowledgement of the program authors must be made in any + * publication of scientific results based in part on use of the + * program. An acceptable form of acknowledgement is citation of + * the article in which the program was described (Matthew + * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher + * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented + * Parallel Simulation Engine for Molecular Dynamics," + * J. Comput. Chem. 26, pp. 252-271 (2005)) + * + * 2. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 3. Redistributions in binary form must reproduce the above copyright + * notice, this list of conditions and the following disclaimer in the + * documentation and/or other materials provided with the + * distribution. + * + * This software is provided "AS IS," without a warranty of any + * kind. All express or implied conditions, representations and + * warranties, including any implied warranty of merchantability, + * fitness for a particular purpose or non-infringement, are hereby + * excluded. The University of Notre Dame and its licensors shall not + * be liable for any damages suffered by licensee as a result of + * using, modifying or distributing the software or its + * derivatives. In no event will the University of Notre Dame or its + * licensors be liable for any lost revenue, profit or data, or for + * direct, indirect, special, consequential, incidental or punitive + * damages, however caused and regardless of the theory of liability, + * arising out of the use of or inability to use software, even if the + * University of Notre Dame has been advised of the possibility of + * such damages. + */ + +/** + * @file Molecule.cpp + * @author tlin + * @date 10/28/2004 + * @version 1.0 + */ +#include +#include #include "primitives/Molecule.hpp" +#include "utils/MemoryUtils.hpp" #include "utils/simError.h" +namespace oopse { + Molecule::Molecule(int stampId, int globalIndex, const std::string& molName) + : stampId_(stampId), globalIndex_(globalIndex), moleculeName_(molName) { + } -Molecule::Molecule( void ){ + Molecule::~Molecule() { - myAtoms = NULL; - myBonds = NULL; - myBends = NULL; - myTorsions = NULL; -} + MemoryUtils::deletePointers(atoms_); + MemoryUtils::deletePointers(bonds_); + MemoryUtils::deletePointers(bends_); + MemoryUtils::deletePointers(torsions_); + MemoryUtils::deletePointers(rigidBodies_); + MemoryUtils::deletePointers(cutoffGroups_); + MemoryUtils::deletePointers(constraintPairs_); + MemoryUtils::deletePointers(constraintElems_); + //integrableObjects_ don't own the objects + integrableObjects_.clear(); + + } -Molecule::~Molecule( void ){ - int i; - CutoffGroup* cg; - vector::iterator iter; - - if( myAtoms != NULL ){ - for(i=0; i::iterator iterCutoff; - Atom* cutoffAtom; - vector::iterator iterAtom; - int atomIndex; - - nAtoms = theInit.nAtoms; - nMembers = nAtoms; - nBonds = theInit.nBonds; - nBends = theInit.nBends; - nTorsions = theInit.nTorsions; - nRigidBodies = theInit.nRigidBodies; - nOriented = theInit.nOriented; + void Molecule::addConstraintPair(ConstraintPair* cp) { + if (std::find(constraintPairs_.begin(), constraintPairs_.end(), cp) == constraintPairs_.end()) { + constraintPairs_.push_back(cp); + } - myAtoms = theInit.myAtoms; - myBonds = theInit.myBonds; - myBends = theInit.myBends; - myTorsions = theInit.myTorsions; - myRigidBodies = theInit.myRigidBodies; + } - myIntegrableObjects = theInit.myIntegrableObjects; + void Molecule::addConstraintElem(ConstraintElem* cp) { + if (std::find(constraintElems_.begin(), constraintElems_.end(), cp) == constraintElems_.end()) { + constraintElems_.push_back(cp); + } - for (int i = 0; i < myRigidBodies.size(); i++) - myRigidBodies[i]->calcRefCoords(); + } - myCutoffGroups = theInit.myCutoffGroups; - nCutoffGroups = myCutoffGroups.size(); + void Molecule::complete() { + + std::set rigidAtoms; + RigidBody* rb; + std::vector::iterator rbIter; -} + + for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { + rigidAtoms.insert(rb->getBeginAtomIter(), rb->getEndAtomIter()); + } -void Molecule::calcForces( void ){ - - int i; - double com[3]; + Atom* atom; + AtomIterator ai; + for (atom = beginAtom(ai); atom != NULL; atom = nextAtom(ai)) { + + if (rigidAtoms.find(*ai) == rigidAtoms.end()) { + //if an atom does not belong to a rigid body, it is an integrable object + integrableObjects_.push_back(*ai); + } + } - for(i=0; iupdateAtoms(); - } + //find all free atoms (which do not belong to rigid bodies) + //performs the "difference" operation from set theory, the output range contains a copy of every + //element that is contained in [allAtoms.begin(), allAtoms.end()) and not contained in + //[rigidAtoms.begin(), rigidAtoms.end()). + //std::set_difference(allAtoms.begin(), allAtoms.end(), rigidAtoms.begin(), rigidAtoms.end(), + // std::back_inserter(integrableObjects_)); - for(i=0; icalc_forces(); + //if (integrableObjects_.size() != allAtoms.size() - rigidAtoms.size()) { + // //Some atoms in rigidAtoms are not in allAtoms, something must be wrong + // sprintf(painCave.errMsg, "Atoms in rigidbody are not in the atom list of the same molecule"); + // + // painCave.isFatal = 1; + // simError(); + //} + for (rb = beginRigidBody(rbIter); rb != NULL; rb = nextRigidBody(rbIter)) { + integrableObjects_.push_back(rb); + } + //integrableObjects_.insert(integrableObjects_.end(), rigidBodies_.begin(), rigidBodies_.end()); } - for(i=0; icalc_forces(); - } + RealType Molecule::getMass() { + StuntDouble* sd; + std::vector::iterator i; + RealType mass = 0.0; - for(i=0; icalc_forces(); - } - - // Rigid Body forces and torques are done after the fortran force loop - -} - + for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ + mass += sd->getMass(); + } -double Molecule::getPotential( void ){ - - int i; - double myPot = 0.0; + return mass; - for(i=0; iupdateAtoms(); } - - for(i=0; iget_potential(); - } - for(i=0; iget_potential(); - } + Vector3d Molecule::getCom() { + StuntDouble* sd; + std::vector::iterator i; + Vector3d com; + RealType totalMass = 0; + RealType mass; + + for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ + mass = sd->getMass(); + totalMass += mass; + com += sd->getPos() * mass; + } - for(i=0; iget_potential(); - } + com /= totalMass; - return myPot; -} - -void Molecule::printMe( void ){ - - int i; - - for(i=0; iprintMe(); + return com; } - for(i=0; iprintMe(); - } + void Molecule::moveCom(const Vector3d& delta) { + StuntDouble* sd; + std::vector::iterator i; + + for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ + sd->setPos(sd->getPos() + delta); + } - for(i=0; iprintMe(); } -} - -void Molecule::moveCOM(double delta[3]){ - double aPos[3]; - int i, j; - - for(i=0; igetPos( aPos ); - - for (j=0; j< 3; j++) - aPos[j] += delta[j]; - - myIntegrableObjects[i]->setPos( aPos ); + Vector3d Molecule::getComVel() { + StuntDouble* sd; + std::vector::iterator i; + Vector3d velCom; + RealType totalMass = 0; + RealType mass; + + for (sd = beginIntegrableObject(i); sd != NULL; sd = nextIntegrableObject(i)){ + mass = sd->getMass(); + totalMass += mass; + velCom += sd->getVel() * mass; } - } - for(i=0; igetPos( aPos ); - - for (j=0; j< 3; j++) - aPos[j] += delta[j]; - - myRigidBodies[i]->setPos( aPos ); - } -} - -void Molecule::atoms2rigidBodies( void ) { - int i; - for (i = 0; i < myRigidBodies.size(); i++) { - myRigidBodies[i]->calcForcesAndTorques(); + return velCom; } -} -void Molecule::getCOM( double COM[3] ) { + RealType Molecule::getPotential() { - double mass, mtot; - double aPos[3]; - int i, j; + Bond* bond; + Bend* bend; + Torsion* torsion; + Molecule::BondIterator bondIter;; + Molecule::BendIterator bendIter; + Molecule::TorsionIterator torsionIter; - for (j=0; j<3; j++) - COM[j] = 0.0; + RealType potential = 0.0; - mtot = 0.0; + for (bond = beginBond(bondIter); bond != NULL; bond = nextBond(bondIter)) { + potential += bond->getPotential(); + } - for (i=0; i < myIntegrableObjects.size(); i++) { - if (myIntegrableObjects[i] != NULL) { + for (bend = beginBend(bendIter); bend != NULL; bend = nextBend(bendIter)) { + potential += bend->getPotential(); + } - mass = myIntegrableObjects[i]->getMass(); - mtot += mass; - - myIntegrableObjects[i]->getPos( aPos ); - - for( j = 0; j < 3; j++) - COM[j] += aPos[j] * mass; - + for (torsion = beginTorsion(torsionIter); torsion != NULL; torsion = nextTorsion(torsionIter)) { + potential += torsion->getPotential(); } - } - for (j = 0; j < 3; j++) - COM[j] /= mtot; -} + return potential; -double Molecule::getCOMvel( double COMvel[3] ) { - - double mass, mtot; - double aVel[3]; - int i, j; - - - for (j=0; j<3; j++) - COMvel[j] = 0.0; - - mtot = 0.0; - - for (i=0; i < myIntegrableObjects.size(); i++) { - if (myIntegrableObjects[i] != NULL) { - - mass = myIntegrableObjects[i]->getMass(); - mtot += mass; - - myIntegrableObjects[i]->getVel(aVel); - - for (j=0; j<3; j++) - COMvel[j] += aVel[j]*mass; - - } } - for (j=0; j<3; j++) - COMvel[j] /= mtot; - - return mtot; - -} - -double Molecule::getTotalMass() -{ - - double totalMass; - - totalMass = 0; - for(int i =0; i < myIntegrableObjects.size(); i++){ - totalMass += myIntegrableObjects[i]->getMass(); + std::ostream& operator <<(std::ostream& o, Molecule& mol) { + o << std::endl; + o << "Molecule " << mol.getGlobalIndex() << "has: " << std::endl; + o << mol.getNAtoms() << " atoms" << std::endl; + o << mol.getNBonds() << " bonds" << std::endl; + o << mol.getNBends() << " bends" << std::endl; + o << mol.getNTorsions() << " torsions" << std::endl; + o << mol.getNRigidBodies() << " rigid bodies" << std::endl; + o << mol.getNIntegrableObjects() << "integrable objects" << std::endl; + o << mol.getNCutoffGroups() << "cutoff groups" << std::endl; + o << mol.getNConstraintPairs() << "constraint pairs" << std::endl; + return o; } - return totalMass; -} +}//end namespace oopse