| 54 | 
  | 
#include "math/Vector3.hpp" | 
| 55 | 
  | 
#include "primitives/Molecule.hpp" | 
| 56 | 
  | 
#include "primitives/StuntDouble.hpp" | 
| 57 | 
– | 
#include "UseTheForce/fCutoffPolicy.h" | 
| 58 | 
– | 
#include "UseTheForce/doForces_interface.h" | 
| 59 | 
– | 
#include "UseTheForce/DarkSide/neighborLists_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 | 
  | 
 | 
| 67 | 
– | 
 | 
| 68 | 
– | 
#ifdef IS_MPI | 
| 69 | 
– | 
#include "UseTheForce/mpiComponentPlan.h" | 
| 70 | 
– | 
#include "UseTheForce/DarkSide/simParallel_interface.h" | 
| 71 | 
– | 
#endif  | 
| 72 | 
– | 
 | 
| 64 | 
  | 
using namespace std; | 
| 65 | 
  | 
namespace OpenMD { | 
| 66 | 
  | 
   | 
| 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 | 
| 128 | 
+ | 
    std::cerr << "nGA = " << nGlobalAtoms_ << "\n"; | 
| 129 | 
+ | 
    std::cerr << "nCA = " << nCutoffAtoms << "\n"; | 
| 130 | 
+ | 
    std::cerr << "nG = " << nGroups << "\n"; | 
| 131 | 
+ | 
 | 
| 132 | 
  | 
    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; | 
| 133 | 
+ | 
 | 
| 134 | 
+ | 
    std::cerr << "nGCG = " << nGlobalCutoffGroups_ << "\n"; | 
| 135 | 
  | 
     | 
| 136 | 
  | 
    //every free atom (atom does not belong to rigid bodies) is an  | 
| 137 | 
  | 
    //integrable object therefore the total number of integrable objects  | 
| 655 | 
  | 
  /** | 
| 656 | 
  | 
   * update | 
| 657 | 
  | 
   * | 
| 658 | 
< | 
   *  Performs the global checks and variable settings after the objects have been | 
| 659 | 
< | 
   *  created.  | 
| 658 | 
> | 
   *  Performs the global checks and variable settings after the | 
| 659 | 
> | 
   *  objects have been created. | 
| 660 | 
  | 
   *  | 
| 661 | 
  | 
   */ | 
| 662 | 
< | 
  void SimInfo::update() { | 
| 666 | 
< | 
     | 
| 662 | 
> | 
  void SimInfo::update() {    | 
| 663 | 
  | 
    setupSimVariables(); | 
| 668 | 
– | 
    setupCutoffs(); | 
| 669 | 
– | 
    setupSwitching(); | 
| 670 | 
– | 
    setupElectrostatics(); | 
| 671 | 
– | 
    setupNeighborlists(); | 
| 672 | 
– | 
 | 
| 673 | 
– | 
#ifdef IS_MPI | 
| 674 | 
– | 
    setupFortranParallel(); | 
| 675 | 
– | 
#endif | 
| 676 | 
– | 
    setupFortranSim(); | 
| 677 | 
– | 
    fortranInitialized_ = true; | 
| 678 | 
– | 
 | 
| 664 | 
  | 
    calcNdf(); | 
| 665 | 
  | 
    calcNdfRaw(); | 
| 666 | 
  | 
    calcNdfTrans(); | 
| 667 | 
  | 
  } | 
| 668 | 
  | 
   | 
| 669 | 
+ | 
  /** | 
| 670 | 
+ | 
   * getSimulatedAtomTypes | 
| 671 | 
+ | 
   * | 
| 672 | 
+ | 
   * Returns an STL set of AtomType* that are actually present in this | 
| 673 | 
+ | 
   * simulation.  Must query all processors to assemble this information. | 
| 674 | 
+ | 
   *  | 
| 675 | 
+ | 
   */ | 
| 676 | 
  | 
  set<AtomType*> SimInfo::getSimulatedAtomTypes() { | 
| 677 | 
  | 
    SimInfo::MoleculeIterator mi; | 
| 678 | 
  | 
    Molecule* mol; | 
| 685 | 
  | 
        atomTypes.insert(atom->getAtomType()); | 
| 686 | 
  | 
      }       | 
| 687 | 
  | 
    }     | 
| 696 | 
– | 
    return atomTypes;         | 
| 697 | 
– | 
  } | 
| 688 | 
  | 
 | 
