--- branches/development/src/brains/SimInfo.cpp 2011/09/13 22:05:04 1627 +++ branches/development/src/brains/SimInfo.cpp 2012/07/09 14:15:52 1769 @@ -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,7 +59,7 @@ #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 @@ -71,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; @@ -224,38 +225,53 @@ 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; + StuntDouble* sd; + Atom* atom; ndf_local = 0; + nfq_local = 0; for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { - for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(j)) { + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { + ndf_local += 3; - if (integrableObject->isDirectional()) { - if (integrableObject->isLinear()) { + if (sd->isDirectional()) { + if (sd->isLinear()) { ndf_local += 2; } else { ndf_local += 3; } } - } + + for (atom = mol->beginFluctuatingCharge(k); atom != NULL; + atom = mol->nextFluctuatingCharge(k)) { + if (atom->isFluctuatingCharge()) { + nfq_local++; + } + } } + ndfLocal_ = ndf_local; + // 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 @@ -298,19 +314,20 @@ namespace OpenMD { MoleculeIterator i; vector::iterator j; Molecule* mol; - StuntDouble* integrableObject; + StuntDouble* sd; // 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)) { + for (sd = mol->beginIntegrableObject(j); sd != NULL; + sd = mol->nextIntegrableObject(j)) { + ndfRaw_local += 3; - if (integrableObject->isDirectional()) { - if (integrableObject->isLinear()) { + if (sd->isDirectional()) { + if (sd->isLinear()) { ndfRaw_local += 2; } else { ndfRaw_local += 3; @@ -370,14 +387,13 @@ namespace OpenMD { Molecule::RigidBodyIterator rbIter; RigidBody* rb; Molecule::IntegrableObjectIterator ii; - StuntDouble* integrableObject; + StuntDouble* sd; - for (integrableObject = mol->beginIntegrableObject(ii); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { - if (integrableObject->isRigidBody()) { - rb = static_cast(integrableObject); + if (sd->isRigidBody()) { + rb = static_cast(sd); vector atoms = rb->getAtoms(); set rigidAtoms; for (int i = 0; i < static_cast(atoms.size()); ++i) { @@ -388,8 +404,8 @@ namespace OpenMD { } } else { set oneAtomSet; - oneAtomSet.insert(integrableObject->getGlobalIndex()); - atomGroups.insert(map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + oneAtomSet.insert(sd->getGlobalIndex()); + atomGroups.insert(map >::value_type(sd->getGlobalIndex(), oneAtomSet)); } } @@ -523,14 +539,13 @@ namespace OpenMD { Molecule::RigidBodyIterator rbIter; RigidBody* rb; Molecule::IntegrableObjectIterator ii; - StuntDouble* integrableObject; - - for (integrableObject = mol->beginIntegrableObject(ii); - integrableObject != NULL; - integrableObject = mol->nextIntegrableObject(ii)) { + StuntDouble* sd; + + for (sd = mol->beginIntegrableObject(ii); sd != NULL; + sd = mol->nextIntegrableObject(ii)) { - if (integrableObject->isRigidBody()) { - rb = static_cast(integrableObject); + if (sd->isRigidBody()) { + rb = static_cast(sd); vector atoms = rb->getAtoms(); set rigidAtoms; for (int i = 0; i < static_cast(atoms.size()); ++i) { @@ -541,8 +556,8 @@ namespace OpenMD { } } else { set oneAtomSet; - oneAtomSet.insert(integrableObject->getGlobalIndex()); - atomGroups.insert(map >::value_type(integrableObject->getGlobalIndex(), oneAtomSet)); + oneAtomSet.insert(sd->getGlobalIndex()); + atomGroups.insert(map >::value_type(sd->getGlobalIndex(), oneAtomSet)); } } @@ -776,31 +791,41 @@ namespace OpenMD { set::iterator i; set atomTypes; atomTypes = getSimulatedAtomTypes(); - int usesElectrostatic = 0; - int usesMetallic = 0; - int usesDirectional = 0; + bool usesElectrostatic = false; + bool usesMetallic = false; + bool usesDirectional = false; + bool usesFluctuatingCharges = false; //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; + +#ifdef IS_MPI + bool temp; temp = usesDirectional; - MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, + MPI::LOR); + temp = usesMetallic; - MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, + MPI::LOR); temp = usesElectrostatic; - MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, + MPI::LOR); + + temp = usesFluctuatingCharges; + MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, + MPI::LOR); #else usesDirectionalAtoms_ = usesDirectional; usesMetallicAtoms_ = usesMetallic; usesElectrostaticAtoms_ = usesElectrostatic; + usesFluctuatingCharges_ = usesFluctuatingCharges; #endif @@ -969,204 +994,13 @@ 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); } @@ -1174,44 +1008,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_);