| 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | + | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 |  | */ | 
| 42 |  |  | 
| 43 |  | /** | 
| 48 |  | * @version 1.0 | 
| 49 |  | */ | 
| 50 |  |  | 
| 51 | + | #include <algorithm> | 
| 52 |  | #include "UseTheForce/ForceField.hpp" | 
| 53 |  | #include "utils/simError.h" | 
| 54 | < | #include "UseTheForce/DarkSide/atype_interface.h" | 
| 55 | < | namespace oopse { | 
| 54 | > | #include "utils/Tuple.hpp" | 
| 55 | > | #include "types/LennardJonesAdapter.hpp" | 
| 56 |  |  | 
| 57 | + | namespace OpenMD { | 
| 58 | + |  | 
| 59 |  | ForceField::ForceField() { | 
| 60 | + |  | 
| 61 |  | char* tempPath; | 
| 62 |  | tempPath = getenv("FORCE_PARAM_PATH"); | 
| 63 | < |  | 
| 63 | > |  | 
| 64 |  | if (tempPath == NULL) { | 
| 65 |  | //convert a macro from compiler to a string in c++ | 
| 66 |  | STR_DEFINE(ffPath_, FRC_PATH ); | 
| 69 |  | } | 
| 70 |  | } | 
| 71 |  |  | 
| 72 | < |  | 
| 73 | < | ForceField::~ForceField() { | 
| 74 | < | deleteAtypes(); | 
| 75 | < | } | 
| 76 | < |  | 
| 72 | > | /** | 
| 73 | > | * getAtomType by string | 
| 74 | > | * | 
| 75 | > | * finds the requested atom type in this force field using the string | 
| 76 | > | * name of the atom type. | 
| 77 | > | */ | 
| 78 |  | AtomType* ForceField::getAtomType(const std::string &at) { | 
| 79 |  | std::vector<std::string> keys; | 
| 80 |  | keys.push_back(at); | 
| 81 |  | return atomTypeCont_.find(keys); | 
| 82 |  | } | 
| 83 |  |  | 
| 84 | < | BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) { | 
| 84 | > | /** | 
| 85 | > | * getAtomType by ident | 
| 86 | > | * | 
| 87 | > | * finds the requested atom type in this force field using the | 
| 88 | > | * integer ident instead of the string name of the atom type. | 
| 89 | > | */ | 
| 90 | > | AtomType* ForceField::getAtomType(int ident) { | 
| 91 | > | std::string at = atypeIdentToName.find(ident)->second; | 
| 92 | > | return getAtomType(at); | 
| 93 | > | } | 
| 94 | > |  | 
| 95 | > | BondType* ForceField::getBondType(const std::string &at1, | 
| 96 | > | const std::string &at2) { | 
| 97 |  | std::vector<std::string> keys; | 
| 98 |  | keys.push_back(at1); | 
| 99 |  | keys.push_back(at2); | 
| 103 |  | if (bondType) { | 
| 104 |  | return bondType; | 
| 105 |  | } else { | 
| 106 | < | //if no exact match found, try wild card match | 
| 107 | < | return bondTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 108 | < | } | 
| 106 | > | AtomType* atype1; | 
| 107 | > | AtomType* atype2; | 
| 108 | > | std::vector<std::string> at1key; | 
| 109 | > | at1key.push_back(at1); | 
| 110 | > | atype1 = atomTypeCont_.find(at1key); | 
| 111 | > |  | 
| 112 | > | std::vector<std::string> at2key; | 
| 113 | > | at2key.push_back(at2); | 
| 114 | > | atype2 = atomTypeCont_.find(at2key); | 
| 115 |  |  | 
| 116 | < | } | 
| 116 | > | // query atom types for their chains of responsibility | 
| 117 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 118 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 119 |  |  | 
| 120 | < | BendType* ForceField::getBendType(const std::string &at1, const std::string &at2, | 
| 120 | > | std::vector<AtomType*>::iterator i; | 
| 121 | > | std::vector<AtomType*>::iterator j; | 
| 122 | > |  | 
| 123 | > | int ii = 0; | 
| 124 | > | int jj = 0; | 
| 125 | > | int bondTypeScore; | 
| 126 | > |  | 
| 127 | > | std::vector<std::pair<int, std::vector<std::string> > > foundBonds; | 
| 128 | > |  | 
| 129 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 130 | > | jj = 0; | 
| 131 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 132 | > |  | 
| 133 | > | bondTypeScore = ii + jj; | 
| 134 | > |  | 
| 135 | > | std::vector<std::string> myKeys; | 
| 136 | > | myKeys.push_back((*i)->getName()); | 
| 137 | > | myKeys.push_back((*j)->getName()); | 
| 138 | > |  | 
| 139 | > | BondType* bondType = bondTypeCont_.find(myKeys); | 
| 140 | > | if (bondType) { | 
| 141 | > | foundBonds.push_back(std::make_pair(bondTypeScore, myKeys)); | 
| 142 | > | } | 
| 143 | > | jj++; | 
| 144 | > | } | 
| 145 | > | ii++; | 
| 146 | > | } | 
| 147 | > |  | 
| 148 | > |  | 
| 149 | > | if (foundBonds.size() > 0) { | 
| 150 | > | // sort the foundBonds by the score: | 
| 151 | > | std::sort(foundBonds.begin(), foundBonds.end()); | 
| 152 | > |  | 
| 153 | > | int bestScore = foundBonds[0].first; | 
| 154 | > | std::vector<std::string> theKeys = foundBonds[0].second; | 
| 155 | > |  | 
| 156 | > | BondType* bestType = bondTypeCont_.find(theKeys); | 
| 157 | > |  | 
| 158 | > | return bestType; | 
| 159 | > | } else { | 
| 160 | > | //if no exact match found, try wild card match | 
| 161 | > | return bondTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 162 | > | } | 
| 163 | > | } | 
| 164 | > | } | 
| 165 | > |  | 
| 166 | > | BendType* ForceField::getBendType(const std::string &at1, | 
| 167 | > | const std::string &at2, | 
| 168 |  | const std::string &at3) { | 
| 169 |  | std::vector<std::string> keys; | 
| 170 |  | keys.push_back(at1); | 
| 176 |  | if (bendType) { | 
| 177 |  | return bendType; | 
| 178 |  | } else { | 
| 179 | < | //if no exact match found, try wild card match | 
| 180 | < | return bendTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 179 | > |  | 
| 180 | > | AtomType* atype1; | 
| 181 | > | AtomType* atype2; | 
| 182 | > | AtomType* atype3; | 
| 183 | > | std::vector<std::string> at1key; | 
| 184 | > | at1key.push_back(at1); | 
| 185 | > | atype1 = atomTypeCont_.find(at1key); | 
| 186 | > |  | 
| 187 | > | std::vector<std::string> at2key; | 
| 188 | > | at2key.push_back(at2); | 
| 189 | > | atype2 = atomTypeCont_.find(at2key); | 
| 190 | > |  | 
| 191 | > | std::vector<std::string> at3key; | 
| 192 | > | at3key.push_back(at3); | 
| 193 | > | atype3 = atomTypeCont_.find(at3key); | 
| 194 | > |  | 
| 195 | > | // query atom types for their chains of responsibility | 
| 196 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 197 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 198 | > | std::vector<AtomType*> at3Chain = atype3->allYourBase(); | 
| 199 | > |  | 
| 200 | > | std::vector<AtomType*>::iterator i; | 
| 201 | > | std::vector<AtomType*>::iterator j; | 
| 202 | > | std::vector<AtomType*>::iterator k; | 
| 203 | > |  | 
| 204 | > | int ii = 0; | 
| 205 | > | int jj = 0; | 
| 206 | > | int kk = 0; | 
| 207 | > | int IKscore; | 
| 208 | > |  | 
| 209 | > | std::vector<tuple3<int, int, std::vector<std::string> > > foundBends; | 
| 210 | > |  | 
| 211 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 212 | > | ii = 0; | 
| 213 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 214 | > | kk = 0; | 
| 215 | > | for (k = at3Chain.begin(); k != at3Chain.end(); k++) { | 
| 216 | > |  | 
| 217 | > | IKscore = ii + kk; | 
| 218 | > |  | 
| 219 | > | std::vector<std::string> myKeys; | 
| 220 | > | myKeys.push_back((*i)->getName()); | 
| 221 | > | myKeys.push_back((*j)->getName()); | 
| 222 | > | myKeys.push_back((*k)->getName()); | 
| 223 | > |  | 
| 224 | > | BendType* bendType = bendTypeCont_.find(myKeys); | 
| 225 | > | if (bendType) { | 
| 226 | > | foundBends.push_back( make_tuple3(jj, IKscore, myKeys) ); | 
| 227 | > | } | 
| 228 | > | kk++; | 
| 229 | > | } | 
| 230 | > | ii++; | 
| 231 | > | } | 
| 232 | > | jj++; | 
| 233 | > | } | 
| 234 | > |  | 
| 235 | > | if (foundBends.size() > 0) { | 
| 236 | > | std::sort(foundBends.begin(), foundBends.end()); | 
| 237 | > | int jscore = foundBends[0].first; | 
| 238 | > | int ikscore = foundBends[0].second; | 
| 239 | > | std::vector<std::string> theKeys = foundBends[0].third; | 
| 240 | > |  | 
| 241 | > | BendType* bestType = bendTypeCont_.find(theKeys); | 
| 242 | > | return bestType; | 
| 243 | > | } else { | 
| 244 | > | //if no exact match found, try wild card match | 
| 245 | > | return bendTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 246 | > | } | 
| 247 |  | } | 
| 248 |  | } | 
| 249 |  |  | 
| 250 | < | TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2, | 
| 251 | < | const std::string &at3, const std::string &at4) { | 
| 250 | > | TorsionType* ForceField::getTorsionType(const std::string &at1, | 
| 251 | > | const std::string &at2, | 
| 252 | > | const std::string &at3, | 
| 253 | > | const std::string &at4) { | 
| 254 |  | std::vector<std::string> keys; | 
| 255 |  | keys.push_back(at1); | 
| 256 |  | keys.push_back(at2); | 
| 257 |  | keys.push_back(at3); | 
| 258 |  | keys.push_back(at4); | 
| 259 |  |  | 
| 260 | + |  | 
| 261 | + | //try exact match first | 
| 262 |  | TorsionType* torsionType = torsionTypeCont_.find(keys); | 
| 263 |  | if (torsionType) { | 
| 264 |  | return torsionType; | 
| 265 |  | } 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_); | 
| 266 |  |  | 
| 267 | + | AtomType* atype1; | 
| 268 | + | AtomType* atype2; | 
| 269 | + | AtomType* atype3; | 
| 270 | + | AtomType* atype4; | 
| 271 | + | std::vector<std::string> at1key; | 
| 272 | + | at1key.push_back(at1); | 
| 273 | + | atype1 = atomTypeCont_.find(at1key); | 
| 274 | + |  | 
| 275 | + | std::vector<std::string> at2key; | 
| 276 | + | at2key.push_back(at2); | 
| 277 | + | atype2 = atomTypeCont_.find(at2key); | 
| 278 | + |  | 
| 279 | + | std::vector<std::string> at3key; | 
| 280 | + | at3key.push_back(at3); | 
| 281 | + | atype3 = atomTypeCont_.find(at3key); | 
| 282 | + |  | 
| 283 | + | std::vector<std::string> at4key; | 
| 284 | + | at4key.push_back(at4); | 
| 285 | + | atype4 = atomTypeCont_.find(at4key); | 
| 286 | + |  | 
| 287 | + | // query atom types for their chains of responsibility | 
| 288 | + | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 289 | + | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 290 | + | std::vector<AtomType*> at3Chain = atype3->allYourBase(); | 
| 291 | + | std::vector<AtomType*> at4Chain = atype4->allYourBase(); | 
| 292 | + |  | 
| 293 | + | std::vector<AtomType*>::iterator i; | 
| 294 | + | std::vector<AtomType*>::iterator j; | 
| 295 | + | std::vector<AtomType*>::iterator k; | 
| 296 | + | std::vector<AtomType*>::iterator l; | 
| 297 | + |  | 
| 298 | + | int ii = 0; | 
| 299 | + | int jj = 0; | 
| 300 | + | int kk = 0; | 
| 301 | + | int ll = 0; | 
| 302 | + | int ILscore; | 
| 303 | + | int JKscore; | 
| 304 | + |  | 
| 305 | + | std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions; | 
| 306 | + |  | 
| 307 | + | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 308 | + | kk = 0; | 
| 309 | + | for (k = at3Chain.begin(); k != at3Chain.end(); k++) { | 
| 310 | + | ii = 0; | 
| 311 | + | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 312 | + | ll = 0; | 
| 313 | + | for (l = at4Chain.begin(); l != at4Chain.end(); l++) { | 
| 314 | + |  | 
| 315 | + | ILscore = ii + ll; | 
| 316 | + | JKscore = jj + kk; | 
| 317 | + |  | 
| 318 | + | std::vector<std::string> myKeys; | 
| 319 | + | myKeys.push_back((*i)->getName()); | 
| 320 | + | myKeys.push_back((*j)->getName()); | 
| 321 | + | myKeys.push_back((*k)->getName()); | 
| 322 | + | myKeys.push_back((*l)->getName()); | 
| 323 | + |  | 
| 324 | + | TorsionType* torsionType = torsionTypeCont_.find(myKeys); | 
| 325 | + | if (torsionType) { | 
| 326 | + | foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) ); | 
| 327 | + | } | 
| 328 | + | ll++; | 
| 329 | + | } | 
| 330 | + | ii++; | 
| 331 | + | } | 
| 332 | + | kk++; | 
| 333 | + | } | 
| 334 | + | jj++; | 
| 335 | + | } | 
| 336 | + |  | 
| 337 | + | if (foundTorsions.size() > 0) { | 
| 338 | + | std::sort(foundTorsions.begin(), foundTorsions.end()); | 
| 339 | + | int jkscore = foundTorsions[0].first; | 
| 340 | + | int ilscore = foundTorsions[0].second; | 
| 341 | + | std::vector<std::string> theKeys = foundTorsions[0].third; | 
| 342 | + |  | 
| 343 | + | TorsionType* bestType = torsionTypeCont_.find(theKeys); | 
| 344 | + | return bestType; | 
| 345 | + | } else { | 
| 346 | + | //if no exact match found, try wild card match | 
| 347 | + | return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 348 | + | } | 
| 349 | + | } | 
| 350 |  | } | 
| 351 |  |  | 
| 352 | < | BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){ | 
| 352 | > | InversionType* ForceField::getInversionType(const std::string &at1, | 
| 353 | > | const std::string &at2, | 
| 354 | > | const std::string &at3, | 
| 355 | > | const std::string &at4) { | 
| 356 | > | std::vector<std::string> keys; | 
| 357 | > | keys.push_back(at1); | 
| 358 | > | keys.push_back(at2); | 
| 359 | > | keys.push_back(at3); | 
| 360 | > | keys.push_back(at4); | 
| 361 | > |  | 
| 362 | > | //try exact match first | 
| 363 | > | InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys); | 
| 364 | > | if (inversionType) { | 
| 365 | > | return inversionType; | 
| 366 | > | } else { | 
| 367 | > |  | 
| 368 | > | AtomType* atype1; | 
| 369 | > | AtomType* atype2; | 
| 370 | > | AtomType* atype3; | 
| 371 | > | AtomType* atype4; | 
| 372 | > | std::vector<std::string> at1key; | 
| 373 | > | at1key.push_back(at1); | 
| 374 | > | atype1 = atomTypeCont_.find(at1key); | 
| 375 | > |  | 
| 376 | > | std::vector<std::string> at2key; | 
| 377 | > | at2key.push_back(at2); | 
| 378 | > | atype2 = atomTypeCont_.find(at2key); | 
| 379 | > |  | 
| 380 | > | std::vector<std::string> at3key; | 
| 381 | > | at3key.push_back(at3); | 
| 382 | > | atype3 = atomTypeCont_.find(at3key); | 
| 383 | > |  | 
| 384 | > | std::vector<std::string> at4key; | 
| 385 | > | at4key.push_back(at4); | 
| 386 | > | atype4 = atomTypeCont_.find(at4key); | 
| 387 | > |  | 
| 388 | > | // query atom types for their chains of responsibility | 
| 389 | > | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 390 | > | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 391 | > | std::vector<AtomType*> at3Chain = atype3->allYourBase(); | 
| 392 | > | std::vector<AtomType*> at4Chain = atype4->allYourBase(); | 
| 393 | > |  | 
| 394 | > | std::vector<AtomType*>::iterator i; | 
| 395 | > | std::vector<AtomType*>::iterator j; | 
| 396 | > | std::vector<AtomType*>::iterator k; | 
| 397 | > | std::vector<AtomType*>::iterator l; | 
| 398 | > |  | 
| 399 | > | int ii = 0; | 
| 400 | > | int jj = 0; | 
| 401 | > | int kk = 0; | 
| 402 | > | int ll = 0; | 
| 403 | > | int Iscore; | 
| 404 | > | int JKLscore; | 
| 405 | > |  | 
| 406 | > | std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions; | 
| 407 | > |  | 
| 408 | > | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 409 | > | kk = 0; | 
| 410 | > | for (k = at3Chain.begin(); k != at3Chain.end(); k++) { | 
| 411 | > | ii = 0; | 
| 412 | > | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 413 | > | ll = 0; | 
| 414 | > | for (l = at4Chain.begin(); l != at4Chain.end(); l++) { | 
| 415 | > |  | 
| 416 | > | Iscore = ii; | 
| 417 | > | JKLscore = jj + kk + ll; | 
| 418 | > |  | 
| 419 | > | std::vector<std::string> myKeys; | 
| 420 | > | myKeys.push_back((*i)->getName()); | 
| 421 | > | myKeys.push_back((*j)->getName()); | 
| 422 | > | myKeys.push_back((*k)->getName()); | 
| 423 | > | myKeys.push_back((*l)->getName()); | 
| 424 | > |  | 
| 425 | > | InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys); | 
| 426 | > | if (inversionType) { | 
| 427 | > | foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) ); | 
| 428 | > | } | 
| 429 | > | ll++; | 
| 430 | > | } | 
| 431 | > | ii++; | 
| 432 | > | } | 
| 433 | > | kk++; | 
| 434 | > | } | 
| 435 | > | jj++; | 
| 436 | > | } | 
| 437 | > |  | 
| 438 | > | if (foundInversions.size() > 0) { | 
| 439 | > | std::sort(foundInversions.begin(), foundInversions.end()); | 
| 440 | > | int iscore = foundInversions[0].first; | 
| 441 | > | int jklscore = foundInversions[0].second; | 
| 442 | > | std::vector<std::string> theKeys = foundInversions[0].third; | 
| 443 | > |  | 
| 444 | > | InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys); | 
| 445 | > | return bestType; | 
| 446 | > | } else { | 
| 447 | > | //if no exact match found, try wild card match | 
| 448 | > | return inversionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 449 | > | } | 
| 450 | > | } | 
| 451 | > | } | 
| 452 | > |  | 
| 453 | > | NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { | 
| 454 | > |  | 
| 455 |  | std::vector<std::string> keys; | 
| 456 |  | keys.push_back(at1); | 
| 457 |  | keys.push_back(at2); | 
| 458 | + |  | 
| 459 | + | //try exact match first | 
| 460 | + | NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys); | 
| 461 | + | if (nbiType) { | 
| 462 | + | return nbiType; | 
| 463 | + | } else { | 
| 464 | + | AtomType* atype1; | 
| 465 | + | AtomType* atype2; | 
| 466 | + | std::vector<std::string> at1key; | 
| 467 | + | at1key.push_back(at1); | 
| 468 | + | atype1 = atomTypeCont_.find(at1key); | 
| 469 | + |  | 
| 470 | + | std::vector<std::string> at2key; | 
| 471 | + | at2key.push_back(at2); | 
| 472 | + | atype2 = atomTypeCont_.find(at2key); | 
| 473 | + |  | 
| 474 | + | // query atom types for their chains of responsibility | 
| 475 | + | std::vector<AtomType*> at1Chain = atype1->allYourBase(); | 
| 476 | + | std::vector<AtomType*> at2Chain = atype2->allYourBase(); | 
| 477 | + |  | 
| 478 | + | std::vector<AtomType*>::iterator i; | 
| 479 | + | std::vector<AtomType*>::iterator j; | 
| 480 | + |  | 
| 481 | + | int ii = 0; | 
| 482 | + | int jj = 0; | 
| 483 | + | int nbiTypeScore; | 
| 484 | + |  | 
| 485 | + | std::vector<std::pair<int, std::vector<std::string> > > foundNBI; | 
| 486 | + |  | 
| 487 | + | for (i = at1Chain.begin(); i != at1Chain.end(); i++) { | 
| 488 | + | jj = 0; | 
| 489 | + | for (j = at2Chain.begin(); j != at2Chain.end(); j++) { | 
| 490 | + |  | 
| 491 | + | nbiTypeScore = ii + jj; | 
| 492 | + |  | 
| 493 | + | std::vector<std::string> myKeys; | 
| 494 | + | myKeys.push_back((*i)->getName()); | 
| 495 | + | myKeys.push_back((*j)->getName()); | 
| 496 | + |  | 
| 497 | + | NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys); | 
| 498 | + | if (nbiType) { | 
| 499 | + | foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys)); | 
| 500 | + | } | 
| 501 | + | jj++; | 
| 502 | + | } | 
| 503 | + | ii++; | 
| 504 | + | } | 
| 505 | + |  | 
| 506 | + |  | 
| 507 | + | if (foundNBI.size() > 0) { | 
| 508 | + | // sort the foundNBI by the score: | 
| 509 | + | std::sort(foundNBI.begin(), foundNBI.end()); | 
| 510 | + |  | 
| 511 | + | int bestScore = foundNBI[0].first; | 
| 512 | + | std::vector<std::string> theKeys = foundNBI[0].second; | 
| 513 | + |  | 
| 514 | + | NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys); | 
| 515 | + | return bestType; | 
| 516 | + | } else { | 
| 517 | + | //if no exact match found, try wild card match | 
| 518 | + | return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 519 | + | } | 
| 520 | + | } | 
| 521 | + | } | 
| 522 | + |  | 
| 523 | + | BondType* ForceField::getExactBondType(const std::string &at1, | 
| 524 | + | const std::string &at2){ | 
| 525 | + | std::vector<std::string> keys; | 
| 526 | + | keys.push_back(at1); | 
| 527 | + | keys.push_back(at2); | 
| 528 |  | return bondTypeCont_.find(keys); | 
| 529 |  | } | 
| 530 | < |  | 
| 531 | < | BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2, | 
| 530 | > |  | 
| 531 | > | BendType* ForceField::getExactBendType(const std::string &at1, | 
| 532 | > | const std::string &at2, | 
| 533 |  | const std::string &at3){ | 
| 534 |  | std::vector<std::string> keys; | 
| 535 |  | keys.push_back(at1); | 
| 537 |  | keys.push_back(at3); | 
| 538 |  | return bendTypeCont_.find(keys); | 
| 539 |  | } | 
| 540 | < |  | 
| 541 | < | TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2, | 
| 542 | < | const std::string &at3, const std::string &at4){ | 
| 540 | > |  | 
| 541 | > | TorsionType* ForceField::getExactTorsionType(const std::string &at1, | 
| 542 | > | const std::string &at2, | 
| 543 | > | const std::string &at3, | 
| 544 | > | const std::string &at4){ | 
| 545 |  | std::vector<std::string> keys; | 
| 546 |  | keys.push_back(at1); | 
| 547 |  | keys.push_back(at2); | 
| 549 |  | keys.push_back(at4); | 
| 550 |  | return torsionTypeCont_.find(keys); | 
| 551 |  | } | 
| 552 | + |  | 
| 553 | + | InversionType* ForceField::getExactInversionType(const std::string &at1, | 
| 554 | + | const std::string &at2, | 
| 555 | + | const std::string &at3, | 
| 556 | + | const std::string &at4){ | 
| 557 | + | std::vector<std::string> keys; | 
| 558 | + | keys.push_back(at1); | 
| 559 | + | keys.push_back(at2); | 
| 560 | + | keys.push_back(at3); | 
| 561 | + | keys.push_back(at4); | 
| 562 | + | return inversionTypeCont_.find(keys); | 
| 563 | + | } | 
| 564 | + |  | 
| 565 | + | NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ | 
| 566 | + | std::vector<std::string> keys; | 
| 567 | + | keys.push_back(at1); | 
| 568 | + | keys.push_back(at2); | 
| 569 | + | return nonBondedInteractionTypeCont_.find(keys); | 
| 570 | + | } | 
| 571 | + |  | 
| 572 | + |  | 
| 573 |  | bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { | 
| 574 |  | std::vector<std::string> keys; | 
| 575 |  | keys.push_back(at); | 
| 576 | + | atypeIdentToName[atomType->getIdent()] = at; | 
| 577 |  | return atomTypeCont_.add(keys, atomType); | 
| 578 |  | } | 
| 579 |  |  | 
| 580 | < | bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) { | 
| 580 | > | bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { | 
| 581 |  | std::vector<std::string> keys; | 
| 582 | + | keys.push_back(at); | 
| 583 | + | atypeIdentToName[atomType->getIdent()] = at; | 
| 584 | + | return atomTypeCont_.replace(keys, atomType); | 
| 585 | + | } | 
| 586 | + |  | 
| 587 | + | bool ForceField::addBondType(const std::string &at1, const std::string &at2, | 
| 588 | + | BondType* bondType) { | 
| 589 | + | std::vector<std::string> keys; | 
| 590 |  | keys.push_back(at1); | 
| 591 |  | keys.push_back(at2); | 
| 592 | < | return bondTypeCont_.add(keys, bondType); | 
| 167 | < |  | 
| 592 | > | return bondTypeCont_.add(keys, bondType); | 
| 593 |  | } | 
| 594 | < |  | 
| 594 | > |  | 
| 595 |  | bool ForceField::addBendType(const std::string &at1, const std::string &at2, | 
| 596 |  | const std::string &at3, BendType* bendType) { | 
| 597 |  | std::vector<std::string> keys; | 
| 600 |  | keys.push_back(at3); | 
| 601 |  | return bendTypeCont_.add(keys, bendType); | 
| 602 |  | } | 
| 603 | < |  | 
| 604 | < | bool ForceField::addTorsionType(const std::string &at1, const std::string &at2, | 
| 605 | < | const std::string &at3, const std::string &at4, TorsionType* torsionType) { | 
| 603 | > |  | 
| 604 | > | bool ForceField::addTorsionType(const std::string &at1, | 
| 605 | > | const std::string &at2, | 
| 606 | > | const std::string &at3, | 
| 607 | > | const std::string &at4, | 
| 608 | > | TorsionType* torsionType) { | 
| 609 |  | std::vector<std::string> keys; | 
| 610 |  | keys.push_back(at1); | 
| 611 |  | keys.push_back(at2); | 
| 614 |  | return torsionTypeCont_.add(keys, torsionType); | 
| 615 |  | } | 
| 616 |  |  | 
| 617 | < | double ForceField::getRcutFromAtomType(AtomType* at) { | 
| 618 | < | /**@todo */ | 
| 619 | < | GenericData* data; | 
| 620 | < | double rcut = 0.0; | 
| 621 | < |  | 
| 622 | < | if (at->isLennardJones()) { | 
| 623 | < | data = at->getPropertyByName("LennardJones"); | 
| 624 | < | if (data != NULL) { | 
| 625 | < | LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); | 
| 626 | < |  | 
| 627 | < | if (ljData != NULL) { | 
| 628 | < | LJParam ljParam = ljData->getData(); | 
| 629 | < |  | 
| 630 | < | //by default use 2.5*sigma as cutoff radius | 
| 631 | < | rcut = 2.5 * ljParam.sigma; | 
| 632 | < |  | 
| 633 | < | } else { | 
| 634 | < | sprintf( painCave.errMsg, | 
| 635 | < | "Can not cast GenericData to LJParam\n"); | 
| 636 | < | painCave.severity = OOPSE_ERROR; | 
| 637 | < | painCave.isFatal = 1; | 
| 638 | < | simError(); | 
| 639 | < | } | 
| 640 | < | } else { | 
| 641 | < | sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); | 
| 642 | < | painCave.severity = OOPSE_ERROR; | 
| 643 | < | painCave.isFatal = 1; | 
| 644 | < | simError(); | 
| 217 | < | } | 
| 617 | > | bool ForceField::addInversionType(const std::string &at1, | 
| 618 | > | const std::string &at2, | 
| 619 | > | const std::string &at3, | 
| 620 | > | const std::string &at4, | 
| 621 | > | InversionType* inversionType) { | 
| 622 | > | std::vector<std::string> keys; | 
| 623 | > | keys.push_back(at1); | 
| 624 | > | keys.push_back(at2); | 
| 625 | > | keys.push_back(at3); | 
| 626 | > | keys.push_back(at4); | 
| 627 | > | return inversionTypeCont_.add(keys, inversionType); | 
| 628 | > | } | 
| 629 | > |  | 
| 630 | > | bool ForceField::addNonBondedInteractionType(const std::string &at1, | 
| 631 | > | const std::string &at2, | 
| 632 | > | NonBondedInteractionType* nbiType) { | 
| 633 | > | std::vector<std::string> keys; | 
| 634 | > | keys.push_back(at1); | 
| 635 | > | keys.push_back(at2); | 
| 636 | > | return nonBondedInteractionTypeCont_.add(keys, nbiType); | 
| 637 | > | } | 
| 638 | > |  | 
| 639 | > | RealType ForceField::getRcutFromAtomType(AtomType* at) { | 
| 640 | > | RealType rcut = 0.0; | 
| 641 | > |  | 
| 642 | > | LennardJonesAdapter lja = LennardJonesAdapter(at); | 
| 643 | > | if (lja.isLennardJones()) { | 
| 644 | > | rcut = 2.5 * lja.getSigma(); | 
| 645 |  | } | 
| 219 | – |  | 
| 646 |  | return rcut; | 
| 647 |  | } | 
| 648 | + |  | 
| 649 |  |  | 
| 223 | – |  | 
| 650 |  | ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { | 
| 651 |  | std::string forceFieldFilename(filename); | 
| 652 |  | ifstrstream* ffStream = new ifstrstream(); | 
| 668 |  | "\tHave you tried setting the FORCE_PARAM_PATH environment " | 
| 669 |  | "variable?\n", | 
| 670 |  | forceFieldFilename.c_str() ); | 
| 671 | < | painCave.severity = OOPSE_ERROR; | 
| 671 | > | painCave.severity = OPENMD_ERROR; | 
| 672 |  | painCave.isFatal = 1; | 
| 673 |  | simError(); | 
| 674 |  | } | 
| 675 |  | } | 
| 250 | – |  | 
| 676 |  | return ffStream; | 
| 252 | – |  | 
| 677 |  | } | 
| 678 |  |  | 
| 679 | < | } //end namespace oopse | 
| 679 | > | } //end namespace OpenMD |