| 689 | 
< | 
  /** | 
| 700 | 
< | 
   * setupCutoffs | 
| 701 | 
< | 
   * | 
| 702 | 
< | 
   * Sets the values of cutoffRadius and cutoffMethod | 
| 703 | 
< | 
   * | 
| 704 | 
< | 
   * cutoffRadius : realType | 
| 705 | 
< | 
   *  If the cutoffRadius was explicitly set, use that value. | 
| 706 | 
< | 
   *  If the cutoffRadius was not explicitly set: | 
| 707 | 
< | 
   *      Are there electrostatic atoms?  Use 12.0 Angstroms. | 
| 708 | 
< | 
   *      No electrostatic atoms?  Poll the atom types present in the | 
| 709 | 
< | 
   *      simulation for suggested cutoff values (e.g. 2.5 * sigma). | 
| 710 | 
< | 
   *      Use the maximum suggested value that was found. | 
| 711 | 
< | 
   * | 
| 712 | 
< | 
   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL) | 
| 713 | 
< | 
   *      If cutoffMethod was explicitly set, use that choice. | 
| 714 | 
< | 
   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE | 
| 715 | 
< | 
   */ | 
| 716 | 
< | 
  void SimInfo::setupCutoffs() { | 
| 717 | 
< | 
     | 
| 718 | 
< | 
    if (simParams_->haveCutoffRadius()) { | 
| 719 | 
< | 
      cutoffRadius_ = simParams_->getCutoffRadius(); | 
| 720 | 
< | 
    } else {       | 
| 721 | 
< | 
      if (usesElectrostaticAtoms_) { | 
| 722 | 
< | 
        sprintf(painCave.errMsg, | 
| 723 | 
< | 
                "SimInfo: No value was set for the cutoffRadius.\n" | 
| 724 | 
< | 
                "\tOpenMD will use a default value of 12.0 angstroms" | 
| 725 | 
< | 
                "\tfor the cutoffRadius.\n"); | 
| 726 | 
< | 
        painCave.isFatal = 0; | 
| 727 | 
< | 
        painCave.severity = OPENMD_INFO; | 
| 728 | 
< | 
        simError(); | 
| 729 | 
< | 
        cutoffRadius_ = 12.0; | 
| 730 | 
< | 
      } else { | 
| 731 | 
< | 
        RealType thisCut; | 
| 732 | 
< | 
        set<AtomType*>::iterator i; | 
| 733 | 
< | 
        set<AtomType*> atomTypes; | 
| 734 | 
< | 
        atomTypes = getSimulatedAtomTypes();         | 
| 735 | 
< | 
        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { | 
| 736 | 
< | 
          thisCut = InteractionManager::Instance()->getSuggestedCutoffRadius((*i)); | 
| 737 | 
< | 
          cutoffRadius_ = max(thisCut, cutoffRadius_); | 
| 738 | 
< | 
        } | 
| 739 | 
< | 
        sprintf(painCave.errMsg, | 
| 740 | 
< | 
                "SimInfo: No value was set for the cutoffRadius.\n" | 
| 741 | 
< | 
                "\tOpenMD will use %lf angstroms.\n", | 
| 742 | 
< | 
                cutoffRadius_); | 
| 743 | 
< | 
        painCave.isFatal = 0; | 
| 744 | 
< | 
        painCave.severity = OPENMD_INFO; | 
| 745 | 
< | 
        simError(); | 
| 746 | 
< | 
      }              | 
| 747 | 
< | 
    } | 
| 689 | 
> | 
#ifdef IS_MPI | 
| 690 | 
  | 
 | 
| 691 | 
< | 
    InteractionManager::Instance()->setCutoffRadius(cutoffRadius_); | 
| 692 | 
< | 
 | 
| 751 | 
< | 
    map<string, CutoffMethod> stringToCutoffMethod; | 
| 752 | 
< | 
    stringToCutoffMethod["HARD"] = HARD; | 
| 753 | 
< | 
    stringToCutoffMethod["SWITCHING_FUNCTION"] = SWITCHING_FUNCTION; | 
| 754 | 
< | 
    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;     | 
| 755 | 
< | 
    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE; | 
| 756 | 
< | 
   | 
| 757 | 
< | 
    if (simParams_->haveCutoffMethod()) { | 
| 758 | 
< | 
      string cutMeth = toUpperCopy(simParams_->getCutoffMethod()); | 
| 759 | 
< | 
      map<string, CutoffMethod>::iterator i; | 
| 760 | 
< | 
      i = stringToCutoffMethod.find(cutMeth); | 
| 761 | 
< | 
      if (i == stringToCutoffMethod.end()) { | 
| 762 | 
< | 
        sprintf(painCave.errMsg, | 
| 763 | 
< | 
                "SimInfo: Could not find chosen cutoffMethod %s\n" | 
| 764 | 
< | 
                "\tShould be one of: " | 
| 765 | 
< | 
                "HARD, SWITCHING_FUNCTION, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n", | 
| 766 | 
< | 
                cutMeth.c_str()); | 
| 767 | 
< | 
        painCave.isFatal = 1; | 
| 768 | 
< | 
        painCave.severity = OPENMD_ERROR; | 
| 769 | 
< | 
        simError(); | 
| 770 | 
< | 
      } else { | 
| 771 | 
< | 
        cutoffMethod_ = i->second; | 
| 772 | 
< | 
      } | 
| 773 | 
< | 
    } else { | 
| 774 | 
< | 
      sprintf(painCave.errMsg, | 
| 775 | 
< | 
              "SimInfo: No value was set for the cutoffMethod.\n" | 
| 776 | 
< | 
              "\tOpenMD will use SHIFTED_FORCE.\n"); | 
| 777 | 
< | 
        painCave.isFatal = 0; | 
| 778 | 
< | 
        painCave.severity = OPENMD_INFO; | 
| 779 | 
< | 
        simError(); | 
| 780 | 
< | 
        cutoffMethod_ = SHIFTED_FORCE;         | 
| 781 | 
< | 
    } | 
| 691 | 
> | 
    // loop over the found atom types on this processor, and add their | 
| 692 | 
> | 
    // numerical idents to a vector: | 
| 693 | 
  | 
 | 
| 694 | 
< | 
    InteractionManager::Instance()->setCutoffMethod(cutoffMethod_); | 
| 695 | 
< | 
  } | 
