| 49 | 
  | 
 */ | 
| 50 | 
  | 
   | 
| 51 | 
  | 
#include <algorithm> | 
| 52 | 
< | 
#include "UseTheForce/ForceField.hpp" | 
| 52 | 
> | 
#include "brains/ForceField.hpp" | 
| 53 | 
  | 
#include "utils/simError.h" | 
| 54 | 
< | 
#include "utils/Tuple.hpp" | 
| 54 | 
> | 
 | 
| 55 | 
> | 
#include "io/OptionSectionParser.hpp" | 
| 56 | 
> | 
#include "io/BaseAtomTypesSectionParser.hpp" | 
| 57 | 
> | 
#include "io/DirectionalAtomTypesSectionParser.hpp" | 
| 58 | 
> | 
#include "io/AtomTypesSectionParser.hpp" | 
| 59 | 
> | 
#include "io/BendTypesSectionParser.hpp" | 
| 60 | 
> | 
#include "io/BondTypesSectionParser.hpp" | 
| 61 | 
> | 
#include "io/ChargeAtomTypesSectionParser.hpp" | 
| 62 | 
> | 
#include "io/EAMAtomTypesSectionParser.hpp" | 
| 63 | 
> | 
#include "io/FluctuatingChargeAtomTypesSectionParser.hpp" | 
| 64 | 
> | 
#include "io/GayBerneAtomTypesSectionParser.hpp" | 
| 65 | 
> | 
#include "io/InversionTypesSectionParser.hpp" | 
| 66 | 
> | 
#include "io/LennardJonesAtomTypesSectionParser.hpp" | 
| 67 | 
> | 
#include "io/MultipoleAtomTypesSectionParser.hpp" | 
| 68 | 
> | 
#include "io/NonBondedInteractionsSectionParser.hpp" | 
| 69 | 
> | 
#include "io/PolarizableAtomTypesSectionParser.hpp" | 
| 70 | 
> | 
#include "io/SCAtomTypesSectionParser.hpp" | 
| 71 | 
> | 
#include "io/ShapeAtomTypesSectionParser.hpp" | 
| 72 | 
> | 
#include "io/StickyAtomTypesSectionParser.hpp" | 
| 73 | 
> | 
#include "io/StickyPowerAtomTypesSectionParser.hpp" | 
| 74 | 
> | 
#include "io/TorsionTypesSectionParser.hpp" | 
| 75 | 
> | 
 | 
| 76 | 
  | 
#include "types/LennardJonesAdapter.hpp" | 
| 77 | 
+ | 
#include "types/EAMAdapter.hpp" | 
| 78 | 
+ | 
#include "types/SuttonChenAdapter.hpp" | 
| 79 | 
+ | 
#include "types/GayBerneAdapter.hpp" | 
| 80 | 
+ | 
#include "types/StickyAdapter.hpp" | 
| 81 | 
  | 
 | 
| 82 | 
  | 
