| 54 | 
  | 
#include "math/Vector3.hpp" | 
| 55 | 
  | 
#include "primitives/Molecule.hpp" | 
| 56 | 
  | 
#include "primitives/StuntDouble.hpp" | 
| 57 | 
– | 
#include "UseTheForce/DarkSide/neighborLists_interface.h" | 
| 58 | 
– | 
#include "UseTheForce/doForces_interface.h" | 
| 57 | 
  | 
#include "utils/MemoryUtils.hpp" | 
| 58 | 
  | 
#include "utils/simError.h" | 
| 59 | 
  | 
#include "selection/SelectionManager.hpp" | 
| 61 | 
  | 
#include "UseTheForce/ForceField.hpp" | 
| 62 | 
  | 
#include "nonbonded/SwitchingFunction.hpp" | 
| 63 | 
  | 
 | 
| 66 | 
– | 
#ifdef IS_MPI | 
| 67 | 
– | 
#include "UseTheForce/mpiComponentPlan.h" | 
| 68 | 
– | 
#include "UseTheForce/DarkSide/simParallel_interface.h" | 
| 69 | 
– | 
#endif  | 
| 70 | 
– | 
 | 
| 64 | 
  | 
using namespace std; | 
| 65 | 
  | 
namespace OpenMD { | 
| 66 | 
  | 
   | 
| 71 | 
  | 
    nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), | 
| 72 | 
  | 
    nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nInversions_(0),  | 
| 73 | 
  | 
    nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),  | 
| 74 | 
< | 
    nConstraints_(0), sman_(NULL), fortranInitialized_(false),  | 
| 74 | 
> | 
    nConstraints_(0), sman_(NULL), topologyDone_(false),  | 
| 75 | 
  | 
    calcBoxDipole_(false), useAtomicVirial_(true) {     | 
| 76 | 
  | 
     | 
| 77 | 
  | 
    MoleculeStamp* molStamp; | 
| 125 | 
  | 
    //equal to the total number of atoms minus number of atoms belong to  | 
| 126 | 
  | 
    //cutoff group defined in meta-data file plus the number of cutoff  | 
| 127 | 
  | 
    //groups defined in meta-data file | 
| 135 | 
– | 
    std::cerr << "nGA = " << nGlobalAtoms_ << "\n"; | 
| 136 | 
– | 
    std::cerr << "nCA = " << nCutoffAtoms << "\n"; | 
| 137 | 
– | 
    std::cerr << "nG = " << nGroups << "\n"; | 
| 128 | 
  | 
 | 
| 129 | 
  | 
    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | 
| 140 | 
– | 
 | 
| 141 | 
– | 
    std::cerr << "nGCG = " << nGlobalCutoffGroups_ << "\n"; | 
| 130 | 
  | 
     | 
| 131 | 
  | 
    //every free atom (atom does not belong to rigid bodies) is an  | 
| 132 | 
  | 
    //integrable object therefore the total number of integrable objects  | 
| 268 | 
  | 
    fdf_ = fdf_local; | 
| 269 | 
  | 
#endif | 
| 270 | 
  | 
    return fdf_; | 
| 271 | 
+ | 
  } | 
| 272 | 
+ | 
   | 
| 273 | 
+ | 
  unsigned int SimInfo::getNLocalCutoffGroups(){ | 
| 274 | 
+ | 
    int nLocalCutoffAtoms = 0; | 
| 275 | 
+ | 
    Molecule* mol; | 
| 276 | 
+ | 
    MoleculeIterator mi; | 
| 277 | 
+ | 
    CutoffGroup* cg; | 
| 278 | 
+ | 
    Molecule::CutoffGroupIterator ci; | 
| 279 | 
+ | 
     | 
| 280 | 
+ | 
    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) { | 
| 281 | 
+ | 
       | 
| 282 | 
+ | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL;  | 
| 283 | 
+ | 
           cg = mol->nextCutoffGroup(ci)) { | 
| 284 | 
+ | 
        nLocalCutoffAtoms += cg->getNumAtom(); | 
| 285 | 
+ | 
         | 
| 286 | 
+ | 
      }         | 
| 287 | 
+ | 
    } | 
| 288 | 
+ | 
     | 
| 289 | 
+ | 
    return nAtoms_ - nLocalCutoffAtoms + nCutoffGroups_; | 