| 696 | 
< | 
   | 
| 697 | 
< | 
  /** | 
| 787 | 
< | 
   * setupSwitching | 
| 788 | 
< | 
   * | 
| 789 | 
< | 
   * Sets the values of switchingRadius and  | 
| 790 | 
< | 
   *  If the switchingRadius was explicitly set, use that value (but check it) | 
| 791 | 
< | 
   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_ | 
| 792 | 
< | 
   */ | 
| 793 | 
< | 
  void SimInfo::setupSwitching() { | 
| 794 | 
< | 
     | 
| 795 | 
< | 
    if (simParams_->haveSwitchingRadius()) { | 
| 796 | 
< | 
      switchingRadius_ = simParams_->getSwitchingRadius(); | 
| 797 | 
< | 
      if (switchingRadius_ > cutoffRadius_) {         | 
| 798 | 
< | 
        sprintf(painCave.errMsg, | 
| 799 | 
< | 
                "SimInfo: switchingRadius (%f) is larger than cutoffRadius(%f)\n", | 
| 800 | 
< | 
                switchingRadius_, cutoffRadius_); | 
| 801 | 
< | 
        painCave.isFatal = 1; | 
| 802 | 
< | 
        painCave.severity = OPENMD_ERROR; | 
| 803 | 
< | 
        simError(); | 
| 804 | 
< | 
      } | 
| 805 | 
< | 
    } else {       | 
| 806 | 
< | 
      switchingRadius_ = 0.85 * cutoffRadius_; | 
| 807 | 
< | 
      sprintf(painCave.errMsg, | 
| 808 | 
< | 
              "SimInfo: No value was set for the switchingRadius.\n" | 
| 809 | 
< | 
              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n" | 
| 810 | 
< | 
              "\tswitchingRadius = %f. for this simulation\n", switchingRadius_); | 
| 811 | 
< | 
      painCave.isFatal = 0; | 
| 812 | 
< | 
      painCave.severity = OPENMD_WARNING; | 
| 813 | 
< | 
      simError(); | 
| 814 | 
< | 
    }            | 
| 815 | 
< | 
   | 
| 816 | 
< | 
    InteractionManager::Instance()->setSwitchingRadius(switchingRadius_); | 
| 694 | 
> | 
    vector<int> foundTypes; | 
| 695 | 
> | 
    set<AtomType*>::iterator i; | 
| 696 | 
> | 
    for (i = atomTypes.begin(); i != atomTypes.end(); ++i)  | 
| 697 | 
> | 
      foundTypes.push_back( (*i)->getIdent() ); | 
| 698 | 
  | 
 | 
| 699 | 
< | 
    SwitchingFunctionType ft; | 
| 699 | 
> | 
    // count_local holds the number of found types on this processor | 
| 700 | 
> | 
    int count_local = foundTypes.size(); | 
| 701 | 
> | 
 | 
| 702 | 
> | 
    // count holds the total number of found types on all processors | 
| 703 | 
> | 
    // (some will be redundant with the ones found locally): | 
| 704 | 
> | 
    int count; | 
| 705 | 
> | 
    MPI::COMM_WORLD.Allreduce(&count_local, &count, 1, MPI::INT, MPI::SUM); | 
| 706 | 
> | 
 | 
| 707 | 
> | 
    // create a vector to hold the globally found types, and resize it: | 
| 708 | 
> | 
    vector<int> ftGlobal; | 
| 709 | 
> | 
    ftGlobal.resize(count); | 
| 710 | 
> | 
    vector<int> counts; | 
| 711 | 
> | 
 | 
| 712 | 
> | 
    int nproc = MPI::COMM_WORLD.Get_size(); | 
| 713 | 
> | 
    counts.resize(nproc); | 
| 714 | 
> | 
    vector<int> disps; | 
| 715 | 
> | 
    disps.resize(nproc); | 
| 716 | 
> | 
 | 
| 717 | 
> | 
    // now spray out the foundTypes to all the other processors: | 
| 718 | 
  | 
     | 
| 719 | 
< | 
    if (simParams_->haveSwitchingFunctionType()) { | 
| 720 | 
< | 
      string funcType = simParams_->getSwitchingFunctionType(); | 
| 822 | 
< | 
      toUpper(funcType); | 
| 823 | 
< | 
      if (funcType == "CUBIC") { | 
| 824 | 
< | 
        ft = cubic; | 
| 825 | 
< | 
      } else { | 
| 826 | 
< | 
        if (funcType == "FIFTH_ORDER_POLYNOMIAL") { | 
| 827 | 
< | 
          ft = fifth_order_poly; | 
| 828 | 
< | 
        } else { | 
| 829 | 
< | 
          // throw error         | 
| 830 | 
< | 
          sprintf( painCave.errMsg, | 
| 831 | 
< | 
                   "SimInfo : Unknown switchingFunctionType. (Input file specified %s .)\n" | 
| 832 | 
< | 
                   "\tswitchingFunctionType must be one of: " | 
| 833 | 
< | 
                   "\"cubic\" or \"fifth_order_polynomial\".",  | 
| 834 | 
< | 
                   funcType.c_str() ); | 
| 835 | 
< | 
          painCave.isFatal = 1; | 
| 836 | 
< | 
          painCave.severity = OPENMD_ERROR; | 
| 837 | 
< | 
          simError(); | 
| 838 | 
< | 
        }            | 
| 839 | 
< | 
      } | 
| 840 | 
< | 
    } | 
| 719 | 
> | 
    MPI::COMM_WORLD.Allgatherv(&foundTypes[0], count_local, MPI::INT,  | 
| 720 | 
> | 
                               &ftGlobal[0], &counts[0], &disps[0], MPI::INT); | 
| 721 | 
  | 
 | 
| 722 | 
< | 
    InteractionManager::Instance()->setSwitchingFunctionType(ft); | 
| 722 | 
> | 
    // foundIdents is a stl set, so inserting an already found ident | 
| 723 | 
> | 
    // will have no effect. | 
| 724 | 
> | 
    set<int> foundIdents; | 
| 725 | 
> | 
    vector<int>::iterator j; | 
| 726 | 
> | 
    for (j = ftGlobal.begin(); j != ftGlobal.end(); ++j) | 
| 727 | 
> | 
      foundIdents.insert((*j)); | 
| 728 | 
> | 
     | 
| 729 | 
> | 
    // now iterate over the foundIdents and get the actual atom types  | 
| 730 | 
> | 
    // that correspond to these: | 
| 731 | 
> | 
    set<int>::iterator it; | 
| 732 | 
> | 
    for (it = foundIdents.begin(); it != foundIdents.end(); ++it) | 
| 733 | 
> | 
      atomTypes.insert( forceField_->getAtomType((*it)) ); | 
| 734 | 
> | 
  | 
| 735 | 
> | 
#endif | 
| 736 | 
> | 
     | 
| 737 | 
> | 
    return atomTypes;         | 
| 738 | 
  | 
  } | 