namespace OpenMD { | 
| 83 | 
  | 
 | 
| 84 | 
< | 
  ForceField::ForceField() {  | 
| 84 | 
> | 
  ForceField::ForceField(std::string ffName) {  | 
| 85 | 
  | 
 | 
| 86 | 
  | 
    char* tempPath;  | 
| 87 | 
  | 
    tempPath = getenv("FORCE_PARAM_PATH"); | 
| 92 | 
  | 
    } else { | 
| 93 | 
  | 
      ffPath_ = tempPath; | 
| 94 | 
  | 
    } | 
| 95 | 
+ | 
 | 
| 96 | 
+ | 
    setForceFieldFileName(ffName + ".frc"); | 
| 97 | 
+ | 
 | 
| 98 | 
+ | 
    /** | 
| 99 | 
+ | 
     * The order of adding section parsers is important. | 
| 100 | 
+ | 
     * | 
| 101 | 
+ | 
     * OptionSectionParser must come first to set options for other | 
| 102 | 
+ | 
     * parsers | 
| 103 | 
+ | 
     *  | 
| 104 | 
+ | 
     * DirectionalAtomTypesSectionParser should be added before | 
| 105 | 
+ | 
     * AtomTypesSectionParser, and these two section parsers will | 
| 106 | 
+ | 
     * actually create "real" AtomTypes (AtomTypesSectionParser will | 
| 107 | 
+ | 
     * create AtomType and DirectionalAtomTypesSectionParser will | 
| 108 | 
+ | 
     * create DirectionalAtomType, which is a subclass of AtomType and | 
| 109 | 
+ | 
     * should come first).  | 
| 110 | 
+ | 
     * | 
| 111 | 
+ | 
     * Other AtomTypes Section Parsers will not create the "real" | 
| 112 | 
+ | 
     * AtomType, they only add and set some attributes of the AtomType | 
| 113 | 
+ | 
     * (via the Adapters). Thus ordering of these is not important. | 
| 114 | 
+ | 
     * AtomTypesSectionParser should be added before other atom type | 
| 115 | 
+ | 
     * | 
| 116 | 
+ | 
     * The order of BondTypesSectionParser, BendTypesSectionParser and | 
| 117 | 
+ | 
     * TorsionTypesSectionParser, etc. are not important. | 
| 118 | 
+ | 
     */ | 
| 119 | 
+ | 
 | 
| 120 | 
+ | 
    spMan_.push_back(new OptionSectionParser(forceFieldOptions_)); | 
| 121 | 
+ | 
    spMan_.push_back(new BaseAtomTypesSectionParser()); | 
| 122 | 
+ | 
    spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_)); | 
| 123 | 
+ | 
    spMan_.push_back(new AtomTypesSectionParser()); | 
| 124 | 
+ | 
 | 
| 125 | 
+ | 
    spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_)); | 
| 126 | 
+ | 
    spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_)); | 
| 127 | 
+ | 
    spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_)); | 
| 128 | 
+ | 
    spMan_.push_back(new FluctuatingChargeAtomTypesSectionParser(forceFieldOptions_)); | 
| 129 | 
+ | 
    spMan_.push_back(new PolarizableAtomTypesSectionParser(forceFieldOptions_)); | 
| 130 | 
+ | 
    spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_)); | 
| 131 | 
+ | 
    spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_)); | 
| 132 | 
+ | 
    spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_)); | 
| 133 | 
+ | 
    spMan_.push_back(new ShapeAtomTypesSectionParser(forceFieldOptions_)); | 
| 134 | 
+ | 
    spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_)); | 
| 135 | 
+ | 
    spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_)); | 
| 136 | 
+ | 
 | 
| 137 | 
+ | 
    spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_)); | 
| 138 | 
+ | 
    spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_)); | 
| 139 | 
+ | 
    spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_)); | 
| 140 | 
+ | 
    spMan_.push_back(new InversionTypesSectionParser(forceFieldOptions_)); | 
| 141 | 
+ | 
 | 
| 142 | 
+ | 
    spMan_.push_back(new NonBondedInteractionsSectionParser(forceFieldOptions_));     | 
| 143 | 
+ | 
  } | 
| 144 | 
+ | 
 | 
| 145 | 
+ | 
  void ForceField::parse(const std::string& filename) { | 
| 146 | 
+ | 
    ifstrstream* ffStream; | 
| 147 | 
+ | 
 | 
| 148 | 
+ | 
    ffStream = openForceFieldFile(filename); | 
| 149 | 
+ | 
 | 
| 150 | 
+ | 
    spMan_.parse(*ffStream, *this); | 
| 151 | 
+ | 
 | 
| 152 | 
+ | 
    ForceField::AtomTypeContainer::MapTypeIterator i; | 
| 153 | 
+ | 
    AtomType* at; | 
| 154 | 
+ | 
 | 
| 155 | 
+ | 
    for (at = atomTypeCont_.beginType(i); at != NULL;  | 
| 156 | 
+ | 
         at = atomTypeCont_.nextType(i)) { | 
| 157 | 
+ | 
 | 
| 158 | 
+ | 
      // useBase sets the responsibilities, and these have to be done  | 
| 159 | 
+ | 
      // after the atomTypes and Base types have all been scanned: | 
| 160 | 
+ | 
 | 
| 161 | 
+ | 
      std::vector<AtomType*> ayb = at->allYourBase();       | 
| 162 | 
+ | 
      if (ayb.size() > 1) { | 
| 163 | 
+ | 
        for (int j = ayb.size()-1; j > 0; j--) { | 
| 164 | 
+ | 
           | 
| 165 | 
+ | 
          ayb[j-1]->useBase(ayb[j]); | 
| 166 | 
+ | 
 | 
| 167 | 
+ | 
        } | 
| 168 | 
+ | 
      } | 
| 169 | 
+ | 
    } | 
| 170 | 
+ | 
 | 
| 171 | 
+ | 
    delete ffStream; | 
| 172 | 
  | 
  } | 
