| 35 |  | * | 
| 36 |  | * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | 
| 37 |  | * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | 
| 38 | < | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). | 
| 38 | > | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). | 
| 39 |  | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 |  | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 |  | */ | 
| 47 |  | * @version 1.0 | 
| 48 |  | */ | 
| 49 |  |  | 
| 50 | + | #ifdef IS_MPI | 
| 51 | + | #include <mpi.h> | 
| 52 | + | #endif | 
| 53 |  | #include <algorithm> | 
| 54 |  | #include <set> | 
| 55 |  | #include <map> | 
| 64 |  | #include "io/ForceFieldOptions.hpp" | 
| 65 |  | #include "brains/ForceField.hpp" | 
| 66 |  | #include "nonbonded/SwitchingFunction.hpp" | 
| 64 | – | #ifdef IS_MPI | 
| 65 | – | #include <mpi.h> | 
| 66 | – | #endif | 
| 67 |  |  | 
| 68 |  | using namespace std; | 
| 69 |  | namespace OpenMD { | 
| 72 |  | forceField_(ff), simParams_(simParams), | 
| 73 |  | ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), | 
| 74 |  | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), | 
| 75 | < | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), nGlobalFluctuatingCharges_(0), | 
| 76 | < | nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0), | 
| 77 | < | nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0), | 
| 78 | < | nConstraints_(0), nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), | 
| 75 | > | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 76 | > | nGlobalFluctuatingCharges_(0), nGlobalBonds_(0), nGlobalBends_(0), | 
| 77 | > | nGlobalTorsions_(0), nGlobalInversions_(0), nAtoms_(0), nBonds_(0), | 
| 78 | > | nBends_(0), nTorsions_(0), nInversions_(0), nRigidBodies_(0), | 
| 79 | > | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), | 
| 80 | > | nFluctuatingCharges_(0), sman_(NULL), topologyDone_(false), | 
| 81 |  | calcBoxDipole_(false), useAtomicVirial_(true) { | 
| 82 |  |  | 
| 83 |  | MoleculeStamp* molStamp; | 
| 93 |  | for (vector<Component*>::iterator i = components.begin(); | 
| 94 |  | i !=components.end(); ++i) { | 
| 95 |  | molStamp = (*i)->getMoleculeStamp(); | 
| 96 | + | if ( (*i)->haveRegion() ) { | 
| 97 | + | molStamp->setRegion( (*i)->getRegion() ); | 
| 98 | + | } else { | 
| 99 | + | // set the region to a disallowed value: | 
| 100 | + | molStamp->setRegion( -1 ); | 
| 101 | + | } | 
| 102 | + |  | 
| 103 |  | nMolWithSameStamp = (*i)->getNMol(); | 
| 104 |  |  | 
| 105 |  | addMoleculeStamp(molStamp, nMolWithSameStamp); | 
| 106 |  |  | 
| 107 |  | //calculate atoms in molecules | 
| 108 | < | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; | 
| 108 | > | nGlobalAtoms_ += molStamp->getNAtoms() * nMolWithSameStamp; | 
| 109 | > | nGlobalBonds_ += molStamp->getNBonds() * nMolWithSameStamp; | 
| 110 | > | nGlobalBends_ += molStamp->getNBends() * nMolWithSameStamp; | 
| 111 | > | nGlobalTorsions_ += molStamp->getNTorsions() * nMolWithSameStamp; | 
| 112 | > | nGlobalInversions_ += molStamp->getNInversions() * nMolWithSameStamp; | 
| 113 |  |  | 
| 114 |  | //calculate atoms in cutoff groups | 
| 115 |  | int nAtomsInGroups = 0; | 
| 281 |  | ndf_local -= nConstraints_; | 
| 282 |  |  | 
| 283 |  | #ifdef IS_MPI | 
| 284 | < | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 285 | < | MPI_Allreduce(&nfq_local,&nGlobalFluctuatingCharges_,1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 284 | > | MPI_Allreduce(&ndf_local, &ndf_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 285 | > | MPI_Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1, | 
| 286 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 287 | > | // MPI::COMM_WORLD.Allreduce(&ndf_local, &ndf_, 1, MPI::INT,MPI::SUM); | 
| 288 | > | // MPI::COMM_WORLD.Allreduce(&nfq_local, &nGlobalFluctuatingCharges_, 1, | 
| 289 | > | //                           MPI::INT, MPI::SUM); | 
| 290 |  | #else | 
| 291 |  | ndf_ = ndf_local; | 
| 292 |  | nGlobalFluctuatingCharges_ = nfq_local; | 
| 300 |  |  | 
| 301 |  | int SimInfo::getFdf() { | 
| 302 |  | #ifdef IS_MPI | 
| 303 | < | MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 303 | > | MPI_Allreduce(&fdf_local, &fdf_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 304 | > | // MPI::COMM_WORLD.Allreduce(&fdf_local, &fdf_, 1, MPI::INT, MPI::SUM); | 
| 305 |  | #else | 
| 306 |  | fdf_ = fdf_local; | 
| 307 |  | #endif | 
| 357 |  | } | 
| 358 |  |  | 
| 359 |  | #ifdef IS_MPI | 
| 360 | < | MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 360 | > | MPI_Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 361 | > | // MPI::COMM_WORLD.Allreduce(&ndfRaw_local, &ndfRaw_, 1, MPI::INT, MPI::SUM); | 
| 362 |  | #else | 
| 363 |  | ndfRaw_ = ndfRaw_local; | 
| 364 |  | #endif | 
| 371 |  |  | 
| 372 |  |  | 
| 373 |  | #ifdef IS_MPI | 
| 374 | < | MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); | 
| 374 | > | MPI_Allreduce(&ndfTrans_local, &ndfTrans_, 1, | 
| 375 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 376 | > | // MPI::COMM_WORLD.Allreduce(&ndfTrans_local, &ndfTrans_, 1, | 
| 377 | > | //                           MPI::INT, MPI::SUM); | 
| 378 |  | #else | 
| 379 |  | ndfTrans_ = ndfTrans_local; | 
| 380 |  | #endif | 
| 431 |  | atomGroups.insert(map<int, set<int> >::value_type(sd->getGlobalIndex(), oneAtomSet)); | 
| 432 |  | } | 
| 433 |  | } | 
| 434 | + |  | 
| 435 |  |  | 
| 436 |  | for (bond= mol->beginBond(bondIter); bond != NULL; | 
| 437 |  | bond = mol->nextBond(bondIter)) { | 
| 756 |  | // count_local holds the number of found types on this processor | 
| 757 |  | int count_local = foundTypes.size(); | 
| 758 |  |  | 
| 759 | < | int nproc = MPI::COMM_WORLD.Get_size(); | 
| 759 | > | int nproc; | 
| 760 | > | MPI_Comm_size( MPI_COMM_WORLD, &nproc); | 
| 761 | > | // int nproc = MPI::COMM_WORLD.Get_size(); | 
| 762 |  |  | 
| 763 |  | // we need arrays to hold the counts and displacement vectors for | 
| 764 |  | // all processors | 
| 766 |  | vector<int> disps(nproc, 0); | 
| 767 |  |  | 
| 768 |  | // fill the counts array | 
| 769 | < | MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], | 
| 770 | < | 1, MPI::INT); | 
| 769 | > | MPI_Allgather(&count_local, 1, MPI_INT, &counts[0], | 
| 770 | > | 1, MPI_INT, MPI_COMM_WORLD); | 
| 771 | > | // MPI::COMM_WORLD.Allgather(&count_local, 1, MPI::INT, &counts[0], | 
| 772 | > | //                           1, MPI::INT); | 
| 773 |  |  | 
| 774 |  | // use the processor counts to compute the displacement array | 
| 775 |  | disps[0] = 0; | 
| 783 |  | vector<int> ftGlobal(totalCount); | 
| 784 |  |  | 
| 785 |  | // now spray out the foundTypes to all the other processors: | 
| 786 | < | MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, | 
| 787 | < | &ftGlobal[0], &counts[0], &disps[0], | 
| 788 | < | MPI::INT); | 
| 786 | > | MPI_Allgatherv(&foundTypes[0], count_local, MPI_INT, | 
| 787 | > | &ftGlobal[0], &counts[0], &disps[0], | 
| 788 | > | MPI_INT, MPI_COMM_WORLD); | 
| 789 | > | // MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT, | 
| 790 | > | //                            &ftGlobal[0], &counts[0], &disps[0], | 
| 791 | > | //                            MPI::INT); | 
| 792 |  |  | 
| 793 |  | vector<int>::iterator j; | 
| 794 |  |  | 
| 810 |  | return atomTypes; | 
| 811 |  | } | 
| 812 |  |  | 
| 813 | + |  | 
| 814 | + | int getGlobalCountOfType(AtomType* atype) { | 
| 815 | + | /* | 
| 816 | + | set<AtomType*> atypes = getSimulatedAtomTypes(); | 
| 817 | + | map<AtomType*, int> counts_; | 
| 818 | + |  | 
| 819 | + | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 820 | + | for(atom = mol->beginAtom(ai); atom != NULL; | 
| 821 | + | atom = mol->nextAtom(ai)) { | 
| 822 | + | atom->getAtomType(); | 
| 823 | + | } | 
| 824 | + | } | 
| 825 | + | */ | 
| 826 | + | return 0; | 
| 827 | + | } | 
| 828 | + |  | 
| 829 |  | void SimInfo::setupSimVariables() { | 
| 830 |  | useAtomicVirial_ = simParams_->getUseAtomicVirial(); | 
| 831 |  | // we only call setAccumulateBoxDipole if the accumulateBoxDipole | 
| 852 |  | } | 
| 853 |  |  | 
| 854 |  | #ifdef IS_MPI | 
| 855 | < | bool temp; | 
| 855 | > | int temp; | 
| 856 | > |  | 
| 857 |  | temp = usesDirectional; | 
| 858 | < | MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, | 
| 859 | < | MPI::LOR); | 
| 860 | < |  | 
| 858 | > | MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD); | 
| 859 | > | usesDirectionalAtoms_ = (temp == 0) ? false : true; | 
| 860 | > |  | 
| 861 | > | // MPI::COMM_WORLD.Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI::BOOL, | 
| 862 | > | //                           MPI::LOR); | 
| 863 | > |  | 
| 864 |  | temp = usesMetallic; | 
| 865 | < | MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, | 
| 866 | < | MPI::LOR); | 
| 865 | > | MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD); | 
| 866 | > | usesMetallicAtoms_ = (temp == 0) ? false : true; | 
| 867 | > |  | 
| 868 | > | // MPI::COMM_WORLD.Allreduce(&temp, &usesMetallicAtoms_, 1, MPI::BOOL, | 
| 869 | > | //                           MPI::LOR); | 
| 870 |  |  | 
| 871 |  | temp = usesElectrostatic; | 
| 872 | < | MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, | 
| 873 | < | MPI::LOR); | 
| 872 | > | MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD); | 
| 873 | > | usesElectrostaticAtoms_ = (temp == 0) ? false : true; | 
| 874 |  |  | 
| 875 | + | // MPI::COMM_WORLD.Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI::BOOL, | 
| 876 | + | //                           MPI::LOR); | 
| 877 | + |  | 
| 878 |  | temp = usesFluctuatingCharges; | 
| 879 | < | MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, | 
| 880 | < | MPI::LOR); | 
| 879 | > | MPI_Allreduce(MPI_IN_PLACE, &temp, 1, MPI_INT,  MPI_LOR, MPI_COMM_WORLD); | 
| 880 | > | usesFluctuatingCharges_ = (temp == 0) ? false : true; | 
| 881 | > |  | 
| 882 | > | // MPI::COMM_WORLD.Allreduce(&temp, &usesFluctuatingCharges_, 1, MPI::BOOL, | 
| 883 | > | //                           MPI::LOR); | 
| 884 | > |  | 
| 885 |  | #else | 
| 886 |  |  | 
| 887 |  | usesDirectionalAtoms_ = usesDirectional; | 
| 973 |  | } | 
| 974 |  | } | 
| 975 |  |  | 
| 976 | < | // Build the identArray_ | 
| 976 | > | // Build the identArray_ and regions_ | 
| 977 |  |  | 
| 978 |  | identArray_.clear(); | 
| 979 | < | identArray_.reserve(getNAtoms()); | 
| 980 | < | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 979 | > | identArray_.reserve(getNAtoms()); | 
| 980 | > | regions_.clear(); | 
| 981 | > | regions_.reserve(getNAtoms()); | 
| 982 | > |  | 
| 983 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 984 | > | int reg = mol->getRegion(); | 
| 985 |  | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 986 |  | identArray_.push_back(atom->getIdent()); | 
| 987 | + | regions_.push_back(reg); | 
| 988 |  | } | 
| 989 |  | } | 
| 990 | < |  | 
| 926 | < | //scan topology | 
| 927 | < |  | 
| 928 | < | int* excludeList = excludedInteractions_.getPairList(); | 
| 929 | < | int* oneTwoList = oneTwoInteractions_.getPairList(); | 
| 930 | < | int* oneThreeList = oneThreeInteractions_.getPairList(); | 
| 931 | < | int* oneFourList = oneFourInteractions_.getPairList(); | 
| 932 | < |  | 
| 990 | > |  | 
| 991 |  | topologyDone_ = true; | 
| 992 |  | } | 
| 993 |  |  | 
| 1022 |  | delete sman_; | 
| 1023 |  | sman_ = sman; | 
| 1024 |  |  | 
| 967 | – | Molecule* mol; | 
| 968 | – | RigidBody* rb; | 
| 969 | – | Atom* atom; | 
| 970 | – | CutoffGroup* cg; | 
| 1025 |  | SimInfo::MoleculeIterator mi; | 
| 1026 | + | Molecule::AtomIterator ai; | 
| 1027 |  | Molecule::RigidBodyIterator rbIter; | 
| 973 | – | Molecule::AtomIterator atomIter; | 
| 1028 |  | Molecule::CutoffGroupIterator cgIter; | 
| 1029 | + | Molecule::BondIterator bondIter; | 
| 1030 | + | Molecule::BendIterator bendIter; | 
| 1031 | + | Molecule::TorsionIterator torsionIter; | 
| 1032 | + | Molecule::InversionIterator inversionIter; | 
| 1033 |  |  | 
| 1034 | + | Molecule* mol; | 
| 1035 | + | Atom* atom; | 
| 1036 | + | RigidBody* rb; | 
| 1037 | + | CutoffGroup* cg; | 
| 1038 | + | Bond* bond; | 
| 1039 | + | Bend* bend; | 
| 1040 | + | Torsion* torsion; | 
| 1041 | + | Inversion* inversion; | 
| 1042 | + |  | 
| 1043 |  | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 1044 |  |  | 
| 1045 | < | for (atom = mol->beginAtom(atomIter); atom != NULL; | 
| 1046 | < | atom = mol->nextAtom(atomIter)) { | 
| 1045 | > | for (atom = mol->beginAtom(ai); atom != NULL; | 
| 1046 | > | atom = mol->nextAtom(ai)) { | 
| 1047 |  | atom->setSnapshotManager(sman_); | 
| 1048 | < | } | 
| 982 | < |  | 
| 1048 | > | } | 
| 1049 |  | for (rb = mol->beginRigidBody(rbIter); rb != NULL; | 
| 1050 |  | rb = mol->nextRigidBody(rbIter)) { | 
| 1051 |  | rb->setSnapshotManager(sman_); | 
| 1052 |  | } | 
| 987 | – |  | 
| 1053 |  | for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; | 
| 1054 |  | cg = mol->nextCutoffGroup(cgIter)) { | 
| 1055 |  | cg->setSnapshotManager(sman_); | 
| 1056 |  | } | 
| 1057 | < | } | 
| 1058 | < |  | 
| 1057 | > | for (bond = mol->beginBond(bondIter); bond != NULL; | 
| 1058 | > | bond = mol->nextBond(bondIter)) { | 
| 1059 | > | bond->setSnapshotManager(sman_); | 
| 1060 | > | } | 
| 1061 | > | for (bend = mol->beginBend(bendIter); bend != NULL; | 
| 1062 | > | bend = mol->nextBend(bendIter)) { | 
| 1063 | > | bend->setSnapshotManager(sman_); | 
| 1064 | > | } | 
| 1065 | > | for (torsion = mol->beginTorsion(torsionIter); torsion != NULL; | 
| 1066 | > | torsion = mol->nextTorsion(torsionIter)) { | 
| 1067 | > | torsion->setSnapshotManager(sman_); | 
| 1068 | > | } | 
| 1069 | > | for (inversion = mol->beginInversion(inversionIter); inversion != NULL; | 
| 1070 | > | inversion = mol->nextInversion(inversionIter)) { | 
| 1071 | > | inversion->setSnapshotManager(sman_); | 
| 1072 | > | } | 
| 1073 | > | } | 
| 1074 |  | } | 
| 1075 |  |  | 
| 1076 |  |  | 
| 1081 |  |  | 
| 1082 |  |  | 
| 1083 |  | StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) { | 
| 1084 | < | if (index >= IOIndexToIntegrableObject.size()) { | 
| 1084 | > | if (index >= int(IOIndexToIntegrableObject.size())) { | 
| 1085 |  | sprintf(painCave.errMsg, | 
| 1086 |  | "SimInfo::getIOIndexToIntegrableObject Error: Integrable Object\n" | 
| 1087 |  | "\tindex exceeds number of known objects!\n"); | 
| 1099 |  | int SimInfo::getNGlobalConstraints() { | 
| 1100 |  | int nGlobalConstraints; | 
| 1101 |  | #ifdef IS_MPI | 
| 1102 | < | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, | 
| 1103 | < | MPI_COMM_WORLD); | 
| 1102 | > | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, | 
| 1103 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 1104 | > | // MPI::COMM_WORLD.Allreduce(&nConstraints_, &nGlobalConstraints, 1, | 
| 1105 | > | //                           MPI::INT, MPI::SUM); | 
| 1106 |  | #else | 
| 1107 |  | nGlobalConstraints =  nConstraints_; | 
| 1108 |  | #endif |