| 739 | 
  | 
 | 
| 740 | 
< | 
  /** | 
| 741 | 
< | 
   * setupSkinThickness | 
| 742 | 
< | 
   * | 
| 743 | 
< | 
   *  If the skinThickness was explicitly set, use that value (but check it) | 
| 744 | 
< | 
   *  If the skinThickness was not explicitly set: use 1.0 angstroms | 
| 745 | 
< | 
   */ | 
| 746 | 
< | 
  void SimInfo::setupSkinThickness() {     | 
| 747 | 
< | 
    if (simParams_->haveSkinThickness()) { | 
| 853 | 
< | 
      skinThickness_ = simParams_->getSkinThickness(); | 
| 854 | 
< | 
    } else {       | 
| 855 | 
< | 
      skinThickness_ = 1.0; | 
| 856 | 
< | 
      sprintf(painCave.errMsg, | 
| 857 | 
< | 
              "SimInfo Warning: No value was set for the skinThickness.\n" | 
| 858 | 
< | 
              "\tOpenMD will use a default value of %f Angstroms\n" | 
| 859 | 
< | 
              "\tfor this simulation\n", skinThickness_); | 
| 860 | 
< | 
      painCave.isFatal = 0; | 
| 861 | 
< | 
      simError(); | 
| 862 | 
< | 
    }              | 
| 863 | 
< | 
  } | 
