--- branches/development/src/brains/SimInfo.cpp 2012/05/26 18:13:43 1725 +++ branches/development/src/brains/SimInfo.cpp 2012/07/09 14:15:52 1769 @@ -231,26 +231,28 @@ namespace OpenMD { 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()) { @@ -259,6 +261,8 @@ namespace OpenMD { } } + ndfLocal_ = ndf_local; + // n_constraints is local, so subtract them on each processor ndf_local -= nConstraints_; @@ -310,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; @@ -382,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) { @@ -400,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)); } } @@ -535,14 +539,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) { @@ -553,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)); } } @@ -788,10 +791,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 +802,24 @@ namespace OpenMD { 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_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; @@ -987,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); } @@ -1192,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,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_);