| 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 |  | /** | 
| 175 |  | * getAtomType by string | 
| 176 |  | * | 
| 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 |  |  |