| 740 | 
> | 
  void SimInfo::setupSimVariables() { | 
| 741 | 
> | 
    useAtomicVirial_ = simParams_->getUseAtomicVirial(); | 
| 742 | 
> | 
    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 743 | 
> | 
    calcBoxDipole_ = false; | 
| 744 | 
> | 
    if ( simParams_->haveAccumulateBoxDipole() )  | 
| 745 | 
> | 
      if ( simParams_->getAccumulateBoxDipole() ) { | 
| 746 | 
> | 
        calcBoxDipole_ = true;        | 
| 747 | 
> | 
      } | 
| 748 | 
  | 
 | 
| 865 | 
– | 
  void SimInfo::setupSimType() { | 
| 749 | 
  | 
    set<AtomType*>::iterator i; | 
| 750 | 
  | 
    set<AtomType*> atomTypes; | 
| 751 | 
< | 
    atomTypes = getSimulatedAtomTypes(); | 
| 869 | 
< | 
 | 
| 870 | 
< | 
    useAtomicVirial_ = simParams_->getUseAtomicVirial(); | 
| 871 | 
< | 
 | 
| 751 | 
> | 
    atomTypes = getSimulatedAtomTypes();     | 
| 752 | 
  | 
    int usesElectrostatic = 0; | 
| 753 | 
  | 
    int usesMetallic = 0; | 
| 754 | 
  | 
    int usesDirectional = 0; | 
| 778 | 
  | 
    fInfo_.SIM_uses_AtomicVirial = usesAtomicVirial_; | 
| 779 | 
  | 
  } | 
| 780 | 
  | 
 | 
| 781 | 
< | 
  void SimInfo::setupFortranSim() { | 
| 781 | 
> | 
 | 
| 782 | 
> | 
  vector<int> SimInfo::getGlobalAtomIndices() { | 
| 783 | 
> | 
    SimInfo::MoleculeIterator mi; | 
| 784 | 
> | 
    Molecule* mol; | 
| 785 | 
> | 
    Molecule::AtomIterator ai; | 
| 786 | 
> | 
    Atom* atom; | 
| 787 | 
> | 
 | 
| 788 | 
> | 
    vector<int> GlobalAtomIndices(getNAtoms(), 0); | 
| 789 | 
> | 
     | 
| 790 | 
> | 
    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) { | 
| 791 | 
> | 
       | 
| 792 | 
> | 
      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 793 | 
> | 
        GlobalAtomIndices[atom->getLocalIndex()] = atom->getGlobalIndex(); | 
| 794 | 
> | 
      } | 
| 795 | 
> | 
    } | 
| 796 | 
> | 
    return GlobalAtomIndices; | 
| 797 | 
> | 
  } | 
| 798 | 
> | 
 | 
| 799 | 
> | 
 | 
| 800 | 
> | 
  vector<int> SimInfo::getGlobalGroupIndices() { | 
| 801 | 
> | 
    SimInfo::MoleculeIterator mi; | 
| 802 | 
> | 
    Molecule* mol; | 
| 803 | 
> | 
    Molecule::CutoffGroupIterator ci; | 
| 804 | 
> | 
    CutoffGroup* cg; | 
| 805 | 
> | 
 | 
| 806 | 
> | 
    vector<int> GlobalGroupIndices; | 
| 807 | 
> | 
     | 
| 808 | 
> | 
    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) { | 
| 809 | 
> | 
       | 
| 810 | 
> | 
      //local index of cutoff group is trivial, it only depends on the | 
| 811 | 
> | 
      //order of travesing | 
| 812 | 
> | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL;  | 
| 813 | 
> | 
           cg = mol->nextCutoffGroup(ci)) { | 
| 814 | 
> | 
        GlobalGroupIndices.push_back(cg->getGlobalIndex()); | 
| 815 | 
> | 
      }         | 
| 816 | 
> | 
    } | 
| 817 | 
> | 
    return GlobalGroupIndices; | 
| 818 | 
> | 
  } | 
| 819 | 
> | 
 | 
