--- branches/development/src/brains/SimInfo.cpp 2012/05/26 18:13:43 1725 +++ branches/development/src/brains/SimInfo.cpp 2013/01/09 19:27:52 1825 @@ -88,7 +88,8 @@ namespace OpenMD { vector components = simParams->getComponents(); - for (vector::iterator i = components.begin(); i !=components.end(); ++i) { + for (vector::iterator i = components.begin(); + i !=components.end(); ++i) { molStamp = (*i)->getMoleculeStamp(); nMolWithSameStamp = (*i)->getNMol(); @@ -231,26 +232,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,12 +262,15 @@ namespace OpenMD { } } + 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); + MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM); + MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1, + MPI::INT, MPI::SUM); #else ndf_ = ndf_local; nGlobalFluctuatingCharges_ = nfq_local; @@ -278,7 +284,7 @@ namespace OpenMD { int SimInfo::getFdf() { #ifdef IS_MPI - MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM); #else fdf_ = fdf_local; #endif @@ -310,19 +316,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; @@ -333,7 +340,7 @@ namespace OpenMD { } #ifdef IS_MPI - MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM); #else ndfRaw_ = ndfRaw_local; #endif @@ -346,7 +353,8 @@ namespace OpenMD { #ifdef IS_MPI - MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1, + MPI::INT, MPI::SUM); #else ndfTrans_ = ndfTrans_local; #endif @@ -382,14 +390,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 +407,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 +542,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 +559,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)); } } @@ -778,7 +784,8 @@ namespace OpenMD { void SimInfo::setupSimVariables() { useAtomicVirial_ = simParams_->getUseAtomicVirial(); - // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true + // we only call setAccumulateBoxDipole if the accumulateBoxDipole + // parameter is true calcBoxDipole_ = false; if ( simParams_->haveAccumulateBoxDipole() ) if ( simParams_->getAccumulateBoxDipole() ) { @@ -788,10 +795,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 +806,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; @@ -868,7 +879,6 @@ namespace OpenMD { void SimInfo::prepareTopology() { - int nExclude, nOneTwo, nOneThree, nOneFour; //calculate mass ratio of cutoff group SimInfo::MoleculeIterator mi; @@ -917,11 +927,6 @@ namespace OpenMD { //scan topology - nExclude = excludedInteractions_.getSize(); - nOneTwo = oneTwoInteractions_.getSize(); - nOneThree = oneThreeInteractions_.getSize(); - nOneFour = oneFourInteractions_.getSize(); - int* excludeList = excludedInteractions_.getPairList(); int* oneTwoList = oneTwoInteractions_.getPairList(); int* oneThreeList = oneThreeInteractions_.getPairList(); @@ -972,280 +977,52 @@ namespace OpenMD { for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { - for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { + 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)) { + for (rb = mol->beginRigidBody(rbIter); rb != NULL; + rb = mol->nextRigidBody(rbIter)) { rb->setSnapshotManager(sman_); } - for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { + for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; + cg = mol->nextCutoffGroup(cgIter)) { cg->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; - - } - 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); + if (index >= IOIndexToIntegrableObject.size()) { + sprintf(painCave.errMsg, + "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n" + "\tindex exceeds number of known objects!\n"); + painCave.isFatal = 1; + simError(); + return NULL; + } else + return IOIndexToIntegrableObject.at(index); } 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_); - sdByGlobalIndex_ = v; - } - - StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) { - //assert(index < nAtoms_ + nRigidBodies_); - return sdByGlobalIndex_.at(index); - } -*/ int SimInfo::getNGlobalConstraints() { int nGlobalConstraints; #ifdef IS_MPI - MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, - MPI_COMM_WORLD); + MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1, + MPI::INT, MPI::SUM); #else nGlobalConstraints = nConstraints_; #endif