| 290 | 
  | 
  } | 
| 291 | 
  | 
     | 
| 292 | 
  | 
  void SimInfo::calcNdfRaw() { | 
| 759 | 
  | 
      if ( simParams_->getAccumulateBoxDipole() ) { | 
| 760 | 
  | 
        calcBoxDipole_ = true;        | 
| 761 | 
  | 
      } | 
| 762 | 
< | 
 | 
| 762 | 
> | 
     | 
| 763 | 
  | 
    set<AtomType*>::iterator i; | 
| 764 | 
  | 
    set<AtomType*> atomTypes; | 
| 765 | 
  | 
    atomTypes = getSimulatedAtomTypes();     | 
| 772 | 
  | 
      usesMetallic |= (*i)->isMetal(); | 
| 773 | 
  | 
      usesDirectional |= (*i)->isDirectional(); | 
| 774 | 
  | 
    } | 
| 775 | 
< | 
 | 
| 775 | 
> | 
     | 
| 776 | 
  | 
#ifdef IS_MPI     | 
| 777 | 
  | 
    int temp; | 
| 778 | 
  | 
    temp = usesDirectional; | 
| 779 | 
  | 
    MPI_Allreduce(&temp, &usesDirectionalAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);     | 
| 780 | 
< | 
 | 
| 780 | 
> | 
     | 
| 781 | 
  | 
    temp = usesMetallic; | 
| 782 | 
  | 
    MPI_Allreduce(&temp, &usesMetallicAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);     | 
| 783 | 
< | 
 | 
| 783 | 
> | 
     | 
| 784 | 
  | 
    temp = usesElectrostatic; | 
| 785 | 
  | 
    MPI_Allreduce(&temp, &usesElectrostaticAtoms_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  | 
| 786 | 
+ | 
#else | 
| 787 | 
+ | 
 | 
| 788 | 
+ | 
    usesDirectionalAtoms_ = usesDirectional; | 
| 789 | 
+ | 
    usesMetallicAtoms_ = usesMetallic; | 
| 790 | 
+ | 
    usesElectrostaticAtoms_ = usesElectrostatic; | 
| 791 | 
+ | 
 | 
| 792 | 
  | 
#endif | 
| 793 | 
< | 
    fInfo_.SIM_uses_PBC = usesPeriodicBoundaries_;     | 
| 794 | 
< | 
    fInfo_.SIM_uses_DirectionalAtoms = usesDirectionalAtoms_; | 
| 795 | 
< | 
    fInfo_.SIM_uses_MetallicAtoms = usesMetallicAtoms_; | 
| 796 | 
< | 
    fInfo_.SIM_requires_SkipCorrection = usesElectrostaticAtoms_; | 
| 784 | 
< | 
    fInfo_.SIM_requires_SelfCorrection = usesElectrostaticAtoms_; | 
| 785 | 
< | 
    fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; | 
| 793 | 
> | 
     | 
| 794 | 
> | 
    requiresPrepair_ = usesMetallicAtoms_ ? true : false;  | 
| 795 | 
> | 
    requiresSkipCorrection_ = usesElectrostaticAtoms_ ? true : false; | 
| 796 | 
> | 
    requiresSelfCorrection_ = usesElectrostaticAtoms_ ? true : false;     | 
| 797 | 
  | 
  } | 
| 798 | 
  | 
 | 
| 799 | 
  | 
 | 
| 836 | 
  | 
  } | 
| 837 | 
  | 
 | 
| 838 | 
  | 
 | 