| 820 | 
> | 
 | 
| 821 | 
> | 
  void SimInfo::setupFortran() { | 
| 822 | 
  | 
    int isError; | 
| 823 | 
  | 
    int nExclude, nOneTwo, nOneThree, nOneFour; | 
| 824 | 
  | 
    vector<int> fortranGlobalGroupMembership; | 
| 825 | 
  | 
     | 
| 906 | 
– | 
    notifyFortranSkinThickness(&skinThickness_); | 
| 907 | 
– | 
 | 
| 908 | 
– | 
    int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0; | 
| 909 | 
– | 
    int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0; | 
| 910 | 
– | 
    notifyFortranCutoffs(&cutoffRadius_, &switchingRadius_, &ljsp, &ljsf); | 
| 911 | 
– | 
 | 
| 826 | 
  | 
    isError = 0; | 
| 827 | 
  | 
 | 
| 828 | 
  | 
    //globalGroupMembership_ is filled by SimCreator     | 
| 857 | 
  | 
      }        | 
| 858 | 
  | 
    } | 
| 859 | 
  | 
 | 
| 860 | 
< | 
    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) | 
| 947 | 
< | 
    vector<int> identArray; | 
| 860 | 
> | 
    // Build the identArray_ | 
| 861 | 
  | 
 | 
| 862 | 
< | 
    //to avoid memory reallocation, reserve enough space identArray | 
| 863 | 
< | 
    identArray.reserve(getNAtoms()); | 
| 951 | 
< | 
     | 
| 862 | 
> | 
    identArray_.clear(); | 
| 863 | 
> | 
    identArray_.reserve(getNAtoms());     | 
| 864 | 
  | 
    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {         | 
| 865 | 
  | 
      for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 866 | 
< | 
        identArray.push_back(atom->getIdent()); | 
| 866 | 
> | 
        identArray_.push_back(atom->getIdent()); | 
| 867 | 
  | 
      } | 
| 868 | 
  | 
    }     | 
| 869 | 
  | 
 | 
| 886 | 
  | 
    int* oneThreeList = oneThreeInteractions_.getPairList(); | 
| 887 | 
  | 
    int* oneFourList = oneFourInteractions_.getPairList(); | 
| 888 | 
  | 
 | 
| 889 | 
< | 
    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0],  | 
| 890 | 
< | 
                   &nExclude, excludeList,  | 
| 891 | 
< | 
                   &nOneTwo, oneTwoList, | 
| 892 | 
< | 
                   &nOneThree, oneThreeList, | 
| 893 | 
< | 
                   &nOneFour, oneFourList, | 
| 894 | 
< | 
                   &molMembershipArray[0], &mfact[0], &nCutoffGroups_,  | 
| 895 | 
< | 
                   &fortranGlobalGroupMembership[0], &isError);  | 
| 889 | 
> | 
    //setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray_[0],  | 
| 890 | 
> | 
    //               &nExclude, excludeList,  | 
| 891 | 
> | 
    //               &nOneTwo, oneTwoList, | 
| 892 | 
> | 
    //               &nOneThree, oneThreeList, | 
| 893 | 
> | 
    //               &nOneFour, oneFourList, | 
| 894 | 
> | 
    //               &molMembershipArray[0], &mfact[0], &nCutoffGroups_,  | 
| 895 | 
> | 
    //               &fortranGlobalGroupMembership[0], &isError);  | 
| 896 | 
  | 
     | 
| 897 | 
< | 
    if( isError ){ | 
| 898 | 
< | 
       | 
| 899 | 
< | 
      sprintf( painCave.errMsg, | 
| 900 | 
< | 
               "There was an error setting the simulation information in fortran.\n" ); | 
| 901 | 
< | 
      painCave.isFatal = 1; | 
| 902 | 
< | 
      painCave.severity = OPENMD_ERROR; | 
| 903 | 
< | 
      simError(); | 
| 904 | 
< | 
    } | 
| 897 | 
> | 
    // if( isError ){ | 
| 898 | 
> | 
    //   | 
| 899 | 
> | 
    //  sprintf( painCave.errMsg, | 
| 900 | 
> | 
    //         "There was an error setting the simulation information in fortran.\n" ); | 
| 901 | 
> | 
    //  painCave.isFatal = 1; | 
| 902 | 
> | 
    //  painCave.severity = OPENMD_ERROR; | 
| 903 | 
> | 
    //  simError(); | 
| 904 | 
> | 
    //} | 
| 905 | 
  | 
     | 
| 906 | 
  | 
     | 
| 907 | 
< | 
    sprintf( checkPointMsg, | 
| 908 | 
< | 
             "succesfully sent the simulation information to fortran.\n"); | 
