--- trunk/OOPSE-2.0/src/brains/SimInfo.cpp 2004/09/24 16:27:58 1492 +++ branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp 2004/12/03 20:30:07 1842 @@ -1,624 +1,826 @@ -#include -#include -#include +/* + * Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project + * + * Contact: oopse@oopse.org + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * All we ask is that proper credit is given for our work, which includes + * - but is not limited to - adding the above copyright notice to the beginning + * of your source code files, and to any copyright notice that you may distribute + * with programs based on this work. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + */ -#include -using namespace std; +/** + * @file SimInfo.cpp + * @author tlin + * @date 11/02/2004 + * @version 1.0 + */ +#include +#include #include "brains/SimInfo.hpp" -#define __C -#include "brains/fSimulation.h" +#include "primitives/Molecule.hpp" +#include "UseTheForce/doForces_interface.h" +#include "UseTheForce/notifyCutoffs_interface.h" +#include "utils/MemoryUtils.hpp" #include "utils/simError.h" -#include "UseTheForce/fortranWrappers.hpp" +namespace oopse { -#include "math/MatVec3.h" +SimInfo::SimInfo(std::vector >& molStampPairs, + ForceField* ff, Globals* simParams) : + 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) { -#ifdef IS_MPI -#include "brains/mpiSimulation.hpp" -#endif + + 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); -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; -} + //calculate atoms in molecules + nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; -SimInfo* currentInfo; -SimInfo::SimInfo(){ + //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(); + } - 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; + nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; - haveRcut = 0; - haveRsw = 0; - boxIsInit = 0; - - resetTime = 1e99; + //calculate atoms in rigid bodies + int nAtomsInRigidBodies = 0; + int nRigidBodiesInStamp = molStamp->getNCutoffGroups(); + + for (int j=0; j < nRigidBodiesInStamp; j++) { + rbStamp = molStamp->getRigidBody(j); + nAtomsInRigidBodies += rbStamp->getNMembers(); + } - orthoRhombic = 0; - orthoTolerance = 1E-6; - useInitXSstate = true; + nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; + nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; + + } - usePBC = 0; - useLJ = 0; - useSticky = 0; - useCharges = 0; - useDipoles = 0; - useReactionField = 0; - useGB = 0; - useEAM = 0; - useSolidThermInt = 0; - useLiquidThermInt = 0; + //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; - haveCutoffGroups = false; + //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_; - excludes = Exclude::Instance(); + nGlobalMols_ = molStampIds_.size(); - myConfiguration = new SimState(); +#ifdef IS_MPI + molToProcMap_.resize(nGlobalMols_); +#endif + +} - has_minimizer = false; - the_minimizer =NULL; +SimInfo::~SimInfo() { + //MemoryUtils::deleteVectorOfPointer(molecules_); - ngroup = 0; + MemoryUtils::deleteVectorOfPointer(moleculeStamps_); + + delete sman_; + delete simParams_; + delete forceField_; - wrapMeSimInfo( this ); } -SimInfo::~SimInfo(){ +bool SimInfo::addMolecule(Molecule* mol) { + MoleculeIterator i; - delete myConfiguration; + i = molecules_.find(mol->getGlobalIndex()); + if (i == molecules_.end() ) { - map::iterator i; - - for(i = properties.begin(); i != properties.end(); i++) - delete (*i).second; + 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->getNConstraints(); + return true; + } else { + return false; + } } -void SimInfo::setBox(double newBox[3]) { - - int i, j; - double tempMat[3][3]; +bool SimInfo::removeMolecule(Molecule* mol) { + MoleculeIterator i; + i = molecules_.find(mol->getGlobalIndex()); - for(i=0; i<3; i++) - for (j=0; j<3; j++) tempMat[i][j] = 0.0;; + if (i != molecules_.end() ) { - tempMat[0][0] = newBox[0]; - tempMat[1][1] = newBox[1]; - tempMat[2][2] = newBox[2]; + 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->getNConstraints(); - setBoxM( tempMat ); + molecules_.erase(mol->getGlobalIndex()); -} + delete mol; + + return true; + } else { + return false; + } -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); - if( !boxIsInit ) boxIsInit = 1; +} - for(i=0; i < 3; i++) - for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; - - calcBoxL(); - calcHmatInv(); + +Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { + i = molecules_.begin(); + return i == molecules_.end() ? NULL : i->second; +} - 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]; - } - } - - setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic); - +Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { + ++i; + return i == molecules_.end() ? NULL : i->second; } - -void SimInfo::getBoxM (double theBox[3][3]) { - 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; + ndf_local = 0; + + for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { + for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(j)) { -void SimInfo::scaleBox(double scale) { - double theBox[3][3]; - int i, j; + ndf_local += 3; - // cerr << "Scaling box by " << scale << "\n"; + 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_; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; +#ifdef IS_MPI + MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndf_ = ndf_local; +#endif - setBoxM(theBox); + // 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]; +void SimInfo::calcNdfRaw() { + int ndfRaw_local; - invertMat3( Hmat, HmatInv ); + MoleculeIterator i; + std::vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; - // check to see if Hmat is orthorhombic - - oldOrtho = orthoRhombic; + // 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)) { - 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; + ndfRaw_local += 3; - 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; - } - } + if (integrableObject->isDirectional()) { + if (integrableObject->isLinear()) { + ndfRaw_local += 2; + } else { + ndfRaw_local += 3; + } + } + + } } - } - - 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(); - } - } +#ifdef IS_MPI + MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndfRaw_ = ndfRaw_local; +#endif } -void SimInfo::calcBoxL( void ){ +void SimInfo::calcNdfTrans() { + int ndfTrans_local; - double dx, dy, dz, dsq; + ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; - // boxVol = Determinant of Hmat - boxVol = matDet3( Hmat ); +#ifdef IS_MPI + MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndfTrans_ = ndfTrans_local; +#endif - // 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]; + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; + +} - // 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]; +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); + } + for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); - // 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]; + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(b, c); + } - //calculate the max cutoff - maxCutoff = calcMaxCutOff(); - - checkCutOffs(); + 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(); + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(a, d); + exclude_.addPair(b, c); + exclude_.addPair(b, d); + exclude_.addPair(c, d); + } + + } +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; + + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + exclude_.removePair(a, b); + } -double SimInfo::calcMaxCutOff(){ + for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { + a = bend->getAtomA()->getGlobalIndex(); + b = bend->getAtomB()->getGlobalIndex(); + c = bend->getAtomC()->getGlobalIndex(); - double ri[3], rj[3], rk[3]; - double rij[3], rjk[3], rki[3]; - double minDist; + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(b, c); + } - ri[0] = Hmat[0][0]; - ri[1] = Hmat[1][0]; - ri[2] = Hmat[2][0]; + 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(); - rj[0] = Hmat[0][1]; - rj[1] = Hmat[1][1]; - rj[2] = Hmat[2][1]; + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(a, d); + exclude_.removePair(b, c); + exclude_.removePair(b, d); + exclude_.removePair(c, d); + } - rk[0] = Hmat[0][2]; - rk[1] = Hmat[1][2]; - rk[2] = Hmat[2][2]; - - crossProduct3(ri, rj, rij); - distXY = dotProduct3(rk,rij) / norm3(rij); +} - crossProduct3(rj,rk, rjk); - distYZ = dotProduct3(ri,rjk) / norm3(rjk); - crossProduct3(rk,ri, rki); - distZX = dotProduct3(rj,rki) / norm3(rki); +void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { + int curStampId; - minDist = min(min(distXY, distYZ), distZX); - return minDist/2; - + //index from 0 + curStampId = molStampIds_.size(); + + moleculeStamps_.push_back(molStamp); + molStampIds_.insert(molStampIds_.end(), nmol, curStampId); } -void SimInfo::wrapVector( double thePos[3] ){ +void SimInfo::update() { - int i; - double scaled[3]; + setupSimType(); - if( !orthoRhombic ){ - // calc the scaled coordinates. - +#ifdef IS_MPI + setupFortranParallel(); +#endif - 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); + setupFortranSim(); - } - else{ - // calc the scaled coordinates. + //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(); + } + - 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]; - } - + setupCutoff(); + + calcNdf(); + calcNdfRaw(); + calcNdfTrans(); + + fortranInitialized_ = true; } +std::set SimInfo::getUniqueAtomTypes() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + std::set atomTypes; -int SimInfo::getNDF(){ - int ndf_local; + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - 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; + for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { + atomTypes.insert(atom->getAtomType()); + } + } - } - // n_constraints is local, so subtract them on each processor: + return atomTypes; +} - ndf_local -= n_constraints; +void SimInfo::setupSimType() { + std::set::iterator i; + std::set atomTypes; + atomTypes = getUniqueAtomTypes(); + + 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 useShape = 0; + int useFLARB = 0; //it is not in AtomType yet + int useDirectionalAtom = 0; + int useElectrostatics = 0; + //usePBC and useRF are from simParams + bool usePBC = simParams_->getPBC(); + bool useRF = simParams_->getUseRF(); -#ifdef IS_MPI - MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndf = ndf_local; -#endif + //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(); + useShape |= (*i)->isShape(); + } - // nZconstraints is global, as are the 3 COM translations for the - // entire system: + if (useSticky || useDipole || useGayBerne || useShape) { + useDirectionalAtom = 1; + } - ndf = ndf - 3 - nZconstraints; + if (useCharge || useDipole) { + useElectrostatics = 1; + } - return ndf; -} +#ifdef IS_MPI + int temp; -int SimInfo::getNDFraw() { - int ndfRaw_local; + temp = usePBC; + MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - // Raw degrees of freedom that we have to set - ndfRaw_local = 0; + temp = useDirectionalAtom; + MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - 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; - } - } - -#ifdef IS_MPI - MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfRaw = ndfRaw_local; -#endif + temp = useLennardJones; + MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - return ndfRaw; -} + temp = useElectrostatics; + MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -int SimInfo::getNDFtranslational() { - int ndfTrans_local; + temp = useCharge; + MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - ndfTrans_local = 3 * integrableObjects.size() - n_constraints; + temp = useDipole; + MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + temp = useSticky; + MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -#ifdef IS_MPI - MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - ndfTrans = ndfTrans_local; -#endif + temp = useGayBerne; + MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - ndfTrans = ndfTrans - 3 - nZconstraints; + temp = useEAM; + MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - return ndfTrans; -} + temp = useShape; + MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -int SimInfo::getTotIntegrableObjects() { - int nObjs_local; - int nObjs; + temp = useFLARB; + MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - nObjs_local = integrableObjects.size(); + temp = useRF; + MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + +#endif + 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_GayBerne = useGayBerne; + fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_Shapes = useShape; + fInfo_.SIM_uses_FLARB = useFLARB; + fInfo_.SIM_uses_RF = useRF; -#ifdef IS_MPI - MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - nObjs = nObjs_local; -#endif + if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { + 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; + } - return nObjs; } -void SimInfo::refreshSim(){ +void SimInfo::setupFortranSim() { + int isError; + int nExclude; + std::vector fortranGlobalGroupMembership; + + nExclude = exclude_.getSize(); + isError = 0; - simtype fInfo; - int isError; - int n_global; - int* excl; + //globalGroupMembership_ is filled by SimCreator + for (int i = 0; i < nGlobalAtoms_; i++) { + fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); + } - fInfo.dielect = 0.0; + //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; - if( useDipoles ){ - if( useReactionField )fInfo.dielect = dielectric; - } + //to avoid memory reallocation, reserve enough space for mfact + mfact.reserve(getNCutoffGroups()); + + for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { + for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { - 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; + totalMass = cg->getMass(); + for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { + mfact.push_back(atom->getMass()/totalMass); + } - 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 - setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, - &nGlobalExcludes, globalExcludes, molMembershipArray, - &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); + } + } - if( isError ){ + //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()); - sprintf( painCave.errMsg, - "There was an error setting the simulation information in fortran.\n" ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); - } - + 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()); + } + } + + //fill molMembershipArray + //molMembershipArray is filled by SimCreator + std::vector molMembershipArray(nGlobalAtoms_); + for (int i = 0; i < nGlobalAtoms_; i++) { + molMembershipArray.push_back(globalMolMembership_[i] + 1); + } + + //setup fortran simulation + //gloalExcludes and molMembershipArray should go away (They are never used) + //why the hell fortran need to know molecule? + //OOPSE = Object-Obfuscated Parallel Simulation Engine + 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, + "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(); + 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 ){ +#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; - rSw = theRsw; - setDefaultRcut( theRcut ); -} + 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::checkCutOffs( void ){ - - if( boxIsInit ){ - - //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(); - } - + //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); + } + + } + + //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)); + + //pass mpiSimData struct and index arrays to fortran + setFsimParallel(parallelData, &(parallelData->nAtomsLocal), + &localToGlobalAtomIndex[0], &(parallelData->nGroupsLocal), + &localToGlobalCutoffGroupIndex[0], &isError); + + 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(); + + } -void SimInfo::addProperty(GenericData* prop){ +#endif - 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()){ +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::setupCutoff() { + double rcut_; //cutoff radius + double rsw_; //switching radius - delete (*result).second; - (*result).second = prop; - - } - else{ + 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(); + } - properties[prop->getID()] = prop; + 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_; + } + } + + 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); } -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; +void SimInfo::addProperty(GenericData* genData) { + properties_.addProperty(genData); } +void SimInfo::removeProperty(const std::string& propName) { + properties_.removeProperty(propName); +} -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(); - +void SimInfo::clearProperties() { + properties_.clearProperties(); +} - // 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); - +std::vector SimInfo::getPropertyNames() { + return properties_.getPropertyNames(); +} + +std::vector SimInfo::getProperties() { + return properties_.getProperties(); +} - 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)){ +GenericData* SimInfo::getPropertyByName(const std::string& propName) { + return properties_.getPropertyByName(propName); +} - totalMass = myCutoffGroup->getMass(); - - for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); - cutoffAtom != NULL; - cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ - mfact.push_back(cutoffAtom->getMass()/totalMass); - } - } - } +void SimInfo::setSnapshotManager(SnapshotManager* 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_); + } + } + } + +std::ostream& operator <<(std::ostream& o, SimInfo& info) { + + return o; +} + +}//end namespace oopse +