| 839 | 
< | 
  void SimInfo::setupFortran() { | 
| 829 | 
< | 
    int isError; | 
| 839 | 
> | 
  void SimInfo::prepareTopology() { | 
| 840 | 
  | 
    int nExclude, nOneTwo, nOneThree, nOneFour; | 
| 831 | 
– | 
    vector<int> fortranGlobalGroupMembership; | 
| 832 | 
– | 
     | 
| 833 | 
– | 
    isError = 0; | 
| 841 | 
  | 
 | 
| 835 | 
– | 
    //globalGroupMembership_ is filled by SimCreator     | 
| 836 | 
– | 
    for (int i = 0; i < nGlobalAtoms_; i++) { | 
| 837 | 
– | 
      fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); | 
| 838 | 
– | 
    } | 
| 839 | 
– | 
 | 
| 842 | 
  | 
    //calculate mass ratio of cutoff group | 
| 841 | 
– | 
    vector<RealType> mfact; | 
| 843 | 
  | 
    SimInfo::MoleculeIterator mi; | 
| 844 | 
  | 
    Molecule* mol; | 
| 845 | 
  | 
    Molecule::CutoffGroupIterator ci; | 
| 848 | 
  | 
    Atom* atom; | 
| 849 | 
  | 
    RealType totalMass; | 
| 850 | 
  | 
 | 
| 851 | 
< | 
    //to avoid memory reallocation, reserve enough space for mfact | 
| 852 | 
< | 
    mfact.reserve(getNCutoffGroups()); | 
| 851 | 
> | 
    /** | 
| 852 | 
> | 
     * The mass factor is the relative mass of an atom to the total | 
| 853 | 
> | 
     * mass of the cutoff group it belongs to.  By default, all atoms | 
| 854 | 
> | 
     * are their own cutoff groups, and therefore have mass factors of | 
| 855 | 
> | 
     * 1.  We need some special handling for massless atoms, which | 
| 856 | 
> | 
     * will be treated as carrying the entire mass of the cutoff | 
| 857 | 
> | 
     * group. | 
| 858 | 
> | 
     */ | 
| 859 | 
> | 
    massFactors_.clear(); | 
| 860 | 
> | 
    massFactors_.resize(getNAtoms(), 1.0); | 
| 861 | 
  | 
     | 
| 862 | 
+ | 
    cerr << "mfs in si = " << massFactors_.size() << "\n"; | 
| 863 | 
  | 
    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {         | 
| 864 | 
< | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 864 | 
> | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL;  | 
| 865 | 
> | 
           cg = mol->nextCutoffGroup(ci)) { | 
| 866 | 
  | 
 | 
| 867 | 
  | 
        totalMass = cg->getMass(); | 
| 868 | 
  | 
        for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { | 
| 869 | 
  | 
          // Check for massless groups - set mfact to 1 if true | 
| 870 | 
< | 
          if (totalMass != 0) | 
| 871 | 
< | 
            mfact.push_back(atom->getMass()/totalMass); | 
| 870 | 
> | 
          if (totalMass != 0)  | 
| 871 | 
> | 
            massFactors_[atom->getLocalIndex()] = atom->getMass()/totalMass; | 
| 872 | 
  | 
          else | 
| 873 | 
< | 
            mfact.push_back( 1.0 ); | 
| 873 | 
> | 
            massFactors_[atom->getLocalIndex()] = 1.0; | 
| 874 | 
  | 
        } | 
| 875 | 
  | 
      }        | 
| 876 | 
  | 
    } | 
| 884 | 
  | 
        identArray_.push_back(atom->getIdent()); | 
| 885 | 
  | 
      } | 
| 886 | 
  | 
    }     | 
| 876 | 
– | 
 | 
| 877 | 
– | 
    //fill molMembershipArray | 
| 878 | 
– | 
    //molMembershipArray is filled by SimCreator     | 
| 879 | 
– | 
    vector<int> molMembershipArray(nGlobalAtoms_); | 
| 880 | 
– | 
    for (int i = 0; i < nGlobalAtoms_; i++) { | 
| 881 | 
– | 
      molMembershipArray[i] = globalMolMembership_[i] + 1; | 
| 882 | 
– | 
    } | 
| 887 | 
  | 
     | 
| 888 | 
< | 
    //setup fortran simulation | 
| 888 | 
> | 
    //scan topology  | 
| 889 | 
  | 
 | 
| 890 | 
  | 
    nExclude = excludedInteractions_.getSize(); | 
| 891 | 
  | 
    nOneTwo = oneTwoInteractions_.getSize(); | 
| 897 | 
  | 
    int* oneThreeList = oneThreeInteractions_.getPairList(); | 
| 898 | 
  | 
    int* oneFourList = oneFourInteractions_.getPairList(); | 
| 899 | 
  | 
 | 
| 900 | 
< | 
    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray_[0],  | 
| 897 | 
< | 
                   &nExclude, excludeList,  | 
| 898 | 
< | 
                   &nOneTwo, oneTwoList, | 
| 899 | 
< | 
                   &nOneThree, oneThreeList, | 
