--- branches/development/src/brains/SimInfo.cpp 2012/05/26 18:13:43 1725 +++ branches/development/src/brains/SimInfo.cpp 2012/07/06 22:01:58 1767 @@ -259,6 +259,8 @@ namespace OpenMD { } } + ndfLocal_ = ndf_local; + // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; @@ -788,10 +790,10 @@ namespace OpenMD { set::iterator i; set atomTypes; atomTypes = getSimulatedAtomTypes(); - int usesElectrostatic = 0; - int usesMetallic = 0; - int usesDirectional = 0; - int usesFluctuatingCharges = 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(); @@ -799,20 +801,24 @@ namespace OpenMD { usesDirectional |= (*i)->isDirectional(); usesFluctuatingCharges |= (*i)->isFluctuatingCharge(); } - -#ifdef IS_MPI - int temp; - temp = usesDirectional; - MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); - + +#ifdef IS_MPI + bool temp; + temp = usesDirectional; + 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_Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, + MPI::LOR); #else usesDirectionalAtoms_ = usesDirectional; @@ -986,205 +992,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); } @@ -1192,44 +1007,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,geomCnst)*sqrt(det); - return; - } - - void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){ - Mat3x3d intTensor; - Vector3d dummyAngMom; - RealType sysconstants; - RealType geomCnst; - - geomCnst = 3.0/2.0; - /* Get the inertial tensor and angular momentum for free*/ - getInertiaTensor(intTensor,dummyAngMom); - - detI = intTensor.determinant(); - sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_; - volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,geomCnst)*sqrt(detI); - return; - } /* void SimInfo::setStuntDoubleFromGlobalIndex(vector v) { assert( v.size() == nAtoms_ + nRigidBodies_);