--- branches/development/src/brains/SimInfo.cpp 2011/08/04 20:04:35 1601 +++ branches/development/src/brains/SimInfo.cpp 2012/07/03 18:32:27 1764 @@ -36,7 +36,8 @@ * [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] Vardeman & Gezelter, in progress (2009). + * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). + * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ /** @@ -58,8 +59,11 @@ #include "utils/simError.h" #include "selection/SelectionManager.hpp" #include "io/ForceFieldOptions.hpp" -#include "UseTheForce/ForceField.hpp" +#include "brains/ForceField.hpp" #include "nonbonded/SwitchingFunction.hpp" +#ifdef IS_MPI +#include +#endif using namespace std; namespace OpenMD { @@ -68,10 +72,10 @@ namespace OpenMD { 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), + nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0), nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0), nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), - nConstraints_(0), sman_(NULL), topologyDone_(false), + nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), calcBoxDipole_(false), useAtomicVirial_(true) { MoleculeStamp* molStamp; @@ -221,13 +225,17 @@ namespace OpenMD { void SimInfo::calcNdf() { - int ndf_local; + int ndf_local, nfq_local; MoleculeIterator i; vector::iterator j; + vector::iterator k; + Molecule* mol; StuntDouble* integrableObject; + Atom* atom; ndf_local = 0; + nfq_local = 0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; @@ -242,17 +250,26 @@ namespace OpenMD { ndf_local += 3; } } - + } + for (atom = mol->beginFluctuatingCharge(k); atom != NULL; + atom = mol->nextFluctuatingCharge(k)) { + if (atom->isFluctuatingCharge()) { + nfq_local++; + } } } + ndfLocal_ = ndf_local; + // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; #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 // nZconstraints_ is global, as are the 3 COM translations for the @@ -776,13 +793,15 @@ namespace OpenMD { 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; @@ -793,11 +812,15 @@ namespace OpenMD { temp = usesElectrostatic; MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + + temp = usesFluctuatingCharges; + MPI_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); #else usesDirectionalAtoms_ = usesDirectional; usesMetallicAtoms_ = usesMetallic; usesElectrostaticAtoms_ = usesElectrostatic; + usesFluctuatingCharges_ = usesFluctuatingCharges; #endif @@ -965,205 +988,14 @@ namespace OpenMD { } } - - 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; - - } - 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; - - - 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); } @@ -1171,44 +1003,6 @@ namespace OpenMD { 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,3.0/2.0)*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,3.0/2.0)*sqrt(detI); - return; - } /* void SimInfo::setStuntDoubleFromGlobalIndex(vector v) { assert( v.size() == nAtoms_ + nRigidBodies_);