| 900 | 
< | 
                   &nOneFour, oneFourList, | 
| 901 | 
< | 
                   &molMembershipArray[0], &mfact[0], &nCutoffGroups_,  | 
| 902 | 
< | 
                   &fortranGlobalGroupMembership[0], &isError);  | 
| 903 | 
< | 
     | 
| 904 | 
< | 
    if( isError ){ | 
| 905 | 
< | 
       | 
| 906 | 
< | 
      sprintf( painCave.errMsg, | 
| 907 | 
< | 
               "There was an error setting the simulation information in fortran.\n" ); | 
| 908 | 
< | 
      painCave.isFatal = 1; | 
| 909 | 
< | 
      painCave.severity = OPENMD_ERROR; | 
| 910 | 
< | 
      simError(); | 
| 911 | 
< | 
    } | 
| 912 | 
< | 
     | 
| 913 | 
< | 
     | 
| 914 | 
< | 
    sprintf( checkPointMsg, | 
| 915 | 
< | 
             "succesfully sent the simulation information to fortran.\n"); | 
| 916 | 
< | 
     | 
| 917 | 
< | 
    errorCheckPoint(); | 
| 918 | 
< | 
     | 
| 919 | 
< | 
    // Setup number of neighbors in neighbor list if present | 
| 920 | 
< | 
    if (simParams_->haveNeighborListNeighbors()) { | 
| 921 | 
< | 
      int nlistNeighbors = simParams_->getNeighborListNeighbors(); | 
| 922 | 
< | 
      setNeighbors(&nlistNeighbors); | 
| 923 | 
< | 
    } | 
| 924 | 
< | 
    | 
| 925 | 
< | 
#ifdef IS_MPI     | 
| 926 | 
< | 
    mpiSimData parallelData; | 
| 927 | 
< | 
 | 
| 928 | 
< | 
    //fill up mpiSimData struct | 
| 929 | 
< | 
    parallelData.nMolGlobal = getNGlobalMolecules(); | 
| 930 | 
< | 
    parallelData.nMolLocal = getNMolecules(); | 
| 931 | 
< | 
    parallelData.nAtomsGlobal = getNGlobalAtoms(); | 
| 932 | 
< | 
    parallelData.nAtomsLocal = getNAtoms(); | 
| 933 | 
< | 
    parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); | 
| 934 | 
< | 
    parallelData.nGroupsLocal = getNCutoffGroups(); | 
| 935 | 
< | 
    parallelData.myNode = worldRank; | 
| 936 | 
< | 
    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); | 
| 937 | 
< | 
 | 
| 938 | 
< | 
    //pass mpiSimData struct and index arrays to fortran | 
| 939 | 
< | 
    //setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), | 
| 940 | 
< | 
    //                &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal), | 
| 941 | 
< | 
    //                &localToGlobalCutoffGroupIndex[0], &isError); | 
| 942 | 
< | 
 | 
| 943 | 
< | 
    if (isError) { | 
| 944 | 
< | 
      sprintf(painCave.errMsg, | 
| 945 | 
< | 
              "mpiRefresh errror: fortran didn't like something we gave it.\n"); | 
| 946 | 
< | 
      painCave.isFatal = 1; | 
| 947 | 
< | 
      simError(); | 
| 948 | 
< | 
    } | 
| 949 | 
< | 
 | 
| 950 | 
< | 
    sprintf(checkPointMsg, " mpiRefresh successful.\n"); | 
| 951 | 
< | 
    errorCheckPoint(); | 
| 952 | 
< | 
#endif | 
| 953 | 
< | 
 | 
| 954 | 
< | 
    initFortranFF(&isError); | 
| 955 | 
< | 
    if (isError) { | 
| 956 | 
< | 
      sprintf(painCave.errMsg, | 
| 957 | 
< | 
              "initFortranFF errror: fortran didn't like something we gave it.\n"); | 
| 958 | 
< | 
      painCave.isFatal = 1; | 
| 959 | 
< | 
      simError(); | 
| 960 | 
< | 
    } | 
| 961 | 
< | 
    fortranInitialized_ = true; | 
| 900 | 
> | 
    topologyDone_ = true; | 
| 901 | 
  | 
  } | 
| 902 | 
  | 
 | 
| 903 | 
  | 
  void SimInfo::addProperty(GenericData* genData) { |