| 907 | 
> | 
    // sprintf( checkPointMsg, | 
| 908 | 
> | 
    //          "succesfully sent the simulation information to fortran.\n"); | 
| 909 | 
  | 
     | 
| 910 | 
< | 
    errorCheckPoint(); | 
| 910 | 
> | 
    // errorCheckPoint(); | 
| 911 | 
  | 
     | 
| 912 | 
  | 
    // Setup number of neighbors in neighbor list if present | 
| 913 | 
< | 
    if (simParams_->haveNeighborListNeighbors()) { | 
| 914 | 
< | 
      int nlistNeighbors = simParams_->getNeighborListNeighbors(); | 
| 915 | 
< | 
      setNeighbors(&nlistNeighbors); | 
| 916 | 
< | 
    } | 
| 917 | 
< | 
    | 
| 1006 | 
< | 
 | 
| 1007 | 
< | 
  } | 
| 1008 | 
< | 
 | 
| 1009 | 
< | 
 | 
| 1010 | 
< | 
  void SimInfo::setupFortranParallel() { | 
| 913 | 
> | 
    //if (simParams_->haveNeighborListNeighbors()) { | 
| 914 | 
> | 
    //  int nlistNeighbors = simParams_->getNeighborListNeighbors(); | 
| 915 | 
> | 
    //  setNeighbors(&nlistNeighbors); | 
| 916 | 
> | 
    //} | 
| 917 | 
> | 
    | 
| 918 | 
  | 
#ifdef IS_MPI     | 
| 919 | 
< | 
    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex | 
| 1013 | 
< | 
    vector<int> localToGlobalAtomIndex(getNAtoms(), 0); | 
| 1014 | 
< | 
    vector<int> localToGlobalCutoffGroupIndex; | 
| 1015 | 
< | 
    SimInfo::MoleculeIterator mi; | 
| 1016 | 
< | 
    Molecule::AtomIterator ai; | 
| 1017 | 
< | 
    Molecule::CutoffGroupIterator ci; | 
| 1018 | 
< | 
    Molecule* mol; | 
| 1019 | 
< | 
    Atom* atom; | 
| 1020 | 
< | 
    CutoffGroup* cg; | 
| 1021 | 
< | 
    mpiSimData parallelData; | 
| 1022 | 
< | 
    int isError; | 
| 919 | 
> | 
    // mpiSimData parallelData; | 
| 920 | 
  | 
 | 
| 1024 | 
– | 
    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) { | 
| 1025 | 
– | 
 | 
| 1026 | 
– | 
      //local index(index in DataStorge) of atom is important | 
| 1027 | 
– | 
      for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { | 
| 1028 | 
– | 
        localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; | 
| 1029 | 
– | 
      } | 
| 1030 | 
– | 
 | 
| 1031 | 
– | 
      //local index of cutoff group is trivial, it only depends on the order of travesing | 
| 1032 | 
– | 
      for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { | 
| 1033 | 
– | 
        localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); | 
| 1034 | 
– | 
      }         | 
| 1035 | 
– | 
         | 
| 1036 | 
– | 
    } | 
| 1037 | 
– | 
 | 
| 921 | 
  | 
    //fill up mpiSimData struct | 
| 922 | 
< | 
    parallelData.nMolGlobal = getNGlobalMolecules(); | 
| 923 | 
< | 
    parallelData.nMolLocal = getNMolecules(); | 
| 924 | 
< | 
    parallelData.nAtomsGlobal = getNGlobalAtoms(); | 
| 925 | 
< | 
    parallelData.nAtomsLocal = getNAtoms(); | 
| 926 | 
< | 
    parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); | 
| 927 | 
< | 
    parallelData.nGroupsLocal = getNCutoffGroups(); | 
| 928 | 
< | 
    parallelData.myNode = worldRank; | 
| 929 | 
< | 
    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); | 
| 922 | 
> | 
    // parallelData.nMolGlobal = getNGlobalMolecules(); | 
| 923 | 
> | 
    // parallelData.nMolLocal = getNMolecules(); | 
| 924 | 
> | 
    // parallelData.nAtomsGlobal = getNGlobalAtoms(); | 
| 925 | 
> | 
    // parallelData.nAtomsLocal = getNAtoms(); | 
| 926 | 
> | 
    // parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); | 
| 927 | 
> | 
    // parallelData.nGroupsLocal = getNCutoffGroups(); | 
| 928 | 
> | 
    // parallelData.myNode = worldRank; | 
| 929 | 
> | 
    // MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); | 
| 930 | 
  | 
 | 
