--- trunk/src/brains/SimInfo.cpp 2004/09/24 04:16:43 2 +++ branches/development/src/brains/SimInfo.cpp 2012/06/07 12:53:46 1750 @@ -1,624 +1,1258 @@ -#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. Redistributions of source code must retain the above copyright + * notice, this list of conditions and the following disclaimer. + * + * 2. 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. + * + * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your + * research, please cite the appropriate papers when you publish your + * work. Good starting points are: + * + * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). + * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). + */ + +/** + * @file SimInfo.cpp + * @author tlin + * @date 11/02/2004 + * @version 1.0 + */ -#include -using namespace std; +#include +#include +#include -#include "SimInfo.hpp" -#define __C -#include "fSimulation.h" -#include "simError.h" - -#include "fortranWrappers.hpp" - -#include "MatVec3.h" - +#include "brains/SimInfo.hpp" +#include "math/Vector3.hpp" +#include "primitives/Molecule.hpp" +#include "primitives/StuntDouble.hpp" +#include "utils/MemoryUtils.hpp" +#include "utils/simError.h" +#include "selection/SelectionManager.hpp" +#include "io/ForceFieldOptions.hpp" +#include "brains/ForceField.hpp" +#include "nonbonded/SwitchingFunction.hpp" #ifdef IS_MPI -#include "mpiSimulation.hpp" +#include #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; -} +using namespace std; +namespace OpenMD { + + SimInfo::SimInfo(ForceField* ff, Globals* simParams) : + forceField_(ff), simParams_(simParams), + ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), + nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), + nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0), + nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), + nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), + nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), + calcBoxDipole_(false), useAtomicVirial_(true) { + + 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; + + vector components = simParams->getComponents(); + + for (vector::iterator i = components.begin(); i !=components.end(); ++i) { + molStamp = (*i)->getMoleculeStamp(); + nMolWithSameStamp = (*i)->getNMol(); + + addMoleculeStamp(molStamp, nMolWithSameStamp); + + //calculate atoms in molecules + nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; + + //calculate atoms in cutoff groups + int nAtomsInGroups = 0; + int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); + + for (int j=0; j < nCutoffGroupsInStamp; j++) { + cgStamp = molStamp->getCutoffGroupStamp(j); + nAtomsInGroups += cgStamp->getNMembers(); + } + + nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; + + //calculate atoms in rigid bodies + int nAtomsInRigidBodies = 0; + int nRigidBodiesInStamp = molStamp->getNRigidBodies(); + + for (int j=0; j < nRigidBodiesInStamp; j++) { + rbStamp = molStamp->getRigidBodyStamp(j); + nAtomsInRigidBodies += rbStamp->getNMembers(); + } + + nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; + nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; + + } + + //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 -SimInfo* currentInfo; + nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; + + //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_; + + nGlobalMols_ = molStampIds_.size(); + molToProcMap_.resize(nGlobalMols_); + } + + SimInfo::~SimInfo() { + map::iterator i; + for (i = molecules_.begin(); i != molecules_.end(); ++i) { + delete i->second; + } + molecules_.clear(); + + delete sman_; + delete simParams_; + delete forceField_; + } -SimInfo::SimInfo(){ - 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; - - haveRcut = 0; - haveRsw = 0; - boxIsInit = 0; + bool SimInfo::addMolecule(Molecule* mol) { + MoleculeIterator i; + + i = molecules_.find(mol->getGlobalIndex()); + if (i == molecules_.end() ) { + + molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); + + nAtoms_ += mol->getNAtoms(); + nBonds_ += mol->getNBonds(); + nBends_ += mol->getNBends(); + nTorsions_ += mol->getNTorsions(); + nInversions_ += mol->getNInversions(); + nRigidBodies_ += mol->getNRigidBodies(); + nIntegrableObjects_ += mol->getNIntegrableObjects(); + nCutoffGroups_ += mol->getNCutoffGroups(); + nConstraints_ += mol->getNConstraintPairs(); + + addInteractionPairs(mol); + + return true; + } else { + return false; + } + } - resetTime = 1e99; + bool SimInfo::removeMolecule(Molecule* mol) { + MoleculeIterator i; + i = molecules_.find(mol->getGlobalIndex()); - orthoRhombic = 0; - orthoTolerance = 1E-6; - useInitXSstate = true; + if (i != molecules_.end() ) { - usePBC = 0; - useLJ = 0; - useSticky = 0; - useCharges = 0; - useDipoles = 0; - useReactionField = 0; - useGB = 0; - useEAM = 0; - useSolidThermInt = 0; - useLiquidThermInt = 0; + assert(mol == i->second); + + nAtoms_ -= mol->getNAtoms(); + nBonds_ -= mol->getNBonds(); + nBends_ -= mol->getNBends(); + nTorsions_ -= mol->getNTorsions(); + nInversions_ -= mol->getNInversions(); + nRigidBodies_ -= mol->getNRigidBodies(); + nIntegrableObjects_ -= mol->getNIntegrableObjects(); + nCutoffGroups_ -= mol->getNCutoffGroups(); + nConstraints_ -= mol->getNConstraintPairs(); - haveCutoffGroups = false; + removeInteractionPairs(mol); + molecules_.erase(mol->getGlobalIndex()); - excludes = Exclude::Instance(); - - myConfiguration = new SimState(); + delete mol; + + return true; + } else { + return false; + } + } - has_minimizer = false; - the_minimizer =NULL; + + Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { + i = molecules_.begin(); + return i == molecules_.end() ? NULL : i->second; + } - ngroup = 0; + Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { + ++i; + return i == molecules_.end() ? NULL : i->second; + } - wrapMeSimInfo( this ); -} + void SimInfo::calcNdf() { + int ndf_local, nfq_local; + MoleculeIterator i; + vector::iterator j; + vector::iterator k; -SimInfo::~SimInfo(){ + Molecule* mol; + StuntDouble* integrableObject; + Atom* atom; - delete myConfiguration; + ndf_local = 0; + nfq_local = 0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { - map::iterator i; - - for(i = properties.begin(); i != properties.end(); i++) - delete (*i).second; + ndf_local += 3; -} + if (integrableObject->isDirectional()) { + if (integrableObject->isLinear()) { + ndf_local += 2; + } else { + ndf_local += 3; + } + } + } + for (atom = mol->beginFluctuatingCharge(k); atom != NULL; + atom = mol->nextFluctuatingCharge(k)) { + if (atom->isFluctuatingCharge()) { + nfq_local++; + } + } + } + + ndfLocal_ = ndf_local; -void SimInfo::setBox(double newBox[3]) { - - int i, j; - double tempMat[3][3]; + // n_constraints is local, so subtract them on each processor + ndf_local -= nConstraints_; - for(i=0; i<3; i++) - for (j=0; j<3; j++) tempMat[i][j] = 0.0;; +#ifdef IS_MPI + MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); +#else + ndf_ = ndf_local; + nGlobalFluctuatingCharges_ = nfq_local; +#endif - tempMat[0][0] = newBox[0]; - tempMat[1][1] = newBox[1]; - tempMat[2][2] = newBox[2]; + // nZconstraints_ is global, as are the 3 COM translations for the + // entire system: + ndf_ = ndf_ - 3 - nZconstraint_; - setBoxM( tempMat ); + } -} - -void SimInfo::setBoxM( double theBox[3][3] ){ + int SimInfo::getFdf() { +#ifdef IS_MPI + MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + fdf_ = fdf_local; +#endif + return fdf_; + } - 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); - - if( !boxIsInit ) boxIsInit = 1; - - 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]; + unsigned int SimInfo::getNLocalCutoffGroups(){ + int nLocalCutoffAtoms = 0; + Molecule* mol; + MoleculeIterator mi; + CutoffGroup* cg; + Molecule::CutoffGroupIterator ci; + + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { + nLocalCutoffAtoms += cg->getNumAtom(); + + } } + + return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; } + + void SimInfo::calcNdfRaw() { + int ndfRaw_local; - setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic); - -} - + MoleculeIterator i; + vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; -void SimInfo::getBoxM (double theBox[3][3]) { + // 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)) { - int i, j; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; -} + 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 + } -void SimInfo::scaleBox(double scale) { - double theBox[3][3]; - int i, j; + void SimInfo::calcNdfTrans() { + int ndfTrans_local; - // cerr << "Scaling box by " << scale << "\n"; + ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; - setBoxM(theBox); +#ifdef IS_MPI + MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndfTrans_ = ndfTrans_local; +#endif -} + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; + + } -void SimInfo::calcHmatInv( void ) { - - int oldOrtho; - int i,j; - double smallDiag; - double tol; - double sanity[3][3]; - - invertMat3( Hmat, HmatInv ); - - // check to see if Hmat is orthorhombic - - oldOrtho = orthoRhombic; + void SimInfo::addInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); + vector::iterator bondIter; + vector::iterator bendIter; + vector::iterator torsionIter; + vector::iterator inversionIter; + Bond* bond; + Bend* bend; + Torsion* torsion; + Inversion* inversion; + int a; + int b; + int c; + int d; - 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; + // atomGroups can be used to add special interaction maps between + // groups of atoms that are in two separate rigid bodies. + // However, most site-site interactions between two rigid bodies + // are probably not special, just the ones between the physically + // bonded atoms. Interactions *within* a single rigid body should + // always be excluded. These are done at the bottom of this + // function. - 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; - } + map > atomGroups; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; + + for (integrableObject = mol->beginIntegrableObject(ii); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + vector atoms = rb->getAtoms(); + set rigidAtoms; + for (int i = 0; i < static_cast(atoms.size()); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < static_cast(atoms.size()); ++i) { + atomGroups.insert(map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); } - } - } + } + + for (bond= mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { - if( oldOrtho != orthoRhombic ){ + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); - 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(); + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + } else { + excludedInteractions_.addPair(a, b); + } } - 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::calcBoxL( void ){ + for (bend= mol->beginBend(bendIter); bend != NULL; + bend = mol->nextBend(bendIter)) { - double dx, dy, dz, dsq; + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(b, c); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(b, c); + } - // boxVol = Determinant of Hmat + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(a, c); + } else { + excludedInteractions_.addPair(a, c); + } + } - boxVol = matDet3( Hmat ); + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; + torsion = mol->nextTorsion(torsionIter)) { - // 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]; + a = torsion->getAtomA()->getGlobalIndex(); + b = torsion->getAtomB()->getGlobalIndex(); + c = torsion->getAtomC()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); - // 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]; + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(b, c); + oneTwoInteractions_.addPair(c, d); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(b, c); + excludedInteractions_.addPair(c, d); + } + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(a, c); + oneThreeInteractions_.addPair(b, d); + } else { + excludedInteractions_.addPair(a, c); + excludedInteractions_.addPair(b, d); + } - // 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]; + if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { + oneFourInteractions_.addPair(a, d); + } else { + excludedInteractions_.addPair(a, d); + } + } - //calculate the max cutoff - maxCutoff = calcMaxCutOff(); - - checkCutOffs(); + for (inversion= mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { -} + a = inversion->getAtomA()->getGlobalIndex(); + b = inversion->getAtomB()->getGlobalIndex(); + c = inversion->getAtomC()->getGlobalIndex(); + d = inversion->getAtomD()->getGlobalIndex(); + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.addPair(a, b); + oneTwoInteractions_.addPair(a, c); + oneTwoInteractions_.addPair(a, d); + } else { + excludedInteractions_.addPair(a, b); + excludedInteractions_.addPair(a, c); + excludedInteractions_.addPair(a, d); + } -double SimInfo::calcMaxCutOff(){ + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.addPair(b, c); + oneThreeInteractions_.addPair(b, d); + oneThreeInteractions_.addPair(c, d); + } else { + excludedInteractions_.addPair(b, c); + excludedInteractions_.addPair(b, d); + excludedInteractions_.addPair(c, d); + } + } - double ri[3], rj[3], rk[3]; - double rij[3], rjk[3], rki[3]; - double minDist; + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + vector atoms = rb->getAtoms(); + for (int i = 0; i < static_cast(atoms.size()) -1 ; ++i) { + for (int j = i + 1; j < static_cast(atoms.size()); ++j) { + a = atoms[i]->getGlobalIndex(); + b = atoms[j]->getGlobalIndex(); + excludedInteractions_.addPair(a, b); + } + } + } - 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]; + void SimInfo::removeInteractionPairs(Molecule* mol) { + ForceFieldOptions& options_ = forceField_->getForceFieldOptions(); + vector::iterator bondIter; + vector::iterator bendIter; + vector::iterator torsionIter; + vector::iterator inversionIter; + Bond* bond; + Bend* bend; + Torsion* torsion; + Inversion* inversion; + int a; + int b; + int c; + int d; - rk[0] = Hmat[0][2]; - rk[1] = Hmat[1][2]; - rk[2] = Hmat[2][2]; + map > atomGroups; + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; - crossProduct3(ri, rj, rij); - distXY = dotProduct3(rk,rij) / norm3(rij); + for (integrableObject = mol->beginIntegrableObject(ii); + integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + vector atoms = rb->getAtoms(); + set rigidAtoms; + for (int i = 0; i < static_cast(atoms.size()); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < static_cast(atoms.size()); ++i) { + atomGroups.insert(map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } - crossProduct3(rj,rk, rjk); - distYZ = dotProduct3(ri,rjk) / norm3(rjk); + for (bond= mol->beginBond(bondIter); bond != NULL; + bond = mol->nextBond(bondIter)) { + + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + } else { + excludedInteractions_.removePair(a, b); + } + } - crossProduct3(rk,ri, rki); - distZX = dotProduct3(rj,rki) / norm3(rki); + for (bend= mol->beginBend(bendIter); bend != NULL; + bend = mol->nextBend(bendIter)) { - minDist = min(min(distXY, distYZ), distZX); - return minDist/2; - -} + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); + + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(b, c); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.removePair(b, c); + } -void SimInfo::wrapVector( double thePos[3] ){ + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(a, c); + } else { + excludedInteractions_.removePair(a, c); + } + } - int i; - double scaled[3]; + for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; + torsion = mol->nextTorsion(torsionIter)) { - if( !orthoRhombic ){ - // calc the scaled coordinates. + a = torsion->getAtomA()->getGlobalIndex(); + b = torsion->getAtomB()->getGlobalIndex(); + c = torsion->getAtomC()->getGlobalIndex(); + d = torsion->getAtomD()->getGlobalIndex(); + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(b, c); + oneTwoInteractions_.removePair(c, d); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.removePair(b, c); + excludedInteractions_.removePair(c, d); + } - 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); + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(a, c); + oneThreeInteractions_.removePair(b, d); + } else { + excludedInteractions_.removePair(a, c); + excludedInteractions_.removePair(b, d); + } - } - 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]; - } - -} + if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) { + oneFourInteractions_.removePair(a, d); + } else { + excludedInteractions_.removePair(a, d); + } + } + for (inversion= mol->beginInversion(inversionIter); inversion != NULL; + inversion = mol->nextInversion(inversionIter)) { -int SimInfo::getNDF(){ - int ndf_local; + a = inversion->getAtomA()->getGlobalIndex(); + b = inversion->getAtomB()->getGlobalIndex(); + c = inversion->getAtomC()->getGlobalIndex(); + d = inversion->getAtomD()->getGlobalIndex(); - 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; + if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) { + oneTwoInteractions_.removePair(a, b); + oneTwoInteractions_.removePair(a, c); + oneTwoInteractions_.removePair(a, d); + } else { + excludedInteractions_.removePair(a, b); + excludedInteractions_.removePair(a, c); + excludedInteractions_.removePair(a, d); + } + + if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) { + oneThreeInteractions_.removePair(b, c); + oneThreeInteractions_.removePair(b, d); + oneThreeInteractions_.removePair(c, d); + } else { + excludedInteractions_.removePair(b, c); + excludedInteractions_.removePair(b, d); + excludedInteractions_.removePair(c, d); + } } + + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { + vector atoms = rb->getAtoms(); + for (int i = 0; i < static_cast(atoms.size()) -1 ; ++i) { + for (int j = i + 1; j < static_cast(atoms.size()); ++j) { + a = atoms[i]->getGlobalIndex(); + b = atoms[j]->getGlobalIndex(); + excludedInteractions_.removePair(a, b); + } + } + } + } + + + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { + int curStampId; + + //index from 0 + curStampId = moleculeStamps_.size(); - // n_constraints is local, so subtract them on each processor: + moleculeStamps_.push_back(molStamp); + molStampIds_.insert(molStampIds_.end(), nmol, curStampId); + } - ndf_local -= n_constraints; -#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 - nZconstraints; - - return ndf; -} - -int SimInfo::getNDFraw() { - int ndfRaw_local; - - // Raw degrees of freedom that we have to set - ndfRaw_local = 0; - - 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; - } + /** + * update + * + * Performs the global checks and variable settings after the + * objects have been created. + * + */ + void SimInfo::update() { + setupSimVariables(); + calcNdf(); + calcNdfRaw(); + calcNdfTrans(); } + + /** + * getSimulatedAtomTypes + * + * Returns an STL set of AtomType* that are actually present in this + * simulation. Must query all processors to assemble this information. + * + */ + set SimInfo::getSimulatedAtomTypes() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + set atomTypes; + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for(atom = mol->beginAtom(ai); atom != NULL; + atom = mol->nextAtom(ai)) { + atomTypes.insert(atom->getAtomType()); + } + } + #ifdef IS_MPI - MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfRaw = ndfRaw_local; -#endif - return ndfRaw; -} + // loop over the found atom types on this processor, and add their + // numerical idents to a vector: + + vector foundTypes; + set::iterator i; + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) + foundTypes.push_back( (*i)->getIdent() ); -int SimInfo::getNDFtranslational() { - int ndfTrans_local; + // count_local holds the number of found types on this processor + int count_local = foundTypes.size(); - ndfTrans_local = 3 * integrableObjects.size() - n_constraints; + int nproc = MPI::COMM_WORLD.Get_size(); + // we need arrays to hold the counts and displacement vectors for + // all processors + vector counts(nproc, 0); + vector disps(nproc, 0); -#ifdef IS_MPI - MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfTrans = ndfTrans_local; -#endif + // fill the counts array + MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], + 1, MPI::INT); + + // use the processor counts to compute the displacement array + disps[0] = 0; + int totalCount = counts[0]; + for (int iproc = 1; iproc < nproc; iproc++) { + disps[iproc] = disps[iproc-1] + counts[iproc-1]; + totalCount += counts[iproc]; + } - ndfTrans = ndfTrans - 3 - nZconstraints; + // we need a (possibly redundant) set of all found types: + vector ftGlobal(totalCount); + + // now spray out the foundTypes to all the other processors: + MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, + &ftGlobal[0], &counts[0], &disps[0], + MPI::INT); - return ndfTrans; -} + vector::iterator j; -int SimInfo::getTotIntegrableObjects() { - int nObjs_local; - int nObjs; + // foundIdents is a stl set, so inserting an already found ident + // will have no effect. + set foundIdents; - nObjs_local = integrableObjects.size(); - - -#ifdef IS_MPI - MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - nObjs = nObjs_local; + for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) + foundIdents.insert((*j)); + + // now iterate over the foundIdents and get the actual atom types + // that correspond to these: + set::iterator it; + for (it = foundIdents.begin(); it != foundIdents.end(); ++it) + atomTypes.insert( forceField_->getAtomType((*it)) ); + #endif - - return nObjs; -} - -void SimInfo::refreshSim(){ - - simtype fInfo; - int isError; - int n_global; - int* excl; - - fInfo.dielect = 0.0; - - if( useDipoles ){ - if( useReactionField )fInfo.dielect = dielectric; + return atomTypes; } - fInfo.SIM_uses_PBC = usePBC; - //fInfo.SIM_uses_LJ = 0; - fInfo.SIM_uses_LJ = useLJ; - fInfo.SIM_uses_sticky = useSticky; - //fInfo.SIM_uses_sticky = 0; - fInfo.SIM_uses_charges = useCharges; - fInfo.SIM_uses_dipoles = useDipoles; - //fInfo.SIM_uses_dipoles = 0; - fInfo.SIM_uses_RF = useReactionField; - //fInfo.SIM_uses_RF = 0; - fInfo.SIM_uses_GB = useGB; - fInfo.SIM_uses_EAM = useEAM; + void SimInfo::setupSimVariables() { + useAtomicVirial_ = simParams_->getUseAtomicVirial(); + // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true + calcBoxDipole_ = false; + if ( simParams_->haveAccumulateBoxDipole() ) + if ( simParams_->getAccumulateBoxDipole() ) { + calcBoxDipole_ = true; + } + + set::iterator i; + set atomTypes; + atomTypes = getSimulatedAtomTypes(); + int usesElectrostatic = 0; + int usesMetallic = 0; + int usesDirectional = 0; + int usesFluctuatingCharges = 0; + //loop over all of the atom types + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + usesElectrostatic |= (*i)->isElectrostatic(); + usesMetallic |= (*i)->isMetal(); + usesDirectional |= (*i)->isDirectional(); + usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); + } + +#ifdef IS_MPI + int temp; + temp = usesDirectional; + MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = usesMetallic; + MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = usesElectrostatic; + MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - n_exclude = excludes->getSize(); - excl = excludes->getFortranArray(); - -#ifdef IS_MPI - n_global = mpiSim->getNAtomsGlobal(); + temp = usesFluctuatingCharges; + MPI_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); #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 - setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, - &nGlobalExcludes, globalExcludes, molMembershipArray, - &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); - if( isError ){ + usesDirectionalAtoms_ = usesDirectional; + usesMetallicAtoms_ = usesMetallic; + usesElectrostaticAtoms_ = usesElectrostatic; + usesFluctuatingCharges_ = usesFluctuatingCharges; + +#endif - sprintf( painCave.errMsg, - "There was an error setting the simulation information in fortran.\n" ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); + requiresPrepair_ = usesMetallicAtoms_ ? true : false; + requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; + requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false; } - -#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(); -} -void SimInfo::setDefaultRcut( double theRcut ){ - - haveRcut = 1; - rCut = theRcut; - rList = rCut + 1.0; - - notifyFortranCutOffs( &rCut, &rSw, &rList ); -} -void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ + vector SimInfo::getGlobalAtomIndices() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; - rSw = theRsw; - setDefaultRcut( theRcut ); -} - - -void SimInfo::checkCutOffs( void ){ - - if( boxIsInit ){ + vector GlobalAtomIndices(getNAtoms(), 0); - //we need to check cutOffs against the box - - if( rCut > maxCutoff ){ - 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; - painCave.isFatal = 1; - 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(); + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + + for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); + } + } + return GlobalAtomIndices; } - -} -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()){ + vector SimInfo::getGlobalGroupIndices() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; + + vector GlobalGroupIndices; - delete (*result).second; - (*result).second = prop; + for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + //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)) { + GlobalGroupIndices.push_back(cg->getGlobalIndex()); + } + } + return GlobalGroupIndices; } - else{ - properties[prop->getID()] = prop; - } + void SimInfo::prepareTopology() { + int nExclude, nOneTwo, nOneThree, nOneFour; + + //calculate mass ratio of cutoff group + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; + Molecule::AtomIterator ai; + Atom* atom; + RealType totalMass; + + /** + * The mass factor is the relative mass of an atom to the total + * mass of the cutoff group it belongs to. By default, all atoms + * are their own cutoff groups, and therefore have mass factors of + * 1. We need some special handling for massless atoms, which + * will be treated as carrying the entire mass of the cutoff + * group. + */ + massFactors_.clear(); + massFactors_.resize(getNAtoms(), 1.0); -} + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for (cg = mol->beginCutoffGroup(ci); cg != NULL; + cg = mol->nextCutoffGroup(ci)) { -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; -} + totalMass = cg->getMass(); + for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { + // Check for massless groups - set mfact to 1 if true + if (totalMass != 0) + massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; + else + massFactors_[atom->getLocalIndex()] = 1.0; + } + } + } + // Build the identArray_ -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(); - + identArray_.clear(); + 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()); + } + } + + //scan topology - // Fix the silly fortran indexing problem + nExclude = excludedInteractions_.getSize(); + nOneTwo = oneTwoInteractions_.getSize(); + nOneThree = oneThreeInteractions_.getSize(); + nOneFour = oneFourInteractions_.getSize(); + + int* excludeList = excludedInteractions_.getPairList(); + int* oneTwoList = oneTwoInteractions_.getPairList(); + int* oneThreeList = oneThreeInteractions_.getPairList(); + int* oneFourList = oneFourInteractions_.getPairList(); + + topologyDone_ = true; + } + + void SimInfo::addProperty(GenericData* genData) { + properties_.addProperty(genData); + } + + void SimInfo::removeProperty(const string& propName) { + properties_.removeProperty(propName); + } + + void SimInfo::clearProperties() { + properties_.clearProperties(); + } + + vector SimInfo::getPropertyNames() { + return properties_.getPropertyNames(); + } + + vector SimInfo::getProperties() { + return properties_.getProperties(); + } + + GenericData* SimInfo::getPropertyByName(const 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; + CutoffGroup* cg; + SimInfo::MoleculeIterator mi; + Molecule::RigidBodyIterator rbIter; + Molecule::AtomIterator atomIter; + Molecule::CutoffGroupIterator cgIter; + + 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_); + } + + for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { + cg->setSnapshotManager(sman_); + } + } + + } + + Vector3d SimInfo::getComVel(){ + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d comVel(0.0); + RealType totalMass = 0.0; + + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + RealType mass = mol->getMass(); + totalMass += mass; + comVel += mass * mol->getComVel(); + } + #ifdef IS_MPI - numAtom = mpiSim->getNAtomsGlobal(); -#else - numAtom = n_atoms; + RealType tmpMass = totalMass; + Vector3d tmpComVel(comVel); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); #endif - for (int i = 0; i < numAtom; i++) - FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); - - 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)){ + comVel /= totalMass; - totalMass = myCutoffGroup->getMass(); + return comVel; + } + + Vector3d SimInfo::getCom(){ + SimInfo::MoleculeIterator i; + Molecule* mol; + + Vector3d com(0.0); + RealType totalMass = 0.0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + RealType mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + } + +#ifdef IS_MPI + RealType tmpMass = totalMass; + Vector3d tmpCom(com); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); +#endif + + com /= totalMass; + + return com; + + } + + ostream& operator <<(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; - for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); - cutoffAtom != NULL; - cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ - mfact.push_back(cutoffAtom->getMass()/totalMass); + + RealType totalMass = 0.0; + + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + RealType mass = mol->getMass(); + totalMass += mass; + com += mass * mol->getCom(); + comVel += mass * mol->getComVel(); } + +#ifdef IS_MPI + RealType tmpMass = totalMass; + Vector3d tmpCom(com); + Vector3d tmpComVel(comVel); + MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,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){ + + + RealType xx = 0.0; + RealType yy = 0.0; + RealType zz = 0.0; + RealType xy = 0.0; + RealType xz = 0.0; + RealType 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); + + RealType 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_REALTYPE,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,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); + + RealType 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_REALTYPE,MPI_SUM, MPI_COMM_WORLD); +#endif + + return angularMomentum; + } + + StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { + return IOIndexToIntegrableObject.at(index); + } + + void SimInfo::setIOIndexToIntegrableObject(const vector& v) { + IOIndexToIntegrableObject= v; + } + + /* Returns the Volume of the simulation based on a ellipsoid with semi-axes + based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3 + where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to + V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536. + */ + void SimInfo::getGyrationalVolume(RealType &volume){ + Mat3x3d intTensor; + RealType det; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor,dummyAngMom); + + det = intTensor.determinant(); + sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(det); + return; + } + + void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ + Mat3x3d intTensor; + Vector3d dummyAngMom; + RealType sysconstants; + RealType geomCnst; + + geomCnst = 3.0/2.0; + /* Get the inertial tensor and angular momentum for free*/ + getInertiaTensor(intTensor,dummyAngMom); + + detI = intTensor.determinant(); + sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; + volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI); + return; + } +/* + void SimInfo::setStuntDoubleFromGlobalIndex(vector v) { + assert( v.size() == nAtoms_ + nRigidBodies_); + sdByGlobalIndex_ = v; } + + StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { + //assert(index < nAtoms_ + nRigidBodies_); + return sdByGlobalIndex_.at(index); + } +*/ + 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; } -} +}//end namespace OpenMD +