--- trunk/src/brains/SimInfo.cpp 2004/10/22 22:54:01 143 +++ trunk/src/brains/SimInfo.cpp 2005/05/31 22:31:54 557 @@ -1,639 +1,1092 @@ -#include -#include -#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 SimInfo.cpp + * @author tlin + * @date 11/02/2004 + * @version 1.0 + */ -#include -using namespace std; +#include +#include #include "brains/SimInfo.hpp" -#define __C -#include "brains/fSimulation.h" -#include "utils/simError.h" -#include "UseTheForce/DarkSide/simulation_interface.h" +#include "math/Vector3.hpp" +#include "primitives/Molecule.hpp" +#include "UseTheForce/doForces_interface.h" #include "UseTheForce/notifyCutoffs_interface.h" +#include "utils/MemoryUtils.hpp" +#include "utils/simError.h" +#include "selection/SelectionManager.hpp" -//#include "UseTheForce/fortranWrappers.hpp" - -#include "math/MatVec3.h" - #ifdef IS_MPI -#include "brains/mpiSimulation.hpp" -#endif +#include "UseTheForce/mpiComponentPlan.h" +#include "UseTheForce/DarkSide/simParallel_interface.h" +#endif -inline double roundMe( double x ){ - return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); -} - -inline double min( double a, double b ){ - return (a < b ) ? a : b; -} +namespace oopse { -SimInfo* currentInfo; + SimInfo::SimInfo(MakeStamps* stamps, std::vector >& molStampPairs, + ForceField* ff, Globals* simParams) : + stamps_(stamps), forceField_(ff), simParams_(simParams), + ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), + nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), + nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), + nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), + nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), + sman_(NULL), fortranInitialized_(false) { -SimInfo::SimInfo(){ + + std::vector >::iterator i; + MoleculeStamp* molStamp; + int nMolWithSameStamp; + int nCutoffAtoms = 0; // number of atoms belong to cutoff groups + int nGroups = 0; //total cutoff groups defined in meta-data file + CutoffGroupStamp* cgStamp; + RigidBodyStamp* rbStamp; + int nRigidAtoms = 0; + + for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { + molStamp = i->first; + nMolWithSameStamp = i->second; + + addMoleculeStamp(molStamp, nMolWithSameStamp); - n_constraints = 0; - nZconstraints = 0; - n_oriented = 0; - n_dipoles = 0; - ndf = 0; - ndfRaw = 0; - nZconstraints = 0; - the_integrator = NULL; - setTemp = 0; - thermalTime = 0.0; - currentTime = 0.0; - rCut = 0.0; - rSw = 0.0; + //calculate atoms in molecules + nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; - haveRcut = 0; - haveRsw = 0; - boxIsInit = 0; - - resetTime = 1e99; - orthoRhombic = 0; - orthoTolerance = 1E-6; - useInitXSstate = true; + //calculate atoms in cutoff groups + int nAtomsInGroups = 0; + int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); + + for (int j=0; j < nCutoffGroupsInStamp; j++) { + cgStamp = molStamp->getCutoffGroup(j); + nAtomsInGroups += cgStamp->getNMembers(); + } - usePBC = 0; - useDirectionalAtoms = 0; - useLennardJones = 0; - useElectrostatics = 0; - useCharges = 0; - useDipoles = 0; - useSticky = 0; - useGayBerne = 0; - useEAM = 0; - useShapes = 0; - useFLARB = 0; + nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; - useSolidThermInt = 0; - useLiquidThermInt = 0; + //calculate atoms in rigid bodies + int nAtomsInRigidBodies = 0; + int nRigidBodiesInStamp = molStamp->getNRigidBodies(); + + for (int j=0; j < nRigidBodiesInStamp; j++) { + rbStamp = molStamp->getRigidBody(j); + nAtomsInRigidBodies += rbStamp->getNMembers(); + } - haveCutoffGroups = false; + nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; + nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; + + } - excludes = Exclude::Instance(); + //every free atom (atom does not belong to cutoff groups) is a cutoff group + //therefore the total number of cutoff groups in the system is equal to + //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data + //file plus the number of cutoff groups defined in meta-data file + nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; - myConfiguration = new SimState(); + //every free atom (atom does not belong to rigid bodies) is an integrable object + //therefore the total number of integrable objects in the system is equal to + //the total number of atoms minus number of atoms belong to rigid body defined in meta-data + //file plus the number of rigid bodies defined in meta-data file + nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; - has_minimizer = false; - the_minimizer =NULL; + nGlobalMols_ = molStampIds_.size(); - ngroup = 0; +#ifdef IS_MPI + molToProcMap_.resize(nGlobalMols_); +#endif -} + } + SimInfo::~SimInfo() { + std::map::iterator i; + for (i = molecules_.begin(); i != molecules_.end(); ++i) { + delete i->second; + } + molecules_.clear(); + + delete stamps_; + delete sman_; + delete simParams_; + delete forceField_; + } -SimInfo::~SimInfo(){ + int SimInfo::getNGlobalConstraints() { + int nGlobalConstraints; +#ifdef IS_MPI + MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, + MPI_COMM_WORLD); +#else + nGlobalConstraints = nConstraints_; +#endif + return nGlobalConstraints; + } - delete myConfiguration; + bool SimInfo::addMolecule(Molecule* mol) { + MoleculeIterator i; - map::iterator i; - - for(i = properties.begin(); i != properties.end(); i++) - delete (*i).second; + i = molecules_.find(mol->getGlobalIndex()); + if (i == molecules_.end() ) { -} - -void SimInfo::setBox(double newBox[3]) { - - int i, j; - double tempMat[3][3]; + molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); + + nAtoms_ += mol->getNAtoms(); + nBonds_ += mol->getNBonds(); + nBends_ += mol->getNBends(); + nTorsions_ += mol->getNTorsions(); + nRigidBodies_ += mol->getNRigidBodies(); + nIntegrableObjects_ += mol->getNIntegrableObjects(); + nCutoffGroups_ += mol->getNCutoffGroups(); + nConstraints_ += mol->getNConstraintPairs(); - for(i=0; i<3; i++) - for (j=0; j<3; j++) tempMat[i][j] = 0.0;; + addExcludePairs(mol); + + return true; + } else { + return false; + } + } - tempMat[0][0] = newBox[0]; - tempMat[1][1] = newBox[1]; - tempMat[2][2] = newBox[2]; + bool SimInfo::removeMolecule(Molecule* mol) { + MoleculeIterator i; + i = molecules_.find(mol->getGlobalIndex()); - setBoxM( tempMat ); + if (i != molecules_.end() ) { -} + assert(mol == i->second); + + nAtoms_ -= mol->getNAtoms(); + nBonds_ -= mol->getNBonds(); + nBends_ -= mol->getNBends(); + nTorsions_ -= mol->getNTorsions(); + nRigidBodies_ -= mol->getNRigidBodies(); + nIntegrableObjects_ -= mol->getNIntegrableObjects(); + nCutoffGroups_ -= mol->getNCutoffGroups(); + nConstraints_ -= mol->getNConstraintPairs(); -void SimInfo::setBoxM( double theBox[3][3] ){ - - int i, j; - double FortranHmat[9]; // to preserve compatibility with Fortran the - // ordering in the array is as follows: - // [ 0 3 6 ] - // [ 1 4 7 ] - // [ 2 5 8 ] - double FortranHmatInv[9]; // the inverted Hmat (for Fortran); + removeExcludePairs(mol); + molecules_.erase(mol->getGlobalIndex()); - if( !boxIsInit ) boxIsInit = 1; + delete mol; + + return true; + } else { + return false; + } - for(i=0; i < 3; i++) - for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; - - calcBoxL(); - calcHmatInv(); - for(i=0; i < 3; i++) { - for (j=0; j < 3; j++) { - FortranHmat[3*j + i] = Hmat[i][j]; - FortranHmatInv[3*j + i] = HmatInv[i][j]; - } - } + } - setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic); - -} - + + Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { + i = molecules_.begin(); + return i == molecules_.end() ? NULL : i->second; + } -void SimInfo::getBoxM (double theBox[3][3]) { + Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { + ++i; + return i == molecules_.end() ? NULL : i->second; + } - int i, j; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; -} + void SimInfo::calcNdf() { + int ndf_local; + MoleculeIterator i; + std::vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; -void SimInfo::scaleBox(double scale) { - double theBox[3][3]; - int i, j; + ndf_local = 0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - // cerr << "Scaling box by " << scale << "\n"; + ndf_local += 3; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; + if (integrableObject->isDirectional()) { + if (integrableObject->isLinear()) { + ndf_local += 2; + } else { + ndf_local += 3; + } + } + + }//end for (integrableObject) + }// end for (mol) + + // n_constraints is local, so subtract them on each processor + ndf_local -= nConstraints_; - setBoxM(theBox); +#ifdef IS_MPI + MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndf_ = ndf_local; +#endif -} + // nZconstraints_ is global, as are the 3 COM translations for the + // entire system: + ndf_ = ndf_ - 3 - nZconstraint_; -void SimInfo::calcHmatInv( void ) { - - int oldOrtho; - int i,j; - double smallDiag; - double tol; - double sanity[3][3]; + } - invertMat3( Hmat, HmatInv ); + void SimInfo::calcNdfRaw() { + int ndfRaw_local; - // check to see if Hmat is orthorhombic - - oldOrtho = orthoRhombic; + MoleculeIterator i; + std::vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; - smallDiag = fabs(Hmat[0][0]); - if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); - if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); - tol = smallDiag * orthoTolerance; + // Raw degrees of freedom that we have to set + ndfRaw_local = 0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - orthoRhombic = 1; - - for (i = 0; i < 3; i++ ) { - for (j = 0 ; j < 3; j++) { - if (i != j) { - if (orthoRhombic) { - if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; - } + ndfRaw_local += 3; + + if (integrableObject->isDirectional()) { + if (integrableObject->isLinear()) { + ndfRaw_local += 2; + } else { + ndfRaw_local += 3; + } + } + } } + +#ifdef IS_MPI + MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndfRaw_ = ndfRaw_local; +#endif } - if( oldOrtho != orthoRhombic ){ - - if( orthoRhombic ) { - sprintf( painCave.errMsg, - "OOPSE is switching from the default Non-Orthorhombic\n" - "\tto the faster Orthorhombic periodic boundary computations.\n" - "\tThis is usually a good thing, but if you wan't the\n" - "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" - "\tvariable ( currently set to %G ) smaller.\n", - orthoTolerance); - painCave.severity = OOPSE_INFO; - simError(); - } - else { - sprintf( painCave.errMsg, - "OOPSE is switching from the faster Orthorhombic to the more\n" - "\tflexible Non-Orthorhombic periodic boundary computations.\n" - "\tThis is usually because the box has deformed under\n" - "\tNPTf integration. If you wan't to live on the edge with\n" - "\tthe Orthorhombic computations, make the orthoBoxTolerance\n" - "\tvariable ( currently set to %G ) larger.\n", - orthoTolerance); - painCave.severity = OOPSE_WARNING; - simError(); - } - } -} + void SimInfo::calcNdfTrans() { + int ndfTrans_local; -void SimInfo::calcBoxL( void ){ + ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; - double dx, dy, dz, dsq; - // boxVol = Determinant of Hmat +#ifdef IS_MPI + MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndfTrans_ = ndfTrans_local; +#endif - boxVol = matDet3( Hmat ); + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; + + } - // boxLx - - dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; - dsq = dx*dx + dy*dy + dz*dz; - boxL[0] = sqrt( dsq ); - //maxCutoff = 0.5 * boxL[0]; + void SimInfo::addExcludePairs(Molecule* mol) { + std::vector::iterator bondIter; + std::vector::iterator bendIter; + std::vector::iterator torsionIter; + Bond* bond; + Bend* bend; + Torsion* torsion; + int a; + int b; + int c; + int d; + + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + exclude_.addPair(a, b); + } - // boxLy - - dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; - dsq = dx*dx + dy*dy + dz*dz; - boxL[1] = sqrt( dsq ); - //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1]; + for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(b, c); + } - // boxLz - - dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; - dsq = dx*dx + dy*dy + dz*dz; - boxL[2] = sqrt( dsq ); - //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2]; + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { + a = torsion->getAtomA()->getGlobalIndex(); + b = torsion->getAtomB()->getGlobalIndex(); + c = torsion->getAtomC()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); - //calculate the max cutoff - maxCutoff = calcMaxCutOff(); - - checkCutOffs(); + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(a, d); + exclude_.addPair(b, c); + exclude_.addPair(b, d); + exclude_.addPair(c, d); + } -} + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + std::vector atoms = rb->getAtoms(); + for (int i = 0; i < atoms.size() -1 ; ++i) { + for (int j = i + 1; j < atoms.size(); ++j) { + a = atoms[i]->getGlobalIndex(); + b = atoms[j]->getGlobalIndex(); + exclude_.addPair(a, b); + } + } + } + } -double SimInfo::calcMaxCutOff(){ - - double ri[3], rj[3], rk[3]; - double rij[3], rjk[3], rki[3]; - double minDist; - - ri[0] = Hmat[0][0]; - ri[1] = Hmat[1][0]; - ri[2] = Hmat[2][0]; - - rj[0] = Hmat[0][1]; - rj[1] = Hmat[1][1]; - rj[2] = Hmat[2][1]; - - rk[0] = Hmat[0][2]; - rk[1] = Hmat[1][2]; - rk[2] = Hmat[2][2]; + void SimInfo::removeExcludePairs(Molecule* mol) { + std::vector::iterator bondIter; + std::vector::iterator bendIter; + std::vector::iterator torsionIter; + Bond* bond; + Bend* bend; + Torsion* torsion; + int a; + int b; + int c; + int d; - crossProduct3(ri, rj, rij); - distXY = dotProduct3(rk,rij) / norm3(rij); + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + exclude_.removePair(a, b); + } - crossProduct3(rj,rk, rjk); - distYZ = dotProduct3(ri,rjk) / norm3(rjk); + for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); - crossProduct3(rk,ri, rki); - distZX = dotProduct3(rj,rki) / norm3(rki); + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(b, c); + } - minDist = min(min(distXY, distYZ), distZX); - return minDist/2; - -} + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { + a = torsion->getAtomA()->getGlobalIndex(); + b = torsion->getAtomB()->getGlobalIndex(); + c = torsion->getAtomC()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); -void SimInfo::wrapVector( double thePos[3] ){ + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(a, d); + exclude_.removePair(b, c); + exclude_.removePair(b, d); + exclude_.removePair(c, d); + } - int i; - double scaled[3]; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + std::vector atoms = rb->getAtoms(); + for (int i = 0; i < atoms.size() -1 ; ++i) { + for (int j = i + 1; j < atoms.size(); ++j) { + a = atoms[i]->getGlobalIndex(); + b = atoms[j]->getGlobalIndex(); + exclude_.removePair(a, b); + } + } + } - if( !orthoRhombic ){ - // calc the scaled coordinates. - - - matVecMul3(HmatInv, thePos, scaled); - - for(i=0; i<3; i++) - scaled[i] -= roundMe(scaled[i]); - - // calc the wrapped real coordinates from the wrapped scaled coordinates - - matVecMul3(Hmat, scaled, thePos); - } - else{ - // calc the scaled coordinates. - - for(i=0; i<3; i++) - scaled[i] = thePos[i]*HmatInv[i][i]; - - // wrap the scaled coordinates - - for(i=0; i<3; i++) - scaled[i] -= roundMe(scaled[i]); - - // calc the wrapped real coordinates from the wrapped scaled coordinates - - for(i=0; i<3; i++) - thePos[i] = scaled[i]*Hmat[i][i]; - } - -} -int SimInfo::getNDF(){ - int ndf_local; + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { + int curStampId; - ndf_local = 0; - - for(int i = 0; i < integrableObjects.size(); i++){ - ndf_local += 3; - if (integrableObjects[i]->isDirectional()) { - if (integrableObjects[i]->isLinear()) - ndf_local += 2; - else - ndf_local += 3; - } + //index from 0 + curStampId = moleculeStamps_.size(); + + moleculeStamps_.push_back(molStamp); + molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } - // n_constraints is local, so subtract them on each processor: + void SimInfo::update() { - ndf_local -= n_constraints; + setupSimType(); #ifdef IS_MPI - MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndf = ndf_local; + setupFortranParallel(); #endif - // nZconstraints is global, as are the 3 COM translations for the - // entire system: + setupFortranSim(); - ndf = ndf - 3 - nZconstraints; + //setup fortran force field + /** @deprecate */ + int isError = 0; + initFortranFF( &fInfo_.SIM_uses_RF , &isError ); + if(isError){ + sprintf( painCave.errMsg, + "ForceField error: There was an error initializing the forceField in fortran.\n" ); + painCave.isFatal = 1; + simError(); + } + + + setupCutoff(); - return ndf; -} + calcNdf(); + calcNdfRaw(); + calcNdfTrans(); -int SimInfo::getNDFraw() { - int ndfRaw_local; + fortranInitialized_ = true; + } - // Raw degrees of freedom that we have to set - ndfRaw_local = 0; + std::set SimInfo::getUniqueAtomTypes() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + std::set atomTypes; - for(int i = 0; i < integrableObjects.size(); i++){ - ndfRaw_local += 3; - if (integrableObjects[i]->isDirectional()) { - if (integrableObjects[i]->isLinear()) - ndfRaw_local += 2; - else - ndfRaw_local += 3; + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + atomTypes.insert(atom->getAtomType()); + } + } + + return atomTypes; } + + void SimInfo::setupSimType() { + std::set::iterator i; + std::set atomTypes; + atomTypes = getUniqueAtomTypes(); -#ifdef IS_MPI - MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfRaw = ndfRaw_local; -#endif + int useLennardJones = 0; + int useElectrostatic = 0; + int useEAM = 0; + int useCharge = 0; + int useDirectional = 0; + int useDipole = 0; + int useGayBerne = 0; + int useSticky = 0; + int useStickyPower = 0; + int useShape = 0; + int useFLARB = 0; //it is not in AtomType yet + int useDirectionalAtom = 0; + int useElectrostatics = 0; + //usePBC and useRF are from simParams + int usePBC = simParams_->getPBC(); + int useRF = simParams_->getUseRF(); - return ndfRaw; -} + //loop over all of the atom types + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + useLennardJones |= (*i)->isLennardJones(); + useElectrostatic |= (*i)->isElectrostatic(); + useEAM |= (*i)->isEAM(); + useCharge |= (*i)->isCharge(); + useDirectional |= (*i)->isDirectional(); + useDipole |= (*i)->isDipole(); + useGayBerne |= (*i)->isGayBerne(); + useSticky |= (*i)->isSticky(); + useStickyPower |= (*i)->isStickyPower(); + useShape |= (*i)->isShape(); + } -int SimInfo::getNDFtranslational() { - int ndfTrans_local; + if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { + useDirectionalAtom = 1; + } - ndfTrans_local = 3 * integrableObjects.size() - n_constraints; + if (useCharge || useDipole) { + useElectrostatics = 1; + } +#ifdef IS_MPI + int temp; -#ifdef IS_MPI - MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfTrans = ndfTrans_local; -#endif + temp = usePBC; + MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - ndfTrans = ndfTrans - 3 - nZconstraints; + temp = useDirectionalAtom; + MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - return ndfTrans; -} + temp = useLennardJones; + MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -int SimInfo::getTotIntegrableObjects() { - int nObjs_local; - int nObjs; + temp = useElectrostatics; + MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - nObjs_local = integrableObjects.size(); + temp = useCharge; + MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useDipole; + MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -#ifdef IS_MPI - MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - nObjs = nObjs_local; -#endif + temp = useSticky; + MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useStickyPower; + MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useGayBerne; + MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - return nObjs; -} + temp = useEAM; + MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -void SimInfo::refreshSim(){ + temp = useShape; + MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - simtype fInfo; - int isError; - int n_global; - int* excl; + temp = useFLARB; + MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - fInfo.dielect = 0.0; + temp = useRF; + MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + +#endif - if( useDipoles ){ - if( useReactionField )fInfo.dielect = dielectric; - } + fInfo_.SIM_uses_PBC = usePBC; + fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; + fInfo_.SIM_uses_LennardJones = useLennardJones; + fInfo_.SIM_uses_Electrostatics = useElectrostatics; + fInfo_.SIM_uses_Charges = useCharge; + fInfo_.SIM_uses_Dipoles = useDipole; + fInfo_.SIM_uses_Sticky = useSticky; + fInfo_.SIM_uses_StickyPower = useStickyPower; + fInfo_.SIM_uses_GayBerne = useGayBerne; + fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_Shapes = useShape; + fInfo_.SIM_uses_FLARB = useFLARB; + fInfo_.SIM_uses_RF = useRF; - fInfo.SIM_uses_PBC = usePBC; + if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { - if (useSticky || useDipoles || useGayBerne || useShapes) { - useDirectionalAtoms = 1; - fInfo.SIM_uses_DirectionalAtoms = useDirectionalAtoms; - } + if (simParams_->haveDielectric()) { + fInfo_.dielect = simParams_->getDielectric(); + } else { + sprintf(painCave.errMsg, + "SimSetup Error: No Dielectric constant was set.\n" + "\tYou are trying to use Reaction Field without" + "\tsetting a dielectric constant!\n"); + painCave.isFatal = 1; + simError(); + } + + } else { + fInfo_.dielect = 0.0; + } - fInfo.SIM_uses_LennardJones = useLennardJones; - - if (useCharges || useDipoles) { - useElectrostatics = 1; - fInfo.SIM_uses_Electrostatics = useElectrostatics; } - fInfo.SIM_uses_Charges = useCharges; - fInfo.SIM_uses_Dipoles = useDipoles; - fInfo.SIM_uses_Sticky = useSticky; - fInfo.SIM_uses_GayBerne = useGayBerne; - fInfo.SIM_uses_EAM = useEAM; - fInfo.SIM_uses_Shapes = useShapes; - fInfo.SIM_uses_FLARB = useFLARB; - fInfo.SIM_uses_RF = useReactionField; + void SimInfo::setupFortranSim() { + int isError; + int nExclude; + std::vector fortranGlobalGroupMembership; + + nExclude = exclude_.getSize(); + isError = 0; - n_exclude = excludes->getSize(); - excl = excludes->getFortranArray(); - -#ifdef IS_MPI - n_global = mpiSim->getNAtomsGlobal(); -#else - n_global = n_atoms; -#endif - - isError = 0; - - getFortranGroupArrays(this, FglobalGroupMembership, mfact); - //it may not be a good idea to pass the address of first element in vector - //since c++ standard does not require vector to be stored continuously in meomory - //Most of the compilers will organize the memory of vector continuously - setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, - &nGlobalExcludes, globalExcludes, molMembershipArray, - &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); + //globalGroupMembership_ is filled by SimCreator + for (int i = 0; i < nGlobalAtoms_; i++) { + fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); + } - if( isError ){ + //calculate mass ratio of cutoff group + std::vector mfact; + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; + Molecule::AtomIterator ai; + Atom* atom; + double totalMass; + + //to avoid memory reallocation, reserve enough space for mfact + mfact.reserve(getNCutoffGroups()); - sprintf( painCave.errMsg, - "There was an error setting the simulation information in fortran.\n" ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); - } - -#ifdef IS_MPI - sprintf( checkPointMsg, - "succesfully sent the simulation information to fortran.\n"); - MPIcheckPoint(); -#endif // is_mpi - - this->ndf = this->getNDF(); - this->ndfRaw = this->getNDFraw(); - this->ndfTrans = this->getNDFtranslational(); -} + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { -void SimInfo::setDefaultRcut( double theRcut ){ - - haveRcut = 1; - rCut = theRcut; - rList = rCut + 1.0; - - notifyFortranCutoffs( &rCut, &rSw, &rList ); -} + totalMass = cg->getMass(); + for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { + mfact.push_back(atom->getMass()/totalMass); + } -void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ + } + } - rSw = theRsw; - setDefaultRcut( theRcut ); -} + //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) + std::vector identArray; + //to avoid memory reallocation, reserve enough space identArray + identArray.reserve(getNAtoms()); + + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + identArray.push_back(atom->getIdent()); + } + } -void SimInfo::checkCutOffs( void ){ - - if( boxIsInit ){ + //fill molMembershipArray + //molMembershipArray is filled by SimCreator + std::vector molMembershipArray(nGlobalAtoms_); + for (int i = 0; i < nGlobalAtoms_; i++) { + molMembershipArray[i] = globalMolMembership_[i] + 1; + } - //we need to check cutOffs against the box - - if( rCut > maxCutoff ){ + //setup fortran simulation + int nGlobalExcludes = 0; + int* globalExcludes = NULL; + int* excludeList = exclude_.getExcludeList(); + setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , + &nGlobalExcludes, globalExcludes, &molMembershipArray[0], + &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); + + if( isError ){ + sprintf( painCave.errMsg, - "cutoffRadius is too large for the current periodic box.\n" - "\tCurrent Value of cutoffRadius = %G at time %G\n " - "\tThis is larger than half of at least one of the\n" - "\tperiodic box vectors. Right now, the Box matrix is:\n" - "\n" - "\t[ %G %G %G ]\n" - "\t[ %G %G %G ]\n" - "\t[ %G %G %G ]\n", - rCut, currentTime, - Hmat[0][0], Hmat[0][1], Hmat[0][2], - Hmat[1][0], Hmat[1][1], Hmat[1][2], - Hmat[2][0], Hmat[2][1], Hmat[2][2]); - painCave.severity = OOPSE_ERROR; + "There was an error setting the simulation information in fortran.\n" ); painCave.isFatal = 1; + painCave.severity = OOPSE_ERROR; simError(); - } - } else { - // initialize this stuff before using it, OK? - sprintf( painCave.errMsg, - "Trying to check cutoffs without a box.\n" - "\tOOPSE should have better programmers than that.\n" ); - painCave.severity = OOPSE_ERROR; - painCave.isFatal = 1; - simError(); - } - -} + } -void SimInfo::addProperty(GenericData* prop){ - - map::iterator result; - result = properties.find(prop->getID()); - - //we can't simply use properties[prop->getID()] = prop, - //it will cause memory leak if we already contain a propery which has the same name of prop - - if(result != properties.end()){ - - delete (*result).second; - (*result).second = prop; - +#ifdef IS_MPI + sprintf( checkPointMsg, + "succesfully sent the simulation information to fortran.\n"); + MPIcheckPoint(); +#endif // is_mpi } - else{ - properties[prop->getID()] = prop; - } +#ifdef IS_MPI + void SimInfo::setupFortranParallel() { -} + //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex + std::vector localToGlobalAtomIndex(getNAtoms(), 0); + std::vector localToGlobalCutoffGroupIndex; + SimInfo::MoleculeIterator mi; + Molecule::AtomIterator ai; + Molecule::CutoffGroupIterator ci; + Molecule* mol; + Atom* atom; + CutoffGroup* cg; + mpiSimData parallelData; + int isError; -GenericData* SimInfo::getProperty(const string& propName){ - - map::iterator result; - - //string lowerCaseName = (); - - result = properties.find(propName); - - if(result != properties.end()) - return (*result).second; - else - return NULL; -} + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + //local index(index in DataStorge) of atom is important + for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; + } -void SimInfo::getFortranGroupArrays(SimInfo* info, - vector& FglobalGroupMembership, - vector& mfact){ - - Molecule* myMols; - Atom** myAtoms; - int numAtom; - double mtot; - int numMol; - int numCutoffGroups; - CutoffGroup* myCutoffGroup; - vector::iterator iterCutoff; - Atom* cutoffAtom; - vector::iterator iterAtom; - int atomIndex; - double totalMass; - - mfact.clear(); - FglobalGroupMembership.clear(); - + //local index of cutoff group is trivial, it only depends on the order of travesing + for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { + localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); + } + + } - // Fix the silly fortran indexing problem -#ifdef IS_MPI - numAtom = mpiSim->getNAtomsGlobal(); -#else - numAtom = n_atoms; -#endif - for (int i = 0; i < numAtom; i++) - FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); - + //fill up mpiSimData struct + parallelData.nMolGlobal = getNGlobalMolecules(); + parallelData.nMolLocal = getNMolecules(); + parallelData.nAtomsGlobal = getNGlobalAtoms(); + parallelData.nAtomsLocal = getNAtoms(); + parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); + parallelData.nGroupsLocal = getNCutoffGroups(); + parallelData.myNode = worldRank; + MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); - myMols = info->molecules; - numMol = info->n_mol; - for(int i = 0; i < numMol; i++){ - numCutoffGroups = myMols[i].getNCutoffGroups(); - for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); - myCutoffGroup != NULL; - myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ + //pass mpiSimData struct and index arrays to fortran + setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), + &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), + &localToGlobalCutoffGroupIndex[0], &isError); - totalMass = myCutoffGroup->getMass(); - - for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); - cutoffAtom != NULL; - cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ - mfact.push_back(cutoffAtom->getMass()/totalMass); - } + if (isError) { + sprintf(painCave.errMsg, + "mpiRefresh errror: fortran didn't like something we gave it.\n"); + painCave.isFatal = 1; + simError(); } + + sprintf(checkPointMsg, " mpiRefresh successful.\n"); + MPIcheckPoint(); + + } -} +#endif + + double SimInfo::calcMaxCutoffRadius() { + + + std::set atomTypes; + std::set::iterator i; + std::vector cutoffRadius; + + //get the unique atom types + atomTypes = getUniqueAtomTypes(); + + //query the max cutoff radius among these atom types + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); + } + + double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); +#ifdef IS_MPI + //pick the max cutoff radius among the processors +#endif + + return maxCutoffRadius; + } + + void SimInfo::getCutoff(double& rcut, double& rsw) { + + if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + + if (!simParams_->haveRcut()){ + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the cutoffRadius.\n" + "\tOOPSE will use a default value of 15.0 angstroms" + "\tfor the cutoffRadius.\n"); + painCave.isFatal = 0; + simError(); + rcut = 15.0; + } else{ + rcut = simParams_->getRcut(); + } + + if (!simParams_->haveRsw()){ + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\n" + "\t0.95 * cutoffRadius for the switchingRadius\n"); + painCave.isFatal = 0; + simError(); + rsw = 0.95 * rcut; + } else{ + rsw = simParams_->getRsw(); + } + + } else { + // if charge, dipole or reaction field is not used and the cutofff radius is not specified in + //meta-data file, the maximum cutoff radius calculated from forcefiled will be used + + if (simParams_->haveRcut()) { + rcut = simParams_->getRcut(); + } else { + //set cutoff radius to the maximum cutoff radius based on atom types in the whole system + rcut = calcMaxCutoffRadius(); + } + + if (simParams_->haveRsw()) { + rsw = simParams_->getRsw(); + } else { + rsw = rcut; + } + + } + } + + void SimInfo::setupCutoff() { + getCutoff(rcut_, rsw_); + double rnblist = rcut_ + 1; // skin of neighbor list + + //Pass these cutoff radius etc. to fortran. This function should be called once and only once + notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); + } + + void SimInfo::addProperty(GenericData* genData) { + properties_.addProperty(genData); + } + + void SimInfo::removeProperty(const std::string& propName) { + properties_.removeProperty(propName); + } + + void SimInfo::clearProperties() { + properties_.clearProperties(); + } + + std::vector SimInfo::getPropertyNames() { + return properties_.getPropertyNames(); + } + + std::vector SimInfo::getProperties() { + return properties_.getProperties(); + } + + GenericData* SimInfo::getPropertyByName(const std::string& propName) { + return properties_.getPropertyByName(propName); + } + + void SimInfo::setSnapshotManager(SnapshotManager* sman) { + if (sman_ == sman) { + return; + } + delete sman_; + sman_ = sman; + + Molecule* mol; + RigidBody* rb; + Atom* atom; + SimInfo::MoleculeIterator mi; + Molecule::RigidBodyIterator rbIter; + Molecule::AtomIterator atomIter;; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { + atom->setSnapshotManager(sman_); + } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { + rb->setSnapshotManager(sman_); + } + } + + } + + Vector3d SimInfo::getComVel(){ + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d comVel(0.0); + double totalMass = 0.0; + + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + double mass = mol->getMass(); + totalMass += mass; + comVel += mass * mol->getComVel(); + } + +#ifdef IS_MPI + double tmpMass = totalMass; + Vector3d tmpComVel(comVel); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + comVel /= totalMass; + + return comVel; + } + + Vector3d SimInfo::getCom(){ + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d com(0.0); + double totalMass = 0.0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + double mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + } + +#ifdef IS_MPI + double tmpMass = totalMass; + Vector3d tmpCom(com); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + com /= totalMass; + + return com; + + } + + std::ostream& operator <<(std::ostream& o, SimInfo& info) { + + return o; + } + + + /* + Returns center of mass and center of mass velocity in one function call. + */ + + void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ + SimInfo::MoleculeIterator i; + Molecule* mol; + + + double totalMass = 0.0; + + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + double mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + comVel += mass * mol->getComVel(); + } + +#ifdef IS_MPI + double tmpMass = totalMass; + Vector3d tmpCom(com); + Vector3d tmpComVel(comVel); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + com /= totalMass; + comVel /= totalMass; + } + + /* + Return intertia tensor for entire system and angular momentum Vector. + + + [ Ixx -Ixy -Ixz ] + J =| -Iyx Iyy -Iyz | + [ -Izx -Iyz Izz ] + */ + + void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ + + + double xx = 0.0; + double yy = 0.0; + double zz = 0.0; + double xy = 0.0; + double xz = 0.0; + double yz = 0.0; + Vector3d com(0.0); + Vector3d comVel(0.0); + + getComAll(com, comVel); + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d thisq(0.0); + Vector3d thisv(0.0); + + double thisMass = 0.0; + + + + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + + thisq = mol->getCom()-com; + thisv = mol->getComVel()-comVel; + thisMass = mol->getMass(); + // Compute moment of intertia coefficients. + xx += thisq[0]*thisq[0]*thisMass; + yy += thisq[1]*thisq[1]*thisMass; + zz += thisq[2]*thisq[2]*thisMass; + + // compute products of intertia + xy += thisq[0]*thisq[1]*thisMass; + xz += thisq[0]*thisq[2]*thisMass; + yz += thisq[1]*thisq[2]*thisMass; + + angularMomentum += cross( thisq, thisv ) * thisMass; + + } + + + inertiaTensor(0,0) = yy + zz; + inertiaTensor(0,1) = -xy; + inertiaTensor(0,2) = -xz; + inertiaTensor(1,0) = -xy; + inertiaTensor(1,1) = xx + zz; + inertiaTensor(1,2) = -yz; + inertiaTensor(2,0) = -xz; + inertiaTensor(2,1) = -yz; + inertiaTensor(2,2) = xx + yy; + +#ifdef IS_MPI + Mat3x3d tmpI(inertiaTensor); + Vector3d tmpAngMom; + MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return; + } + + //Returns the angular momentum of the system + Vector3d SimInfo::getAngularMomentum(){ + + Vector3d com(0.0); + Vector3d comVel(0.0); + Vector3d angularMomentum(0.0); + + getComAll(com,comVel); + + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d thisr(0.0); + Vector3d thisp(0.0); + + double thisMass; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + thisMass = mol->getMass(); + thisr = mol->getCom()-com; + thisp = (mol->getComVel()-comVel)*thisMass; + + angularMomentum += cross( thisr, thisp ); + + } + +#ifdef IS_MPI + Vector3d tmpAngMom; + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return angularMomentum; + } + + +}//end namespace oopse +