| 931 | 
  | 
    //pass mpiSimData struct and index arrays to fortran | 
| 932 | 
< | 
    setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), | 
| 933 | 
< | 
                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal), | 
| 934 | 
< | 
                    &localToGlobalCutoffGroupIndex[0], &isError); | 
| 932 | 
> | 
    //setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), | 
| 933 | 
> | 
    //                &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal), | 
| 934 | 
> | 
    //                &localToGlobalCutoffGroupIndex[0], &isError); | 
| 935 | 
  | 
 | 
| 936 | 
< | 
    if (isError) { | 
| 937 | 
< | 
      sprintf(painCave.errMsg, | 
| 938 | 
< | 
              "mpiRefresh errror: fortran didn't like something we gave it.\n"); | 
| 939 | 
< | 
      painCave.isFatal = 1; | 
| 940 | 
< | 
      simError(); | 
| 941 | 
< | 
    } | 
| 936 | 
> | 
    // if (isError) { | 
| 937 | 
> | 
    //   sprintf(painCave.errMsg, | 
| 938 | 
> | 
    //           "mpiRefresh errror: fortran didn't like something we gave it.\n"); | 
| 939 | 
> | 
    //   painCave.isFatal = 1; | 
| 940 | 
> | 
    //   simError(); | 
| 941 | 
> | 
    // } | 
| 942 | 
  | 
 | 
| 943 | 
< | 
    sprintf(checkPointMsg, " mpiRefresh successful.\n"); | 
| 944 | 
< | 
    errorCheckPoint(); | 
| 1062 | 
< | 
 | 
| 943 | 
> | 
    // sprintf(checkPointMsg, " mpiRefresh successful.\n"); | 
| 944 | 
> | 
    // errorCheckPoint(); | 
| 945 | 
  | 
#endif | 
| 1064 | 
– | 
  } | 
| 946 | 
  | 
 | 
| 947 | 
< | 
 | 
| 948 | 
< | 
  void SimInfo::setupSwitchingFunction() {     | 
| 949 | 
< | 
 | 
| 950 | 
< | 
  } | 
| 951 | 
< | 
 | 
| 952 | 
< | 
  void SimInfo::setupAccumulateBoxDipole() {     | 
| 953 | 
< | 
 | 
| 954 | 
< | 
    // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true | 
| 1074 | 
< | 
    if ( simParams_->haveAccumulateBoxDipole() )  | 
| 1075 | 
< | 
      if ( simParams_->getAccumulateBoxDipole() ) { | 
| 1076 | 
< | 
        calcBoxDipole_ = true; | 
| 1077 | 
< | 
      } | 
| 1078 | 
< | 
 | 
| 947 | 
> | 
    // initFortranFF(&isError); | 
| 948 | 
> | 
    // if (isError) { | 
| 949 | 
> | 
    //   sprintf(painCave.errMsg, | 
| 950 | 
> | 
    //           "initFortranFF errror: fortran didn't like something we gave it.\n"); | 
| 951 | 
> | 
    //   painCave.isFatal = 1; | 
| 952 | 
> | 
    //   simError(); | 
| 953 | 
> | 
    // } | 
| 954 | 
> | 
    // fortranInitialized_ = true; | 
| 955 | 
  | 
  } | 
| 956 | 
  | 
 | 
| 957 | 
  | 
  void SimInfo::addProperty(GenericData* genData) { | 
| 988 | 
  | 
    Molecule* mol; | 
| 989 | 
  | 
    RigidBody* rb; | 
| 990 | 
  | 
    Atom* atom; | 
| 991 | 
+ | 
    CutoffGroup* cg; | 
| 992 | 
  | 
    SimInfo::MoleculeIterator mi; | 
| 993 | 
  | 
    Molecule::RigidBodyIterator rbIter; | 
| 994 | 
< | 
    Molecule::AtomIterator atomIter;; | 
| 994 | 
> | 
    Molecule::AtomIterator atomIter; | 
| 995 | 
> | 
    Molecule::CutoffGroupIterator cgIter; | 
| 996 | 
  | 
  | 
| 997 | 
  | 
    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { | 
| 998 | 
  | 
         | 
| 1003 | 
  | 
      for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { | 
| 1004 | 
  | 
        rb->setSnapshotManager(sman_); | 
| 1005 | 
  | 
      } | 
| 1006 | 
+ | 
 | 
| 1007 | 
+ | 
      for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) { | 
| 1008 | 
+ | 
        cg->setSnapshotManager(sman_); | 
| 1009 | 
+ | 
      } | 
| 1010 | 
  | 
    }     | 
| 1011 | 
  | 
     | 
| 1012 | 
  | 
  } |