--- trunk/src/brains/SimInfo.cpp 2004/09/24 04:16:43 2 +++ trunk/src/brains/SimInfo.cpp 2006/05/17 21:51:42 963 @@ -1,624 +1,1415 @@ -#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 -#include "SimInfo.hpp" -#define __C -#include "fSimulation.h" -#include "simError.h" +#include "brains/SimInfo.hpp" +#include "math/Vector3.hpp" +#include "primitives/Molecule.hpp" +#include "UseTheForce/fCutoffPolicy.h" +#include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" +#include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" +#include "UseTheForce/DarkSide/fSwitchingFunctionType.h" +#include "UseTheForce/doForces_interface.h" +#include "UseTheForce/DarkSide/electrostatic_interface.h" +#include "UseTheForce/DarkSide/switcheroo_interface.h" +#include "utils/MemoryUtils.hpp" +#include "utils/simError.h" +#include "selection/SelectionManager.hpp" +#include "io/ForceFieldOptions.hpp" +#include "UseTheForce/ForceField.hpp" -#include "fortranWrappers.hpp" - -#include "MatVec3.h" - #ifdef IS_MPI -#include "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 { + std::set getRigidSet(int index, std::map >& container) { + std::map >::iterator i = container.find(index); + std::set result; + if (i != container.end()) { + result = i->second; + } -SimInfo* currentInfo; - -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; + return result; + } - resetTime = 1e99; + 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), + nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), + nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), + sman_(NULL), fortranInitialized_(false) { - orthoRhombic = 0; - orthoTolerance = 1E-6; - useInitXSstate = 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; + std::vector components = simParams->getComponents(); + + for (std::vector::iterator i = components.begin(); i !=components.end(); ++i) { + molStamp = (*i)->getMoleculeStamp(); + nMolWithSameStamp = (*i)->getNMol(); + + addMoleculeStamp(molStamp, nMolWithSameStamp); - usePBC = 0; - useLJ = 0; - useSticky = 0; - useCharges = 0; - useDipoles = 0; - useReactionField = 0; - useGB = 0; - useEAM = 0; - useSolidThermInt = 0; - useLiquidThermInt = 0; + //calculate atoms in molecules + nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; - haveCutoffGroups = false; + //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(); + } - excludes = Exclude::Instance(); + nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; - myConfiguration = new SimState(); + nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; - has_minimizer = false; - the_minimizer =NULL; + //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(); + } - ngroup = 0; + nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; + nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; + + } - wrapMeSimInfo( this ); -} + //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; - -SimInfo::~SimInfo(){ - - delete myConfiguration; - - map::iterator i; + //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_; - for(i = properties.begin(); i != properties.end(); i++) - delete (*i).second; + nGlobalMols_ = molStampIds_.size(); -} +#ifdef IS_MPI + molToProcMap_.resize(nGlobalMols_); +#endif -void SimInfo::setBox(double newBox[3]) { - - int i, j; - double tempMat[3][3]; + } - for(i=0; i<3; i++) - for (j=0; j<3; j++) tempMat[i][j] = 0.0;; + SimInfo::~SimInfo() { + std::map::iterator i; + for (i = molecules_.begin(); i != molecules_.end(); ++i) { + delete i->second; + } + molecules_.clear(); + + delete sman_; + delete simParams_; + delete forceField_; + } - tempMat[0][0] = newBox[0]; - tempMat[1][1] = newBox[1]; - tempMat[2][2] = newBox[2]; + 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; + } - setBoxM( tempMat ); + bool SimInfo::addMolecule(Molecule* mol) { + MoleculeIterator i; -} + i = molecules_.find(mol->getGlobalIndex()); + if (i == molecules_.end() ) { -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; + 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++) 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]; + addExcludePairs(mol); + + return true; + } else { + return false; } } - setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic); - -} - + bool SimInfo::removeMolecule(Molecule* mol) { + MoleculeIterator i; + i = molecules_.find(mol->getGlobalIndex()); -void SimInfo::getBoxM (double theBox[3][3]) { + if (i != molecules_.end() ) { - int i, j; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; -} + 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(); + removeExcludePairs(mol); + molecules_.erase(mol->getGlobalIndex()); -void SimInfo::scaleBox(double scale) { - double theBox[3][3]; - int i, j; + delete mol; + + return true; + } else { + return false; + } - // cerr << "Scaling box by " << scale << "\n"; - for(i=0; i<3; i++) - for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; + } - setBoxM(theBox); + + Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { + i = molecules_.begin(); + return i == molecules_.end() ? NULL : i->second; + } -} + Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { + ++i; + return i == molecules_.end() ? NULL : i->second; + } -void SimInfo::calcHmatInv( void ) { - - int oldOrtho; - int i,j; - double smallDiag; - double tol; - double sanity[3][3]; - invertMat3( Hmat, HmatInv ); + void SimInfo::calcNdf() { + int ndf_local; + MoleculeIterator i; + std::vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; - // check to see if Hmat is orthorhombic - - oldOrtho = orthoRhombic; + ndf_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; + ndf_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()) { + ndf_local += 2; + } else { + ndf_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(); - } - } -} + // n_constraints is local, so subtract them on each processor + ndf_local -= nConstraints_; -void SimInfo::calcBoxL( void ){ - - double dx, dy, dz, dsq; - - // boxVol = Determinant of Hmat - - boxVol = matDet3( Hmat ); - - // 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]; - - // 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]; - - - // 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]; - - //calculate the max cutoff - maxCutoff = calcMaxCutOff(); - - checkCutOffs(); - -} - - -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]; - - 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); - - minDist = min(min(distXY, distYZ), distZX); - return minDist/2; - -} - -void SimInfo::wrapVector( double thePos[3] ){ - - int i; - double scaled[3]; - - if( !orthoRhombic ){ - // calc the scaled coordinates. - +#ifdef IS_MPI + MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); +#else + ndf_ = ndf_local; +#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); + // nZconstraints_ is global, as are the 3 COM translations for the + // entire system: + ndf_ = ndf_ - 3 - nZconstraint_; } - 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; - - 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; - } - } - - // n_constraints is local, so subtract them on each processor: - - ndf_local -= n_constraints; - + int SimInfo::getFdf() { #ifdef IS_MPI - MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); #else - ndf = ndf_local; + fdf_ = fdf_local; #endif + return fdf_; + } + + void SimInfo::calcNdfRaw() { + int ndfRaw_local; - // nZconstraints is global, as are the 3 COM translations for the - // entire system: + MoleculeIterator i; + std::vector::iterator j; + Molecule* mol; + StuntDouble* integrableObject; - ndf = ndf - 3 - nZconstraints; + // 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)) { - return ndf; -} + ndfRaw_local += 3; -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; + 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); + MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); #else - ndfRaw = ndfRaw_local; + ndfRaw_ = ndfRaw_local; #endif + } - return ndfRaw; -} + void SimInfo::calcNdfTrans() { + int ndfTrans_local; -int SimInfo::getNDFtranslational() { - int ndfTrans_local; + ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; - ndfTrans_local = 3 * integrableObjects.size() - n_constraints; - #ifdef IS_MPI - MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); #else - ndfTrans = ndfTrans_local; + ndfTrans_ = ndfTrans_local; #endif - ndfTrans = ndfTrans - 3 - nZconstraints; + ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; + + } - return ndfTrans; -} + 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; -int SimInfo::getTotIntegrableObjects() { - int nObjs_local; - int nObjs; + std::map > atomGroups; - nObjs_local = integrableObjects.size(); + 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); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < atoms.size(); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < atoms.size(); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } -#ifdef IS_MPI - MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); -#else - nObjs = nObjs_local; -#endif + + + 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(); + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); - return nObjs; -} + exclude_.addPairs(rigidSetA, rigidSetB); + exclude_.addPairs(rigidSetA, rigidSetC); + exclude_.addPairs(rigidSetB, rigidSetC); + + //exclude_.addPair(a, b); + //exclude_.addPair(a, c); + //exclude_.addPair(b, c); + } -void SimInfo::refreshSim(){ + 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(); + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); + std::set rigidSetD = getRigidSet(d, atomGroups); - simtype fInfo; - int isError; - int n_global; - int* excl; + exclude_.addPairs(rigidSetA, rigidSetB); + exclude_.addPairs(rigidSetA, rigidSetC); + exclude_.addPairs(rigidSetA, rigidSetD); + exclude_.addPairs(rigidSetB, rigidSetC); + exclude_.addPairs(rigidSetB, rigidSetD); + exclude_.addPairs(rigidSetC, rigidSetD); - fInfo.dielect = 0.0; + /* + exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); + exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.addPairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.addPairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.addPairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); + + + exclude_.addPair(a, b); + exclude_.addPair(a, c); + exclude_.addPair(a, d); + exclude_.addPair(b, c); + exclude_.addPair(b, d); + exclude_.addPair(c, d); + */ + } - if( useDipoles ){ - if( useReactionField )fInfo.dielect = dielectric; + 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); + } + } + } + } - 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::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; - 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); + std::map > atomGroups; - if( isError ){ + Molecule::RigidBodyIterator rbIter; + RigidBody* rb; + Molecule::IntegrableObjectIterator ii; + StuntDouble* integrableObject; - sprintf( painCave.errMsg, - "There was an error setting the simulation information in fortran.\n" ); - painCave.isFatal = 1; - painCave.severity = OOPSE_ERROR; - simError(); + for (integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; + integrableObject = mol->nextIntegrableObject(ii)) { + + if (integrableObject->isRigidBody()) { + rb = static_cast(integrableObject); + std::vector atoms = rb->getAtoms(); + std::set rigidAtoms; + for (int i = 0; i < atoms.size(); ++i) { + rigidAtoms.insert(atoms[i]->getGlobalIndex()); + } + for (int i = 0; i < atoms.size(); ++i) { + atomGroups.insert(std::map >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms)); + } + } else { + std::set oneAtomSet; + oneAtomSet.insert(integrableObject->getGlobalIndex()); + atomGroups.insert(std::map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + } + } + + + for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { + a = bond->getAtomA()->getGlobalIndex(); + b = bond->getAtomB()->getGlobalIndex(); + exclude_.removePair(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(); + + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); + + exclude_.removePairs(rigidSetA, rigidSetB); + exclude_.removePairs(rigidSetA, rigidSetC); + exclude_.removePairs(rigidSetB, rigidSetC); + + //exclude_.removePair(a, b); + //exclude_.removePair(a, c); + //exclude_.removePair(b, c); + } + + 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(); + + std::set rigidSetA = getRigidSet(a, atomGroups); + std::set rigidSetB = getRigidSet(b, atomGroups); + std::set rigidSetC = getRigidSet(c, atomGroups); + std::set rigidSetD = getRigidSet(d, atomGroups); + + exclude_.removePairs(rigidSetA, rigidSetB); + exclude_.removePairs(rigidSetA, rigidSetC); + exclude_.removePairs(rigidSetA, rigidSetD); + exclude_.removePairs(rigidSetB, rigidSetC); + exclude_.removePairs(rigidSetB, rigidSetD); + exclude_.removePairs(rigidSetC, rigidSetD); + + /* + exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetB.begin(), rigidSetB.end()); + exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.removePairs(rigidSetA.begin(), rigidSetA.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetC.begin(), rigidSetC.end()); + exclude_.removePairs(rigidSetB.begin(), rigidSetB.end(), rigidSetD.begin(), rigidSetD.end()); + exclude_.removePairs(rigidSetC.begin(), rigidSetC.end(), rigidSetD.begin(), rigidSetD.end()); + + + exclude_.removePair(a, b); + exclude_.removePair(a, c); + exclude_.removePair(a, d); + exclude_.removePair(b, c); + exclude_.removePair(b, d); + exclude_.removePair(c, d); + */ + } + + 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); + } + } + } + } - + + + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { + int curStampId; + + //index from 0 + curStampId = moleculeStamps_.size(); + + moleculeStamps_.push_back(molStamp); + molStampIds_.insert(molStampIds_.end(), nmol, curStampId); + } + + void SimInfo::update() { + + setupSimType(); + #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(); -} + setupFortranParallel(); +#endif -void SimInfo::setDefaultRcut( double theRcut ){ + setupFortranSim(); + + //setup fortran force field + /** @deprecate */ + int isError = 0; + + setupElectrostaticSummationMethod( isError ); + setupSwitchingFunction(); + + if(isError){ + sprintf( painCave.errMsg, + "ForceField error: There was an error initializing the forceField in fortran.\n" ); + painCave.isFatal = 1; + simError(); + } - haveRcut = 1; - rCut = theRcut; - rList = rCut + 1.0; - - notifyFortranCutOffs( &rCut, &rSw, &rList ); -} + + setupCutoff(); -void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ + calcNdf(); + calcNdfRaw(); + calcNdfTrans(); - rSw = theRsw; - setDefaultRcut( theRcut ); -} + fortranInitialized_ = true; + } + std::set SimInfo::getUniqueAtomTypes() { + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::AtomIterator ai; + Atom* atom; + std::set atomTypes; -void SimInfo::checkCutOffs( void ){ - - if( boxIsInit ){ + 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(); - //we need to check cutOffs against the box + int useLennardJones = 0; + int useElectrostatic = 0; + int useEAM = 0; + int useSC = 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_->getUsePeriodicBoundaryConditions(); + int useRF; + int useSF; + std::string myMethod; + + // set the useRF logical + useRF = 0; + useSF = 0; + + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "REACTION_FIELD") { + useRF=1; + } else { + if (myMethod == "SHIFTED_FORCE") { + useSF = 1; + } + } + } + + //loop over all of the atom types + for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { + useLennardJones |= (*i)->isLennardJones(); + useElectrostatic |= (*i)->isElectrostatic(); + useEAM |= (*i)->isEAM(); + useSC |= (*i)->isSC(); + useCharge |= (*i)->isCharge(); + useDirectional |= (*i)->isDirectional(); + useDipole |= (*i)->isDipole(); + useGayBerne |= (*i)->isGayBerne(); + useSticky |= (*i)->isSticky(); + useStickyPower |= (*i)->isStickyPower(); + useShape |= (*i)->isShape(); + } + + if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { + useDirectionalAtom = 1; + } + + if (useCharge || useDipole) { + useElectrostatics = 1; + } + +#ifdef IS_MPI + int temp; + + temp = usePBC; + MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useDirectionalAtom; + MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useLennardJones; + MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useElectrostatics; + MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + 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); + + 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); - 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(); - } - -} + temp = useGayBerne; + MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); -void SimInfo::addProperty(GenericData* prop){ + temp = useEAM; + MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - 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()){ + temp = useSC; + MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - delete (*result).second; - (*result).second = prop; + temp = useShape; + MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useFLARB; + MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useRF; + MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = useSF; + MPI_Allreduce(&temp, &useSF, 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_StickyPower = useStickyPower; + fInfo_.SIM_uses_GayBerne = useGayBerne; + fInfo_.SIM_uses_EAM = useEAM; + fInfo_.SIM_uses_SC = useSC; + fInfo_.SIM_uses_Shapes = useShape; + fInfo_.SIM_uses_FLARB = useFLARB; + fInfo_.SIM_uses_RF = useRF; + fInfo_.SIM_uses_SF = useSF; + + if( myMethod == "REACTION_FIELD") { + 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{ - properties[prop->getID()] = prop; + void SimInfo::setupFortranSim() { + int isError; + int nExclude; + std::vector fortranGlobalGroupMembership; + + nExclude = exclude_.getSize(); + isError = 0; - } + //globalGroupMembership_ is filled by SimCreator + for (int i = 0; i < nGlobalAtoms_; i++) { + fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); + } + + //calculate mass ratio of cutoff group + std::vector mfact; + SimInfo::MoleculeIterator mi; + Molecule* mol; + Molecule::CutoffGroupIterator ci; + CutoffGroup* cg; + Molecule::AtomIterator ai; + Atom* atom; + RealType totalMass; + + //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)) { -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) + mfact.push_back(atom->getMass()/totalMass); + else + mfact.push_back( 1.0 ); + } + } + } -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(); - + //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) + std::vector identArray; - // Fix the silly fortran indexing problem + //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()); + } + } + + //fill molMembershipArray + //molMembershipArray is filled by SimCreator + std::vector molMembershipArray(nGlobalAtoms_); + for (int i = 0; i < nGlobalAtoms_; i++) { + molMembershipArray[i] = globalMolMembership_[i] + 1; + } + + //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, + "There was an error setting the simulation information in fortran.\n" ); + painCave.isFatal = 1; + painCave.severity = OOPSE_ERROR; + simError(); + } + #ifdef IS_MPI - numAtom = mpiSim->getNAtomsGlobal(); -#else - numAtom = n_atoms; + sprintf( checkPointMsg, + "succesfully sent the simulation information to fortran.\n"); + MPIcheckPoint(); +#endif // is_mpi + } + + +#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; + + 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; + } + + //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(¶llelData, &(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(); + + + } + #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)){ + void SimInfo::setupCutoff() { + + ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions(); - totalMass = myCutoffGroup->getMass(); - - for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); - cutoffAtom != NULL; - cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ - mfact.push_back(cutoffAtom->getMass()/totalMass); - } + // Check the cutoff policy + int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default + + std::string myPolicy; + if (forceFieldOptions_.haveCutoffPolicy()){ + myPolicy = forceFieldOptions_.getCutoffPolicy(); + }else if (simParams_->haveCutoffPolicy()) { + myPolicy = simParams_->getCutoffPolicy(); } + + if (!myPolicy.empty()){ + toUpper(myPolicy); + if (myPolicy == "MIX") { + cp = MIX_CUTOFF_POLICY; + } else { + if (myPolicy == "MAX") { + cp = MAX_CUTOFF_POLICY; + } else { + if (myPolicy == "TRADITIONAL") { + cp = TRADITIONAL_CUTOFF_POLICY; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + } + notifyFortranCutoffPolicy(&cp); + + // Check the Skin Thickness for neighborlists + RealType skin; + if (simParams_->haveSkinThickness()) { + skin = simParams_->getSkinThickness(); + notifyFortranSkinThickness(&skin); + } + + // Check if the cutoff was set explicitly: + if (simParams_->haveCutoffRadius()) { + rcut_ = simParams_->getCutoffRadius(); + if (simParams_->haveSwitchingRadius()) { + rsw_ = simParams_->getSwitchingRadius(); + } else { + if (fInfo_.SIM_uses_Charges | + fInfo_.SIM_uses_Dipoles | + fInfo_.SIM_uses_RF) { + + rsw_ = 0.85 * rcut_; + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the switchingRadius.\n" + "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n" + "\tswitchingRadius = %f. for this simulation\n", rsw_); + painCave.isFatal = 0; + simError(); + } else { + rsw_ = rcut_; + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for the switchingRadius.\n" + "\tOOPSE will use the same value as the cutoffRadius.\n" + "\tswitchingRadius = %f. for this simulation\n", rsw_); + painCave.isFatal = 0; + simError(); + } + } + + notifyFortranCutoffs(&rcut_, &rsw_); + + } else { + + // For electrostatic atoms, we'll assume a large safe value: + if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { + 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; + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") { + if (simParams_->haveSwitchingRadius()){ + sprintf(painCave.errMsg, + "SimInfo Warning: A value was set for the switchingRadius\n" + "\teven though the electrostaticSummationMethod was\n" + "\tset to %s\n", myMethod.c_str()); + painCave.isFatal = 1; + simError(); + } + } + } + + if (simParams_->haveSwitchingRadius()){ + rsw_ = simParams_->getSwitchingRadius(); + } else { + sprintf(painCave.errMsg, + "SimCreator Warning: No value was set for switchingRadius.\n" + "\tOOPSE will use a default value of\n" + "\t0.85 * cutoffRadius for the switchingRadius\n"); + painCave.isFatal = 0; + simError(); + rsw_ = 0.85 * rcut_; + } + notifyFortranCutoffs(&rcut_, &rsw_); + } else { + // We didn't set rcut explicitly, and we don't have electrostatic atoms, so + // We'll punt and let fortran figure out the cutoffs later. + + notifyFortranYouAreOnYourOwn(); + + } + } } -} + void SimInfo::setupElectrostaticSummationMethod( int isError ) { + + int errorOut; + int esm = NONE; + int sm = UNDAMPED; + RealType alphaVal; + RealType dielectric; + + errorOut = isError; + alphaVal = simParams_->getDampingAlpha(); + dielectric = simParams_->getDielectric(); + + if (simParams_->haveElectrostaticSummationMethod()) { + std::string myMethod = simParams_->getElectrostaticSummationMethod(); + toUpper(myMethod); + if (myMethod == "NONE") { + esm = NONE; + } else { + if (myMethod == "SWITCHING_FUNCTION") { + esm = SWITCHING_FUNCTION; + } else { + if (myMethod == "SHIFTED_POTENTIAL") { + esm = SHIFTED_POTENTIAL; + } else { + if (myMethod == "SHIFTED_FORCE") { + esm = SHIFTED_FORCE; + } else { + if (myMethod == "REACTION_FIELD") { + esm = REACTION_FIELD; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticSummationMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticSummationMethod must be one of: \"none\",\n" + "\t\"shifted_potential\", \"shifted_force\", or \n" + "\t\"reaction_field\".\n", myMethod.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + } + } + } + + if (simParams_->haveElectrostaticScreeningMethod()) { + std::string myScreen = simParams_->getElectrostaticScreeningMethod(); + toUpper(myScreen); + if (myScreen == "UNDAMPED") { + sm = UNDAMPED; + } else { + if (myScreen == "DAMPED") { + sm = DAMPED; + if (!simParams_->haveDampingAlpha()) { + //throw error + sprintf( painCave.errMsg, + "SimInfo warning: dampingAlpha was not specified in the input file.\n" + "\tA default value of %f (1/ang) will be used.\n", alphaVal); + painCave.isFatal = 0; + simError(); + } + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown electrostaticScreeningMethod.\n" + "\t(Input file specified %s .)\n" + "\telectrostaticScreeningMethod must be one of: \"undamped\"\n" + "or \"damped\".\n", myScreen.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + + // let's pass some summation method variables to fortran + setElectrostaticSummationMethod( &esm ); + setFortranElectrostaticMethod( &esm ); + setScreeningMethod( &sm ); + setDampingAlpha( &alphaVal ); + setReactionFieldDielectric( &dielectric ); + initFortranFF( &errorOut ); + } + + void SimInfo::setupSwitchingFunction() { + int ft = CUBIC; + + if (simParams_->haveSwitchingFunctionType()) { + std::string funcType = simParams_->getSwitchingFunctionType(); + toUpper(funcType); + if (funcType == "CUBIC") { + ft = CUBIC; + } else { + if (funcType == "FIFTH_ORDER_POLYNOMIAL") { + ft = FIFTH_ORDER_POLY; + } else { + // throw error + sprintf( painCave.errMsg, + "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); + painCave.isFatal = 1; + simError(); + } + } + } + + // send switching function notification to switcheroo + setFunctionType(&ft); + + } + + 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); + 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 + 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 + + comVel /= totalMass; + + 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; + + } + + 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; + + + 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; + } + + +}//end namespace oopse +