| 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 | 
  | 
/** | 
| 51 | 
  | 
#include "UseTheForce/ForceField.hpp" | 
| 52 | 
  | 
#include "utils/simError.h" | 
| 53 | 
  | 
#include "utils/Tuple.hpp" | 
| 54 | 
< | 
#include "UseTheForce/DarkSide/atype_interface.h" | 
| 55 | 
< | 
#include "UseTheForce/DarkSide/fForceOptions_interface.h" | 
| 56 | 
< | 
#include "UseTheForce/DarkSide/switcheroo_interface.h" | 
| 57 | 
< | 
namespace oopse { | 
| 54 | 
> | 
namespace OpenMD { | 
| 55 | 
  | 
 | 
| 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 ); | 
| 66 | 
  | 
    } | 
| 67 | 
  | 
  } | 
| 68 | 
  | 
 | 
| 71 | 
– | 
 | 
| 72 | 
– | 
  ForceField::~ForceField() { | 
| 73 | 
– | 
    deleteAtypes(); | 
| 74 | 
– | 
    deleteSwitch(); | 
| 75 | 
– | 
  } | 
| 76 | 
– | 
 | 
| 69 | 
  | 
  AtomType* ForceField::getAtomType(const std::string &at) { | 
| 70 | 
  | 
    std::vector<std::string> keys; | 
| 71 | 
  | 
    keys.push_back(at); | 
| 125 | 
  | 
        ii++; | 
| 126 | 
  | 
      } | 
| 127 | 
  | 
 | 
| 136 | 
– | 
      // sort the foundBonds by the score: | 
| 128 | 
  | 
 | 
| 129 | 
< | 
      std::sort(foundBonds.begin(), foundBonds.end()); | 
| 130 | 
< | 
       | 
| 131 | 
< | 
      int bestScore = foundBonds[0].first; | 
| 132 | 
< | 
      std::vector<std::string> theKeys = foundBonds[0].second; | 
| 133 | 
< | 
 | 
| 134 | 
< | 
      std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";       | 
| 135 | 
< | 
      BondType* bestType = bondTypeCont_.find(theKeys); | 
| 136 | 
< | 
      if (bestType)  | 
| 137 | 
< | 
        return bestType; | 
| 138 | 
< | 
      else { | 
| 139 | 
< | 
        //if no exact match found, try wild card match | 
| 140 | 
< | 
        return bondTypeCont_.find(keys, wildCardAtomTypeName_);       | 
| 129 | 
> | 
      if (foundBonds.size() > 0) { | 
| 130 | 
> | 
        // sort the foundBonds by the score: | 
| 131 | 
> | 
        std::sort(foundBonds.begin(), foundBonds.end()); | 
| 132 | 
> | 
      | 
| 133 | 
> | 
        int bestScore = foundBonds[0].first; | 
| 134 | 
> | 
        std::vector<std::string> theKeys = foundBonds[0].second; | 
| 135 | 
> | 
         | 
| 136 | 
> | 
        BondType* bestType = bondTypeCont_.find(theKeys); | 
| 137 | 
> | 
         | 
| 138 | 
> | 
        return bestType; | 
| 139 | 
> | 
      } else { | 
| 140 | 
> | 
        //if no exact match found, try wild card match | 
| 141 | 
> | 
        return bondTypeCont_.find(keys, wildCardAtomTypeName_);       | 
| 142 | 
  | 
      } | 
| 143 | 
  | 
    } | 
| 144 | 
  | 
  } | 
| 211 | 
  | 
        } | 
| 212 | 
  | 
        jj++; | 
| 213 | 
  | 
      } | 
| 222 | 
– | 
       | 
| 223 | 
– | 
      std::sort(foundBends.begin(), foundBends.end()); | 
| 224 | 
– | 
 | 
| 225 | 
– | 
      int jscore = foundBends[0].first; | 
| 226 | 
– | 
      int ikscore = foundBends[0].second; | 
| 227 | 
– | 
      std::vector<std::string> theKeys = foundBends[0].third; | 
| 228 | 
– | 
 | 
| 229 | 
– | 
      std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";       | 
| 230 | 
– | 
 | 
| 231 | 
– | 
      BendType* bestType = bendTypeCont_.find(theKeys);   | 
| 232 | 
– | 
      if (bestType)  | 
| 233 | 
– | 
        return bestType; | 
| 234 | 
– | 
      else { | 
| 214 | 
  | 
       | 
| 215 | 
+ | 
      if (foundBends.size() > 0) { | 
| 216 | 
+ | 
        std::sort(foundBends.begin(), foundBends.end()); | 
| 217 | 
+ | 
        int jscore = foundBends[0].first; | 
| 218 | 
+ | 
        int ikscore = foundBends[0].second; | 
| 219 | 
+ | 
        std::vector<std::string> theKeys = foundBends[0].third;        | 
| 220 | 
+ | 
         | 
| 221 | 
+ | 
        BendType* bestType = bendTypeCont_.find(theKeys);   | 
| 222 | 
+ | 
        return bestType; | 
| 223 | 
+ | 
      } else {         | 
| 224 | 
  | 
        //if no exact match found, try wild card match | 
| 225 | 
  | 
        return bendTypeCont_.find(keys, wildCardAtomTypeName_);       | 
| 226 | 
  | 
      } | 
| 227 | 
  | 
    } | 
| 228 | 
  | 
  } | 
| 229 | 
  | 
 | 
| 242 | 
– | 
 | 
| 230 | 
  | 
  TorsionType* ForceField::getTorsionType(const std::string &at1,  | 
| 231 | 
  | 
                                          const std::string &at2, | 
| 232 | 
  | 
                                          const std::string &at3,  | 
| 314 | 
  | 
        jj++; | 
| 315 | 
  | 
      } | 
| 316 | 
  | 
       | 
| 317 | 
< | 
      std::sort(foundTorsions.begin(), foundTorsions.end()); | 
| 318 | 
< | 
 | 
| 319 | 
< | 
      int jkscore = foundTorsions[0].first; | 
| 320 | 
< | 
      int ilscore = foundTorsions[0].second; | 
| 321 | 
< | 
      std::vector<std::string> theKeys = foundTorsions[0].third; | 
| 322 | 
< | 
 | 
| 323 | 
< | 
      std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n"; | 
| 324 | 
< | 
 | 
| 338 | 
< | 
       | 
| 339 | 
< | 
      TorsionType* bestType = torsionTypeCont_.find(theKeys); | 
| 340 | 
< | 
      if (bestType) { | 
| 341 | 
< | 
        return bestType; | 
| 317 | 
> | 
      if (foundTorsions.size() > 0) { | 
| 318 | 
> | 
        std::sort(foundTorsions.begin(), foundTorsions.end()); | 
| 319 | 
> | 
        int jkscore = foundTorsions[0].first; | 
| 320 | 
> | 
        int ilscore = foundTorsions[0].second; | 
| 321 | 
> | 
        std::vector<std::string> theKeys = foundTorsions[0].third; | 
| 322 | 
> | 
         | 
| 323 | 
> | 
        TorsionType* bestType = torsionTypeCont_.find(theKeys); | 
| 324 | 
> | 
        return bestType; | 
| 325 | 
  | 
      } else { | 
| 326 | 
  | 
        //if no exact match found, try wild card match | 
| 327 | 
  | 
        return torsionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 340 | 
  | 
    keys.push_back(at4);     | 
| 341 | 
  | 
 | 
| 342 | 
  | 
    //try exact match first | 
| 343 | 
< | 
    InversionType* inversionType = inversionTypeCont_.find(keys); | 
| 343 | 
> | 
    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys); | 
| 344 | 
  | 
    if (inversionType) { | 
| 345 | 
  | 
      return inversionType; | 
| 346 | 
  | 
    } else { | 
| 402 | 
  | 
              myKeys.push_back((*k)->getName()); | 
| 403 | 
  | 
              myKeys.push_back((*l)->getName()); | 
| 404 | 
  | 
               | 
| 405 | 
< | 
              InversionType* inversionType = inversionTypeCont_.find(myKeys); | 
| 405 | 
> | 
              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys); | 
| 406 | 
  | 
              if (inversionType) {  | 
| 407 | 
  | 
                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) ); | 
| 408 | 
  | 
              } | 
| 414 | 
  | 
        } | 
| 415 | 
  | 
        jj++; | 
| 416 | 
  | 
      } | 
| 417 | 
< | 
       | 
| 418 | 
< | 
      std::sort(foundInversions.begin(), foundInversions.end()); | 
| 419 | 
< | 
       | 
| 420 | 
< | 
      int iscore = foundInversions[0].first; | 
| 421 | 
< | 
      int jklscore = foundInversions[0].second; | 
| 422 | 
< | 
      std::vector<std::string> theKeys = foundInversions[0].third; | 
| 423 | 
< | 
       | 
| 424 | 
< | 
      std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n"; | 
| 425 | 
< | 
       | 
| 443 | 
< | 
       | 
| 444 | 
< | 
      InversionType* bestType = inversionTypeCont_.find(theKeys); | 
| 445 | 
< | 
      if (bestType) { | 
| 446 | 
< | 
        return bestType; | 
| 417 | 
> | 
          | 
| 418 | 
> | 
      if (foundInversions.size() > 0) { | 
| 419 | 
> | 
        std::sort(foundInversions.begin(), foundInversions.end()); | 
| 420 | 
> | 
        int iscore = foundInversions[0].first; | 
| 421 | 
> | 
        int jklscore = foundInversions[0].second; | 
| 422 | 
> | 
        std::vector<std::string> theKeys = foundInversions[0].third; | 
| 423 | 
> | 
         | 
| 424 | 
> | 
        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys); | 
| 425 | 
> | 
        return bestType; | 
| 426 | 
  | 
      } else { | 
| 427 | 
  | 
        //if no exact match found, try wild card match | 
| 428 | 
  | 
        return inversionTypeCont_.find(keys, wildCardAtomTypeName_); | 
| 501 | 
  | 
    return atomTypeCont_.add(keys, atomType); | 
| 502 | 
  | 
  } | 
| 503 | 
  | 
 | 
| 504 | 
+ | 
  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { | 
| 505 | 
+ | 
    std::vector<std::string> keys; | 
| 506 | 
+ | 
    keys.push_back(at); | 
| 507 | 
+ | 
    return atomTypeCont_.replace(keys, atomType); | 
| 508 | 
+ | 
  } | 
| 509 | 
+ | 
 | 
| 510 | 
  | 
  bool ForceField::addBondType(const std::string &at1, const std::string &at2, | 
| 511 | 
  | 
                               BondType* bondType) { | 
| 512 | 
  | 
    std::vector<std::string> keys; | 
| 578 | 
  | 
        } else { | 
| 579 | 
  | 
          sprintf( painCave.errMsg, | 
| 580 | 
  | 
                   "Can not cast GenericData to LJParam\n"); | 
| 581 | 
< | 
          painCave.severity = OOPSE_ERROR; | 
| 581 | 
> | 
          painCave.severity = OPENMD_ERROR; | 
| 582 | 
  | 
          painCave.isFatal = 1; | 
| 583 | 
  | 
          simError();           | 
| 584 | 
  | 
        }             | 
| 585 | 
  | 
      } else { | 
| 586 | 
  | 
        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); | 
| 587 | 
< | 
        painCave.severity = OOPSE_ERROR; | 
| 587 | 
> | 
        painCave.severity = OPENMD_ERROR; | 
| 588 | 
  | 
        painCave.isFatal = 1; | 
| 589 | 
  | 
        simError();           | 
| 590 | 
  | 
      } | 
| 614 | 
  | 
                 "\tHave you tried setting the FORCE_PARAM_PATH environment " | 
| 615 | 
  | 
                 "variable?\n", | 
| 616 | 
  | 
                 forceFieldFilename.c_str() ); | 
| 617 | 
< | 
        painCave.severity = OOPSE_ERROR; | 
| 617 | 
> | 
        painCave.severity = OPENMD_ERROR; | 
| 618 | 
  | 
        painCave.isFatal = 1; | 
| 619 | 
  | 
        simError(); | 
| 620 | 
  | 
      } | 
| 622 | 
  | 
    return ffStream; | 
| 623 | 
  | 
  } | 
| 624 | 
  | 
 | 
| 625 | 
< | 
  void ForceField::setFortranForceOptions(){ | 
| 641 | 
< | 
    ForceOptions theseFortranOptions; | 
| 642 | 
< | 
    forceFieldOptions_.makeFortranOptions(theseFortranOptions); | 
| 643 | 
< | 
    setfForceOptions(&theseFortranOptions); | 
| 644 | 
< | 
  } | 
| 645 | 
< | 
} //end namespace oopse | 
| 625 | 
> | 
} //end namespace OpenMD |