| 173 | 
  | 
 | 
| 174 | 
  | 
  /** | 
| 739 | 
  | 
  } | 
| 740 | 
  | 
   | 
| 741 | 
  | 
  RealType ForceField::getRcutFromAtomType(AtomType* at) { | 
| 742 | 
< | 
    RealType rcut = 0.0; | 
| 742 | 
> | 
    RealType rcut(0.0); | 
| 743 | 
  | 
     | 
| 744 | 
  | 
    LennardJonesAdapter lja = LennardJonesAdapter(at); | 
| 745 | 
  | 
    if (lja.isLennardJones()) { | 
| 746 | 
  | 
      rcut = 2.5 * lja.getSigma(); | 
| 747 | 
+ | 
    } | 
| 748 | 
+ | 
    EAMAdapter ea = EAMAdapter(at); | 
| 749 | 
+ | 
    if (ea.isEAM()) { | 
| 750 | 
+ | 
      rcut = max(rcut, ea.getRcut()); | 
| 751 | 
+ | 
    } | 
| 752 | 
+ | 
    SuttonChenAdapter sca = SuttonChenAdapter(at); | 
| 753 | 
+ | 
    if (sca.isSuttonChen()) { | 
| 754 | 
+ | 
      rcut = max(rcut, 2.0 * sca.getAlpha()); | 
| 755 | 
  | 
    } | 
| 756 | 
+ | 
    GayBerneAdapter gba = GayBerneAdapter(at); | 
| 757 | 
+ | 
    if (gba.isGayBerne()) { | 
| 758 | 
+ | 
      rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL())); | 
| 759 | 
+ | 
    } | 
| 760 | 
+ | 
    StickyAdapter sa = StickyAdapter(at); | 
| 761 | 
+ | 
    if (sa.isSticky()) { | 
| 762 | 
+ | 
      rcut = max(rcut, max(sa.getRu(), sa.getRup())); | 
| 763 | 
+ | 
    } | 
| 764 | 
+ | 
 | 
| 765 | 
  | 
    return rcut;     | 
| 766 | 
  | 
  } | 
| 767 | 
  | 
   | 
| 771 | 
  | 
    ifstrstream* ffStream = new ifstrstream(); | 
| 772 | 
  | 
     | 
| 773 | 
  | 
    //try to open the force filed file in current directory first     | 
| 774 | 
< | 
    ffStream->open(forceFieldFilename.c_str()); | 
| 774 | 
> | 
    ffStream->open(forceFieldFilename.c_str(), ifstream::in | ifstream::binary); | 
| 775 | 
> | 
 | 
| 776 | 
  | 
    if(!ffStream->is_open()){ | 
| 777 | 
  | 
 | 
| 778 | 
  | 
      forceFieldFilename = ffPath_ + "/" + forceFieldFilename; | 
| 779 | 
< | 
      ffStream->open( forceFieldFilename.c_str() ); | 
| 779 | 
> | 
      ffStream->open( forceFieldFilename.c_str(),   | 
| 780 | 
> | 
                      ifstream::in | ifstream::binary ); | 
| 781 | 
  | 
 | 
| 782 | 
  | 
      //if current directory does not contain the force field file, | 
| 783 | 
  | 
      //try to open it in the path         |