| 1 | < | /* | 
| 1 | > | /* | 
| 2 |  | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. | 
| 3 |  | * | 
| 4 |  | * The University of Notre Dame grants you ("Licensee") a | 
| 6 |  | * redistribute this software in source and binary code form, provided | 
| 7 |  | * that the following conditions are met: | 
| 8 |  | * | 
| 9 | < | * 1. Acknowledgement of the program authors must be made in any | 
| 10 | < | *    publication of scientific results based in part on use of the | 
| 11 | < | *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | < | *    the article in which the program was described (Matthew | 
| 13 | < | *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | < | *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | < | *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | < | *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | < | * | 
| 18 | < | * 2. Redistributions of source code must retain the above copyright | 
| 9 | > | * 1. Redistributions of source code must retain the above copyright | 
| 10 |  | *    notice, this list of conditions and the following disclaimer. | 
| 11 |  | * | 
| 12 | < | * 3. Redistributions in binary form must reproduce the above copyright | 
| 12 | > | * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 |  | *    notice, this list of conditions and the following disclaimer in the | 
| 14 |  | *    documentation and/or other materials provided with the | 
| 15 |  | *    distribution. | 
| 28 |  | * arising out of the use of or inability to use software, even if the | 
| 29 |  | * University of Notre Dame has been advised of the possibility of | 
| 30 |  | * such damages. | 
| 31 | + | * | 
| 32 | + | * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your | 
| 33 | + | * research, please cite the appropriate papers when you publish your | 
| 34 | + | * work.  Good starting points are: | 
| 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). | 
| 39 | + | * [4]  Vardeman & Gezelter, in progress (2009). | 
| 40 |  | */ | 
| 41 |  |  | 
| 42 | < | /** | 
| 43 | < | * @file ForceField.cpp | 
| 44 | < | * @author tlin | 
| 45 | < | * @date 11/04/2004 | 
| 46 | < | * @time 22:51am | 
| 47 | < | * @version 1.0 | 
| 48 | < | */ | 
| 42 | > | /** | 
| 43 | > | * @file ForceField.cpp | 
| 44 | > | * @author tlin | 
| 45 | > | * @date 11/04/2004 | 
| 46 | > | * @time 22:51am | 
| 47 | > | * @version 1.0 | 
| 48 | > | */ | 
| 49 |  |  | 
| 50 | + | #include <algorithm> | 
| 51 |  | #include "UseTheForce/ForceField.hpp" | 
| 52 |  | #include "utils/simError.h" | 
| 53 | < | #include "UseTheForce/DarkSide/atype_interface.h" | 
| 54 | < | namespace oopse { | 
| 53 | > | #include "utils/Tuple.hpp" | 
| 54 | > | namespace OpenMD { | 
| 55 |  |  | 
| 56 | < | ForceField::ForceField() { | 
| 56 | > | ForceField::ForceField() { | 
| 57 | > |  | 
| 58 |  | char* tempPath; | 
| 59 |  | tempPath = getenv("FORCE_PARAM_PATH"); | 
| 60 | < |  | 
| 60 | > |  | 
| 61 |  | if (tempPath == NULL) { | 
| 62 | < | //convert a macro from compiler to a string in c++ | 
| 63 | < | STR_DEFINE(ffPath_, FRC_PATH ); | 
| 62 | > | //convert a macro from compiler to a string in c++ | 
| 63 | > | STR_DEFINE(ffPath_, FRC_PATH ); | 
| 64 |  | } else { | 
| 65 | < | ffPath_ = tempPath; | 
| 65 | > | ffPath_ = tempPath; | 
| 66 |  | } | 
| 67 | < | } | 
| 67 | > | } | 
| 68 |  |  | 
| 69 | < |  | 
| 70 | < | ForceField::~ForceField() { | 
| 71 | < | deleteAtypes(); | 
| 72 | < | } | 
| 73 | < |  | 
| 74 | < | AtomType* ForceField::getAtomType(const std::string &at) { | 
| 69 | > | /** | 
| 70 | > | * getAtomType by string | 
| 71 | > | * | 
| 72 | > | * finds the requested atom type in this force field using the string | 
| 73 | > | * name of the atom type. | 
| 74 | > | */ | 
| 75 | > | AtomType* ForceField::getAtomType(const std::string &at) { | 
| 76 |  | std::vector<std::string> keys; | 
| 77 |  | keys.push_back(at); | 
| 78 |  | return atomTypeCont_.find(keys); | 
| 79 | < | } | 
| 79 | > | } | 
| 80 |  |  | 
| 81 | < | BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) { | 
| 81 | > | /** | 
| 82 | > | * getAtomType by ident | 
| 83 | > | * | 
| 84 | > | * finds the requested atom type in this force field using the | 
| 85 | > | * integer ident instead of the string name of the atom type. | 
| 86 | > | */ | 
| 87 | > | AtomType* ForceField::getAtomType(int ident) { | 
| 88 | > | std::string at = atypeIdentToName.find(ident)->second; | 
| 89 | > | return getAtomType(at); | 
| 90 | > | } | 
| 91 | > |  | 
| 92 | > | BondType* ForceField::getBondType(const std::string &at1, | 
| 93 | > | const std::string &at2) { | 
| 94 |  | std::vector<std::string> keys; | 
| 95 |  | keys.push_back(at1); | 
| 96 |  | keys.push_back(at2); | 
| 98 |  | //try exact match first | 
| 99 |  | BondType* bondType = bondTypeCont_.find(keys); | 
| 100 |  | if (bondType) { | 
| 101 | < | return bondType; | 
| 101 | > | return bondType; | 
| 102 |  | } else { | 
| 103 | < | //if no exact match found, try wild card match | 
| 104 | < | return bondTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 105 | < | } | 
| 103 | > | AtomType* atype1; | 
| 104 | > | AtomType* atype2; | 
| 105 | > | std::vector<std::string> at1key; | 
| 106 | > | at1key.push_back(at1); | 
| 107 | > | atype1 = atomTypeCont_.find(at1key); | 
| 108 | > |  | 
| 109 | > | std::vector<std::string> at2key; | 
| 110 | > | at2key.push_back(at2); | 
| 111 | > | atype2 = atomTypeCont_.find(at2key); | 
| 112 |  |  | 
| 113 | < | } | 
| 113 | > | // query atom types for their chains of responsibility | 
| 114 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 115 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 116 |  |  | 
| 117 | < | BendType* ForceField::getBendType(const std::string &at1, const std::string &at2, | 
| 118 | < | const std::string &at3) { | 
| 117 | > | std::vector<AtomType*>::iterator i; | 
| 118 | > | std::vector<AtomType*>::iterator j; | 
| 119 | > |  | 
| 120 | > | int ii = 0; | 
| 121 | > | int jj = 0; | 
| 122 | > | int bondTypeScore; | 
| 123 | > |  | 
| 124 | > | std::vector<std::pair<int, std::vector<std::string> > > foundBonds; | 
| 125 | > |  | 
| 126 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 127 | > | jj = 0; | 
| 128 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 129 | > |  | 
| 130 | > | bondTypeScore = ii + jj; | 
| 131 | > |  | 
| 132 | > | std::vector<std::string> myKeys; | 
| 133 | > | myKeys.push_back((*i)->getName()); | 
| 134 | > | myKeys.push_back((*j)->getName()); | 
| 135 | > |  | 
| 136 | > | BondType* bondType = bondTypeCont_.find(myKeys); | 
| 137 | > | if (bondType) { | 
| 138 | > | foundBonds.push_back(std::make_pair(bondTypeScore, myKeys)); | 
| 139 | > | } | 
| 140 | > | jj++; | 
| 141 | > | } | 
| 142 | > | ii++; | 
| 143 | > | } | 
| 144 | > |  | 
| 145 | > |  | 
| 146 | > | if (foundBonds.size() > 0) { | 
| 147 | > | // sort the foundBonds by the score: | 
| 148 | > | std::sort(foundBonds.begin(), foundBonds.end()); | 
| 149 | > |  | 
| 150 | > | int bestScore = foundBonds[0].first; | 
| 151 | > | std::vector<std::string> theKeys = foundBonds[0].second; | 
| 152 | > |  | 
| 153 | > | BondType* bestType = bondTypeCont_.find(theKeys); | 
| 154 | > |  | 
| 155 | > | return bestType; | 
| 156 | > | } else { | 
| 157 | > | //if no exact match found, try wild card match | 
| 158 | > | return bondTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 159 | > | } | 
| 160 | > | } | 
| 161 | > | } | 
| 162 | > |  | 
| 163 | > | BendType* ForceField::getBendType(const std::string &at1, | 
| 164 | > | const std::string &at2, | 
| 165 | > | const std::string &at3) { | 
| 166 |  | std::vector<std::string> keys; | 
| 167 |  | keys.push_back(at1); | 
| 168 |  | keys.push_back(at2); | 
| 171 |  | //try exact match first | 
| 172 |  | BendType* bendType = bendTypeCont_.find(keys); | 
| 173 |  | if (bendType) { | 
| 174 | < | return bendType; | 
| 174 | > | return bendType; | 
| 175 |  | } else { | 
| 176 | < | //if no exact match found, try wild card match | 
| 177 | < | return bendTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 176 | > |  | 
| 177 | > | AtomType* atype1; | 
| 178 | > | AtomType* atype2; | 
| 179 | > | AtomType* atype3; | 
| 180 | > | std::vector<std::string> at1key; | 
| 181 | > | at1key.push_back(at1); | 
| 182 | > | atype1 = atomTypeCont_.find(at1key); | 
| 183 | > |  | 
| 184 | > | std::vector<std::string> at2key; | 
| 185 | > | at2key.push_back(at2); | 
| 186 | > | atype2 = atomTypeCont_.find(at2key); | 
| 187 | > |  | 
| 188 | > | std::vector<std::string> at3key; | 
| 189 | > | at3key.push_back(at3); | 
| 190 | > | atype3 = atomTypeCont_.find(at3key); | 
| 191 | > |  | 
| 192 | > | // query atom types for their chains of responsibility | 
| 193 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 194 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 195 | > | std::vector<AtomType*> at3Chain = atype3->allYourBase(); | 
| 196 | > |  | 
| 197 | > | std::vector<AtomType*>::iterator i; | 
| 198 | > | std::vector<AtomType*>::iterator j; | 
| 199 | > | std::vector<AtomType*>::iterator k; | 
| 200 | > |  | 
| 201 | > | int ii = 0; | 
| 202 | > | int jj = 0; | 
| 203 | > | int kk = 0; | 
| 204 | > | int IKscore; | 
| 205 | > |  | 
| 206 | > | std::vector<tuple3<int, int, std::vector<std::string> > > foundBends; | 
| 207 | > |  | 
| 208 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 209 | > | ii = 0; | 
| 210 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 211 | > | kk = 0; | 
| 212 | > | for (k = at3Chain.begin(); k != at3Chain.end(); k++) { | 
| 213 | > |  | 
| 214 | > | IKscore = ii + kk; | 
| 215 | > |  | 
| 216 | > | std::vector<std::string> myKeys; | 
| 217 | > | myKeys.push_back((*i)->getName()); | 
| 218 | > | myKeys.push_back((*j)->getName()); | 
| 219 | > | myKeys.push_back((*k)->getName()); | 
| 220 | > |  | 
| 221 | > | BendType* bendType = bendTypeCont_.find(myKeys); | 
| 222 | > | if (bendType) { | 
| 223 | > | foundBends.push_back( make_tuple3(jj, IKscore, myKeys) ); | 
| 224 | > | } | 
| 225 | > | kk++; | 
| 226 | > | } | 
| 227 | > | ii++; | 
| 228 | > | } | 
| 229 | > | jj++; | 
| 230 | > | } | 
| 231 | > |  | 
| 232 | > | if (foundBends.size() > 0) { | 
| 233 | > | std::sort(foundBends.begin(), foundBends.end()); | 
| 234 | > | int jscore = foundBends[0].first; | 
| 235 | > | int ikscore = foundBends[0].second; | 
| 236 | > | std::vector<std::string> theKeys = foundBends[0].third; | 
| 237 | > |  | 
| 238 | > | BendType* bestType = bendTypeCont_.find(theKeys); | 
| 239 | > | return bestType; | 
| 240 | > | } else { | 
| 241 | > | //if no exact match found, try wild card match | 
| 242 | > | return bendTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 243 | > | } | 
| 244 |  | } | 
| 245 | < | } | 
| 245 | > | } | 
| 246 |  |  | 
| 247 | < | TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2, | 
| 248 | < | const std::string &at3, const std::string &at4) { | 
| 247 | > | TorsionType* ForceField::getTorsionType(const std::string &at1, | 
| 248 | > | const std::string &at2, | 
| 249 | > | const std::string &at3, | 
| 250 | > | const std::string &at4) { | 
| 251 |  | std::vector<std::string> keys; | 
| 252 |  | keys.push_back(at1); | 
| 253 |  | keys.push_back(at2); | 
| 254 |  | keys.push_back(at3); | 
| 255 |  | keys.push_back(at4); | 
| 256 |  |  | 
| 257 | + |  | 
| 258 | + | //try exact match first | 
| 259 |  | TorsionType* torsionType = torsionTypeCont_.find(keys); | 
| 260 |  | if (torsionType) { | 
| 261 | < | return torsionType; | 
| 261 | > | return torsionType; | 
| 262 |  | } else { | 
| 123 | – | //if no exact match found, try wild card match | 
| 124 | – | return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 125 | – | } | 
| 126 | – |  | 
| 127 | – | return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 263 |  |  | 
| 264 | < | } | 
| 264 | > | AtomType* atype1; | 
| 265 | > | AtomType* atype2; | 
| 266 | > | AtomType* atype3; | 
| 267 | > | AtomType* atype4; | 
| 268 | > | std::vector<std::string> at1key; | 
| 269 | > | at1key.push_back(at1); | 
| 270 | > | atype1 = atomTypeCont_.find(at1key); | 
| 271 | > |  | 
| 272 | > | std::vector<std::string> at2key; | 
| 273 | > | at2key.push_back(at2); | 
| 274 | > | atype2 = atomTypeCont_.find(at2key); | 
| 275 |  |  | 
| 276 | < | BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){ | 
| 276 | > | std::vector<std::string> at3key; | 
| 277 | > | at3key.push_back(at3); | 
| 278 | > | atype3 = atomTypeCont_.find(at3key); | 
| 279 | > |  | 
| 280 | > | std::vector<std::string> at4key; | 
| 281 | > | at4key.push_back(at4); | 
| 282 | > | atype4 = atomTypeCont_.find(at4key); | 
| 283 | > |  | 
| 284 | > | // query atom types for their chains of responsibility | 
| 285 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 286 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 287 | > | std::vector<AtomType*> at3Chain = atype3->allYourBase(); | 
| 288 | > | std::vector<AtomType*> at4Chain = atype4->allYourBase(); | 
| 289 | > |  | 
| 290 | > | std::vector<AtomType*>::iterator i; | 
| 291 | > | std::vector<AtomType*>::iterator j; | 
| 292 | > | std::vector<AtomType*>::iterator k; | 
| 293 | > | std::vector<AtomType*>::iterator l; | 
| 294 | > |  | 
| 295 | > | int ii = 0; | 
| 296 | > | int jj = 0; | 
| 297 | > | int kk = 0; | 
| 298 | > | int ll = 0; | 
| 299 | > | int ILscore; | 
| 300 | > | int JKscore; | 
| 301 | > |  | 
| 302 | > | std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions; | 
| 303 | > |  | 
| 304 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 305 | > | kk = 0; | 
| 306 | > | for (k = at3Chain.begin(); k != at3Chain.end(); k++) { | 
| 307 | > | ii = 0; | 
| 308 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 309 | > | ll = 0; | 
| 310 | > | for (l = at4Chain.begin(); l != at4Chain.end(); l++) { | 
| 311 | > |  | 
| 312 | > | ILscore = ii + ll; | 
| 313 | > | JKscore = jj + kk; | 
| 314 | > |  | 
| 315 | > | std::vector<std::string> myKeys; | 
| 316 | > | myKeys.push_back((*i)->getName()); | 
| 317 | > | myKeys.push_back((*j)->getName()); | 
| 318 | > | myKeys.push_back((*k)->getName()); | 
| 319 | > | myKeys.push_back((*l)->getName()); | 
| 320 | > |  | 
| 321 | > | TorsionType* torsionType = torsionTypeCont_.find(myKeys); | 
| 322 | > | if (torsionType) { | 
| 323 | > | foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) ); | 
| 324 | > | } | 
| 325 | > | ll++; | 
| 326 | > | } | 
| 327 | > | ii++; | 
| 328 | > | } | 
| 329 | > | kk++; | 
| 330 | > | } | 
| 331 | > | jj++; | 
| 332 | > | } | 
| 333 | > |  | 
| 334 | > | if (foundTorsions.size() > 0) { | 
| 335 | > | std::sort(foundTorsions.begin(), foundTorsions.end()); | 
| 336 | > | int jkscore = foundTorsions[0].first; | 
| 337 | > | int ilscore = foundTorsions[0].second; | 
| 338 | > | std::vector<std::string> theKeys = foundTorsions[0].third; | 
| 339 | > |  | 
| 340 | > | TorsionType* bestType = torsionTypeCont_.find(theKeys); | 
| 341 | > | return bestType; | 
| 342 | > | } else { | 
| 343 | > | //if no exact match found, try wild card match | 
| 344 | > | return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 345 | > | } | 
| 346 | > | } | 
| 347 | > | } | 
| 348 | > |  | 
| 349 | > | InversionType* ForceField::getInversionType(const std::string &at1, | 
| 350 | > | const std::string &at2, | 
| 351 | > | const std::string &at3, | 
| 352 | > | const std::string &at4) { | 
| 353 |  | std::vector<std::string> keys; | 
| 354 |  | keys.push_back(at1); | 
| 355 |  | keys.push_back(at2); | 
| 356 | < | return bondTypeCont_.find(keys); | 
| 357 | < | } | 
| 356 | > | keys.push_back(at3); | 
| 357 | > | keys.push_back(at4); | 
| 358 |  |  | 
| 359 | < | BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2, | 
| 360 | < | const std::string &at3){ | 
| 359 | > | //try exact match first | 
| 360 | > | InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys); | 
| 361 | > | if (inversionType) { | 
| 362 | > | return inversionType; | 
| 363 | > | } else { | 
| 364 | > |  | 
| 365 | > | AtomType* atype1; | 
| 366 | > | AtomType* atype2; | 
| 367 | > | AtomType* atype3; | 
| 368 | > | AtomType* atype4; | 
| 369 | > | std::vector<std::string> at1key; | 
| 370 | > | at1key.push_back(at1); | 
| 371 | > | atype1 = atomTypeCont_.find(at1key); | 
| 372 | > |  | 
| 373 | > | std::vector<std::string> at2key; | 
| 374 | > | at2key.push_back(at2); | 
| 375 | > | atype2 = atomTypeCont_.find(at2key); | 
| 376 | > |  | 
| 377 | > | std::vector<std::string> at3key; | 
| 378 | > | at3key.push_back(at3); | 
| 379 | > | atype3 = atomTypeCont_.find(at3key); | 
| 380 | > |  | 
| 381 | > | std::vector<std::string> at4key; | 
| 382 | > | at4key.push_back(at4); | 
| 383 | > | atype4 = atomTypeCont_.find(at4key); | 
| 384 | > |  | 
| 385 | > | // query atom types for their chains of responsibility | 
| 386 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 387 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 388 | > | std::vector<AtomType*> at3Chain = atype3->allYourBase(); | 
| 389 | > | std::vector<AtomType*> at4Chain = atype4->allYourBase(); | 
| 390 | > |  | 
| 391 | > | std::vector<AtomType*>::iterator i; | 
| 392 | > | std::vector<AtomType*>::iterator j; | 
| 393 | > | std::vector<AtomType*>::iterator k; | 
| 394 | > | std::vector<AtomType*>::iterator l; | 
| 395 | > |  | 
| 396 | > | int ii = 0; | 
| 397 | > | int jj = 0; | 
| 398 | > | int kk = 0; | 
| 399 | > | int ll = 0; | 
| 400 | > | int Iscore; | 
| 401 | > | int JKLscore; | 
| 402 | > |  | 
| 403 | > | std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions; | 
| 404 | > |  | 
| 405 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 406 | > | kk = 0; | 
| 407 | > | for (k = at3Chain.begin(); k != at3Chain.end(); k++) { | 
| 408 | > | ii = 0; | 
| 409 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 410 | > | ll = 0; | 
| 411 | > | for (l = at4Chain.begin(); l != at4Chain.end(); l++) { | 
| 412 | > |  | 
| 413 | > | Iscore = ii; | 
| 414 | > | JKLscore = jj + kk + ll; | 
| 415 | > |  | 
| 416 | > | std::vector<std::string> myKeys; | 
| 417 | > | myKeys.push_back((*i)->getName()); | 
| 418 | > | myKeys.push_back((*j)->getName()); | 
| 419 | > | myKeys.push_back((*k)->getName()); | 
| 420 | > | myKeys.push_back((*l)->getName()); | 
| 421 | > |  | 
| 422 | > | InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys); | 
| 423 | > | if (inversionType) { | 
| 424 | > | foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) ); | 
| 425 | > | } | 
| 426 | > | ll++; | 
| 427 | > | } | 
| 428 | > | ii++; | 
| 429 | > | } | 
| 430 | > | kk++; | 
| 431 | > | } | 
| 432 | > | jj++; | 
| 433 | > | } | 
| 434 | > |  | 
| 435 | > | if (foundInversions.size() > 0) { | 
| 436 | > | std::sort(foundInversions.begin(), foundInversions.end()); | 
| 437 | > | int iscore = foundInversions[0].first; | 
| 438 | > | int jklscore = foundInversions[0].second; | 
| 439 | > | std::vector<std::string> theKeys = foundInversions[0].third; | 
| 440 | > |  | 
| 441 | > | InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys); | 
| 442 | > | return bestType; | 
| 443 | > | } else { | 
| 444 | > | //if no exact match found, try wild card match | 
| 445 | > | return inversionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 446 | > | } | 
| 447 | > | } | 
| 448 | > | } | 
| 449 | > |  | 
| 450 | > | NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { | 
| 451 |  | std::vector<std::string> keys; | 
| 452 |  | keys.push_back(at1); | 
| 453 |  | keys.push_back(at2); | 
| 454 | + |  | 
| 455 | + | //try exact match first | 
| 456 | + | NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys); | 
| 457 | + | if (nbiType) { | 
| 458 | + | return nbiType; | 
| 459 | + | } else { | 
| 460 | + | //if no exact match found, try wild card match | 
| 461 | + | return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 462 | + | } | 
| 463 | + | } | 
| 464 | + |  | 
| 465 | + | BondType* ForceField::getExactBondType(const std::string &at1, | 
| 466 | + | const std::string &at2){ | 
| 467 | + | std::vector<std::string> keys; | 
| 468 | + | keys.push_back(at1); | 
| 469 | + | keys.push_back(at2); | 
| 470 | + | return bondTypeCont_.find(keys); | 
| 471 | + | } | 
| 472 | + |  | 
| 473 | + | BendType* ForceField::getExactBendType(const std::string &at1, | 
| 474 | + | const std::string &at2, | 
| 475 | + | const std::string &at3){ | 
| 476 | + | std::vector<std::string> keys; | 
| 477 | + | keys.push_back(at1); | 
| 478 | + | keys.push_back(at2); | 
| 479 |  | keys.push_back(at3); | 
| 480 |  | return bendTypeCont_.find(keys); | 
| 481 | < | } | 
| 482 | < |  | 
| 483 | < | TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2, | 
| 484 | < | const std::string &at3, const std::string &at4){ | 
| 481 | > | } | 
| 482 | > |  | 
| 483 | > | TorsionType* ForceField::getExactTorsionType(const std::string &at1, | 
| 484 | > | const std::string &at2, | 
| 485 | > | const std::string &at3, | 
| 486 | > | const std::string &at4){ | 
| 487 |  | std::vector<std::string> keys; | 
| 488 |  | keys.push_back(at1); | 
| 489 |  | keys.push_back(at2); | 
| 490 |  | keys.push_back(at3); | 
| 491 |  | keys.push_back(at4); | 
| 492 |  | return torsionTypeCont_.find(keys); | 
| 493 | < | } | 
| 494 | < | bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { | 
| 493 | > | } | 
| 494 | > |  | 
| 495 | > | InversionType* ForceField::getExactInversionType(const std::string &at1, | 
| 496 | > | const std::string &at2, | 
| 497 | > | const std::string &at3, | 
| 498 | > | const std::string &at4){ | 
| 499 |  | std::vector<std::string> keys; | 
| 500 | + | keys.push_back(at1); | 
| 501 | + | keys.push_back(at2); | 
| 502 | + | keys.push_back(at3); | 
| 503 | + | keys.push_back(at4); | 
| 504 | + | return inversionTypeCont_.find(keys); | 
| 505 | + | } | 
| 506 | + |  | 
| 507 | + | NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ | 
| 508 | + | std::vector<std::string> keys; | 
| 509 | + | keys.push_back(at1); | 
| 510 | + | keys.push_back(at2); | 
| 511 | + | return nonBondedInteractionTypeCont_.find(keys); | 
| 512 | + | } | 
| 513 | + |  | 
| 514 | + |  | 
| 515 | + | bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { | 
| 516 | + | std::vector<std::string> keys; | 
| 517 |  | keys.push_back(at); | 
| 518 | + | atypeIdentToName[atomType->getIdent()] = at; | 
| 519 |  | return atomTypeCont_.add(keys, atomType); | 
| 520 | < | } | 
| 520 | > | } | 
| 521 |  |  | 
| 522 | < | bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) { | 
| 522 | > | bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { | 
| 523 |  | std::vector<std::string> keys; | 
| 524 | + | keys.push_back(at); | 
| 525 | + | atypeIdentToName[atomType->getIdent()] = at; | 
| 526 | + | return atomTypeCont_.replace(keys, atomType); | 
| 527 | + | } | 
| 528 | + |  | 
| 529 | + | bool ForceField::addBondType(const std::string &at1, const std::string &at2, | 
| 530 | + | BondType* bondType) { | 
| 531 | + | std::vector<std::string> keys; | 
| 532 |  | keys.push_back(at1); | 
| 533 |  | keys.push_back(at2); | 
| 534 | < | return bondTypeCont_.add(keys, bondType); | 
| 535 | < |  | 
| 536 | < | } | 
| 537 | < |  | 
| 538 | < | bool ForceField::addBendType(const std::string &at1, const std::string &at2, | 
| 171 | < | const std::string &at3, BendType* bendType) { | 
| 534 | > | return bondTypeCont_.add(keys, bondType); | 
| 535 | > | } | 
| 536 | > |  | 
| 537 | > | bool ForceField::addBendType(const std::string &at1, const std::string &at2, | 
| 538 | > | const std::string &at3, BendType* bendType) { | 
| 539 |  | std::vector<std::string> keys; | 
| 540 |  | keys.push_back(at1); | 
| 541 |  | keys.push_back(at2); | 
| 542 |  | keys.push_back(at3); | 
| 543 |  | return bendTypeCont_.add(keys, bendType); | 
| 544 | < | } | 
| 545 | < |  | 
| 546 | < | bool ForceField::addTorsionType(const std::string &at1, const std::string &at2, | 
| 547 | < | const std::string &at3, const std::string &at4, TorsionType* torsionType) { | 
| 544 | > | } | 
| 545 | > |  | 
| 546 | > | bool ForceField::addTorsionType(const std::string &at1, | 
| 547 | > | const std::string &at2, | 
| 548 | > | const std::string &at3, | 
| 549 | > | const std::string &at4, | 
| 550 | > | TorsionType* torsionType) { | 
| 551 |  | std::vector<std::string> keys; | 
| 552 |  | keys.push_back(at1); | 
| 553 |  | keys.push_back(at2); | 
| 554 |  | keys.push_back(at3); | 
| 555 |  | keys.push_back(at4); | 
| 556 |  | return torsionTypeCont_.add(keys, torsionType); | 
| 557 | < | } | 
| 557 | > | } | 
| 558 |  |  | 
| 559 | < | double ForceField::getRcutFromAtomType(AtomType* at) { | 
| 559 | > | bool ForceField::addInversionType(const std::string &at1, | 
| 560 | > | const std::string &at2, | 
| 561 | > | const std::string &at3, | 
| 562 | > | const std::string &at4, | 
| 563 | > | InversionType* inversionType) { | 
| 564 | > | std::vector<std::string> keys; | 
| 565 | > | keys.push_back(at1); | 
| 566 | > | keys.push_back(at2); | 
| 567 | > | keys.push_back(at3); | 
| 568 | > | keys.push_back(at4); | 
| 569 | > | return inversionTypeCont_.add(keys, inversionType); | 
| 570 | > | } | 
| 571 | > |  | 
| 572 | > | bool ForceField::addNonBondedInteractionType(const std::string &at1, | 
| 573 | > | const std::string &at2, | 
| 574 | > | NonBondedInteractionType* nbiType) { | 
| 575 | > | std::vector<std::string> keys; | 
| 576 | > | keys.push_back(at1); | 
| 577 | > | keys.push_back(at2); | 
| 578 | > | return nonBondedInteractionTypeCont_.add(keys, nbiType); | 
| 579 | > | } | 
| 580 | > |  | 
| 581 | > | RealType ForceField::getRcutFromAtomType(AtomType* at) { | 
| 582 |  | /**@todo */ | 
| 583 |  | GenericData* data; | 
| 584 | < | double rcut = 0.0; | 
| 585 | < |  | 
| 584 | > | RealType rcut = 0.0; | 
| 585 | > |  | 
| 586 |  | if (at->isLennardJones()) { | 
| 587 | < | data = at->getPropertyByName("LennardJones"); | 
| 588 | < | if (data != NULL) { | 
| 589 | < | LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); | 
| 590 | < |  | 
| 591 | < | if (ljData != NULL) { | 
| 592 | < | LJParam ljParam = ljData->getData(); | 
| 593 | < |  | 
| 594 | < | //by default use 2.5*sigma as cutoff radius | 
| 595 | < | rcut = 2.5 * ljParam.sigma; | 
| 596 | < |  | 
| 597 | < | } else { | 
| 598 | < | sprintf( painCave.errMsg, | 
| 599 | < | "Can not cast GenericData to LJParam\n"); | 
| 600 | < | painCave.severity = OOPSE_ERROR; | 
| 601 | < | painCave.isFatal = 1; | 
| 602 | < | simError(); | 
| 603 | < | } | 
| 604 | < | } else { | 
| 605 | < | sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); | 
| 606 | < | painCave.severity = OOPSE_ERROR; | 
| 607 | < | painCave.isFatal = 1; | 
| 608 | < | simError(); | 
| 609 | < | } | 
| 587 | > | data = at->getPropertyByName("LennardJones"); | 
| 588 | > | if (data != NULL) { | 
| 589 | > | LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); | 
| 590 | > |  | 
| 591 | > | if (ljData != NULL) { | 
| 592 | > | LJParam ljParam = ljData->getData(); | 
| 593 | > |  | 
| 594 | > | //by default use 2.5*sigma as cutoff radius | 
| 595 | > | rcut = 2.5 * ljParam.sigma; | 
| 596 | > |  | 
| 597 | > | } else { | 
| 598 | > | sprintf( painCave.errMsg, | 
| 599 | > | "Can not cast GenericData to LJParam\n"); | 
| 600 | > | painCave.severity = OPENMD_ERROR; | 
| 601 | > | painCave.isFatal = 1; | 
| 602 | > | simError(); | 
| 603 | > | } | 
| 604 | > | } else { | 
| 605 | > | sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); | 
| 606 | > | painCave.severity = OPENMD_ERROR; | 
| 607 | > | painCave.isFatal = 1; | 
| 608 | > | simError(); | 
| 609 | > | } | 
| 610 |  | } | 
| 219 | – |  | 
| 611 |  | return rcut; | 
| 612 | < | } | 
| 612 | > | } | 
| 613 | > |  | 
| 614 |  |  | 
| 615 | < |  | 
| 224 | < | ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { | 
| 615 | > | ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { | 
| 616 |  | std::string forceFieldFilename(filename); | 
| 617 |  | ifstrstream* ffStream = new ifstrstream(); | 
| 618 |  |  | 
| 620 |  | ffStream->open(forceFieldFilename.c_str()); | 
| 621 |  | if(!ffStream->is_open()){ | 
| 622 |  |  | 
| 623 | < | forceFieldFilename = ffPath_ + "/" + forceFieldFilename; | 
| 624 | < | ffStream->open( forceFieldFilename.c_str() ); | 
| 623 | > | forceFieldFilename = ffPath_ + "/" + forceFieldFilename; | 
| 624 | > | ffStream->open( forceFieldFilename.c_str() ); | 
| 625 |  |  | 
| 626 | < | //if current directory does not contain the force field file, | 
| 627 | < | //try to open it in the path | 
| 628 | < | if(!ffStream->is_open()){ | 
| 626 | > | //if current directory does not contain the force field file, | 
| 627 | > | //try to open it in the path | 
| 628 | > | if(!ffStream->is_open()){ | 
| 629 |  |  | 
| 630 | < | sprintf( painCave.errMsg, | 
| 631 | < | "Error opening the force field parameter file:\n" | 
| 632 | < | "\t%s\n" | 
| 633 | < | "\tHave you tried setting the FORCE_PARAM_PATH environment " | 
| 634 | < | "variable?\n", | 
| 635 | < | forceFieldFilename.c_str() ); | 
| 636 | < | painCave.severity = OOPSE_ERROR; | 
| 637 | < | painCave.isFatal = 1; | 
| 638 | < | simError(); | 
| 639 | < | } | 
| 630 | > | sprintf( painCave.errMsg, | 
| 631 | > | "Error opening the force field parameter file:\n" | 
| 632 | > | "\t%s\n" | 
| 633 | > | "\tHave you tried setting the FORCE_PARAM_PATH environment " | 
| 634 | > | "variable?\n", | 
| 635 | > | forceFieldFilename.c_str() ); | 
| 636 | > | painCave.severity = OPENMD_ERROR; | 
| 637 | > | painCave.isFatal = 1; | 
| 638 | > | simError(); | 
| 639 | > | } | 
| 640 |  | } | 
| 250 | – |  | 
| 641 |  | return ffStream; | 
| 642 | + | } | 
| 643 |  |  | 
| 644 | < | } | 
| 254 | < |  | 
| 255 | < | } //end namespace oopse | 
| 644 | > | } //end namespace OpenMD |