src/brains/ForceField.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */
 
/**
 * @file ForceField.cpp
 * @author tlin
 * @date 11/04/2004
 * @time 22:51am
 * @version 1.0
 */
  
#include <algorithm>
#include "UseTheForce/ForceField.hpp"
#include "utils/simError.h"
#include "utils/Tuple.hpp"
#include "types/LennardJonesAdapter.hpp"

namespace OpenMD {

  ForceField::ForceField() { 

    char* tempPath; 
    tempPath = getenv("FORCE_PARAM_PATH");
    
    if (tempPath == NULL) {
      //convert a macro from compiler to a string in c++
      STR_DEFINE(ffPath_, FRC_PATH );
    } else {
      ffPath_ = tempPath;
    }
  }

  /**
   * getAtomType by string
   *
   * finds the requested atom type in this force field using the string
   * name of the atom type.
   */
  AtomType* ForceField::getAtomType(const std::string &at) {
    std::vector<std::string> keys;
    keys.push_back(at);
    return atomTypeCont_.find(keys);
  }

  /**
   * getAtomType by ident
   *
   * finds the requested atom type in this force field using the
   * integer ident instead of the string name of the atom type.
   */
  AtomType* ForceField::getAtomType(int ident) {   
    std::string at = atypeIdentToName.find(ident)->second;
    return getAtomType(at);
  }

  BondType* ForceField::getBondType(const std::string &at1, 
                                    const std::string &at2) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    

    //try exact match first
    BondType* bondType = bondTypeCont_.find(keys);
    if (bondType) {
      return bondType;
    } else {
      AtomType* atype1;
      AtomType* atype2;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
  
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;

      int ii = 0;
      int jj = 0;
      int bondTypeScore;

      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;

      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
        jj = 0;
        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {

          bondTypeScore = ii + jj;

          std::vector<std::string> myKeys;
          myKeys.push_back((*i)->getName());
          myKeys.push_back((*j)->getName());

          BondType* bondType = bondTypeCont_.find(myKeys);
          if (bondType) {
            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
          }
          jj++;
        }
        ii++;
      }


      if (foundBonds.size() > 0) {
        // sort the foundBonds by the score:
        std::sort(foundBonds.begin(), foundBonds.end());
     
        int bestScore = foundBonds[0].first;
        std::vector<std::string> theKeys = foundBonds[0].second;
        
        BondType* bestType = bondTypeCont_.find(theKeys);
        
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
      }
    }
  }
  
  BendType* ForceField::getBendType(const std::string &at1, 
                                    const std::string &at2,
                                    const std::string &at3) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    

    //try exact match first
    BendType* bendType = bendTypeCont_.find(keys);
    if (bendType) {
      return bendType;
    } else {

      AtomType* atype1;
      AtomType* atype2;
      AtomType* atype3;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
  
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);

      std::vector<std::string> at3key;
      at3key.push_back(at3);
      atype3 = atomTypeCont_.find(at3key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      std::vector<AtomType*> at3Chain = atype3->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      std::vector<AtomType*>::iterator k;

      int ii = 0;
      int jj = 0;
      int kk = 0;
      int IKscore;

      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;

      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
        ii = 0;
        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
          kk = 0;
          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
          
            IKscore = ii + kk;

            std::vector<std::string> myKeys;
            myKeys.push_back((*i)->getName());
            myKeys.push_back((*j)->getName());
            myKeys.push_back((*k)->getName());

            BendType* bendType = bendTypeCont_.find(myKeys);
            if (bendType) { 
              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
            }
            kk++;
          }
          ii++;
        }
        jj++;
      }
      
      if (foundBends.size() > 0) {
        std::sort(foundBends.begin(), foundBends.end());
        int jscore = foundBends[0].first;
        int ikscore = foundBends[0].second;
        std::vector<std::string> theKeys = foundBends[0].third;       
        
        BendType* bestType = bendTypeCont_.find(theKeys);  
        return bestType;
      } else {        
        //if no exact match found, try wild card match
        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
      }
    }
  }

  TorsionType* ForceField::getTorsionType(const std::string &at1, 
                                          const std::string &at2,
                                          const std::string &at3, 
                                          const std::string &at4) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    


    //try exact match first
    TorsionType* torsionType = torsionTypeCont_.find(keys);
    if (torsionType) {
      return torsionType;
    } else {

      AtomType* atype1;
      AtomType* atype2;
      AtomType* atype3;
      AtomType* atype4;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
  
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);

      std::vector<std::string> at3key;
      at3key.push_back(at3);
      atype3 = atomTypeCont_.find(at3key);

      std::vector<std::string> at4key;
      at4key.push_back(at4);
      atype4 = atomTypeCont_.find(at4key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      std::vector<AtomType*> at3Chain = atype3->allYourBase();
      std::vector<AtomType*> at4Chain = atype4->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      std::vector<AtomType*>::iterator k;
      std::vector<AtomType*>::iterator l;

      int ii = 0;
      int jj = 0;
      int kk = 0;
      int ll = 0;
      int ILscore;
      int JKscore;

      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;

      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
        kk = 0;
        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
          ii = 0;       
          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
            ll = 0;
            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
          
              ILscore = ii + ll;
              JKscore = jj + kk;

              std::vector<std::string> myKeys;
              myKeys.push_back((*i)->getName());
              myKeys.push_back((*j)->getName());
              myKeys.push_back((*k)->getName());
              myKeys.push_back((*l)->getName());

              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
              if (torsionType) { 
                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
              }
              ll++;
            }
            ii++;
          }
          kk++;
        }
        jj++;
      }
      
      if (foundTorsions.size() > 0) {
        std::sort(foundTorsions.begin(), foundTorsions.end());
        int jkscore = foundTorsions[0].first;
        int ilscore = foundTorsions[0].second;
        std::vector<std::string> theKeys = foundTorsions[0].third;
        
        TorsionType* bestType = torsionTypeCont_.find(theKeys);
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
      }
    }
  }

  InversionType* ForceField::getInversionType(const std::string &at1, 
                                              const std::string &at2,
                                              const std::string &at3, 
                                              const std::string &at4) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    

    //try exact match first
    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
    if (inversionType) {
      return inversionType;
    } else {
      
      AtomType* atype1;
      AtomType* atype2;
      AtomType* atype3;
      AtomType* atype4;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
      
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);
      
      std::vector<std::string> at3key;
      at3key.push_back(at3);
      atype3 = atomTypeCont_.find(at3key);
      
      std::vector<std::string> at4key;
      at4key.push_back(at4);
      atype4 = atomTypeCont_.find(at4key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      std::vector<AtomType*> at3Chain = atype3->allYourBase();
      std::vector<AtomType*> at4Chain = atype4->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      std::vector<AtomType*>::iterator k;
      std::vector<AtomType*>::iterator l;

      int ii = 0;
      int jj = 0;
      int kk = 0;
      int ll = 0;
      int Iscore;
      int JKLscore;
      
      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
      
      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
        kk = 0;
        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
          ii = 0;       
          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
            ll = 0;
            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
              
              Iscore = ii;
              JKLscore = jj + kk + ll;
              
              std::vector<std::string> myKeys;
              myKeys.push_back((*i)->getName());
              myKeys.push_back((*j)->getName());
              myKeys.push_back((*k)->getName());
              myKeys.push_back((*l)->getName());
              
              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
              if (inversionType) { 
                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
              }
              ll++;
            }
            ii++;
          }
          kk++;
        }
        jj++;
      }
         
      if (foundInversions.size() > 0) {
        std::sort(foundInversions.begin(), foundInversions.end());
        int iscore = foundInversions[0].first;
        int jklscore = foundInversions[0].second;
        std::vector<std::string> theKeys = foundInversions[0].third;
        
        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
      }
    }
  }
  
  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
    
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    
    //try exact match first
    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
    if (nbiType) {
      return nbiType;
    } else {
      AtomType* atype1;
      AtomType* atype2;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
      
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);
      
      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      
      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      
      int ii = 0;
      int jj = 0;
      int nbiTypeScore;
      
      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
      
      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
        jj = 0;
        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
          
          nbiTypeScore = ii + jj;
          
          std::vector<std::string> myKeys;
          myKeys.push_back((*i)->getName());
          myKeys.push_back((*j)->getName());
          
          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
          if (nbiType) {
            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
          }
          jj++;
        }
        ii++;
      }
      
      
      if (foundNBI.size() > 0) {
        // sort the foundNBI by the score:
        std::sort(foundNBI.begin(), foundNBI.end());
        
        int bestScore = foundNBI[0].first;
        std::vector<std::string> theKeys = foundNBI[0].second;
        
        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
      }
    }
  }
  
  BondType* ForceField::getExactBondType(const std::string &at1, 
                                         const std::string &at2){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return bondTypeCont_.find(keys);
  }
  
  BendType* ForceField::getExactBendType(const std::string &at1, 
                                         const std::string &at2,
                                         const std::string &at3){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    return bendTypeCont_.find(keys);
  }
  
  TorsionType* ForceField::getExactTorsionType(const std::string &at1, 
                                               const std::string &at2,
                                               const std::string &at3, 
                                               const std::string &at4){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);   
    return torsionTypeCont_.find(keys);
  }
  
  InversionType* ForceField::getExactInversionType(const std::string &at1, 
                                                   const std::string &at2,
                                                   const std::string &at3, 
                                                   const std::string &at4){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);   
    return inversionTypeCont_.find(keys);
  }
  
  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return nonBondedInteractionTypeCont_.find(keys);
  }
  

  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
    std::vector<std::string> keys;
    keys.push_back(at);
    atypeIdentToName[atomType->getIdent()] = at;
    return atomTypeCont_.add(keys, atomType);
  }

  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
    std::vector<std::string> keys;
    keys.push_back(at);
    atypeIdentToName[atomType->getIdent()] = at;
    return atomTypeCont_.replace(keys, atomType);
  }

  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
                               BondType* bondType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return bondTypeCont_.add(keys, bondType);    
  }
  
  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
                               const std::string &at3, BendType* bendType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    return bendTypeCont_.add(keys, bendType);
  }
  
  bool ForceField::addTorsionType(const std::string &at1, 
                                  const std::string &at2,
                                  const std::string &at3, 
                                  const std::string &at4, 
                                  TorsionType* torsionType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    
    return torsionTypeCont_.add(keys, torsionType);
  }

  bool ForceField::addInversionType(const std::string &at1, 
                                    const std::string &at2,
                                    const std::string &at3, 
                                    const std::string &at4, 
                                    InversionType* inversionType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    
    return inversionTypeCont_.add(keys, inversionType);
  }
  
  bool ForceField::addNonBondedInteractionType(const std::string &at1, 
                                               const std::string &at2, 
                                               NonBondedInteractionType* nbiType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return nonBondedInteractionTypeCont_.add(keys, nbiType);
  }
  
  RealType ForceField::getRcutFromAtomType(AtomType* at) {
    RealType rcut = 0.0;
    
    LennardJonesAdapter lja = LennardJonesAdapter(at);
    if (lja.isLennardJones()) {
      rcut = 2.5 * lja.getSigma();
    }
    return rcut;    
  }
  

  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
    std::string forceFieldFilename(filename);
    ifstrstream* ffStream = new ifstrstream();
    
    //try to open the force filed file in current directory first    
    ffStream->open(forceFieldFilename.c_str());
    if(!ffStream->is_open()){

      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
      ffStream->open( forceFieldFilename.c_str() );

      //if current directory does not contain the force field file,
      //try to open it in the path        
      if(!ffStream->is_open()){

        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file:\n"
                 "\t%s\n"
                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
                 "variable?\n",
                 forceFieldFilename.c_str() );
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();
      }
    }  
    return ffStream;
  }

} //end namespace OpenMD
Revision:	1710
Committed:	Fri May 18 21:44:02 2012 UTC (13 years, 5 months ago) by gezelter
Original Path:	*branches/development/src/UseTheForce/ForceField.cpp*
File size:	20941 byte(s)
Log Message:	Added an adapter layer between the AtomType and the rest of the code to handle the bolt-on capabilities of new types. Fixed a long-standing bug in how storageLayout was being set to the maximum possible value. Started to add infrastructure for Polarizable and fluc-Q calculations.
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
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.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
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	/**
44	* @file ForceField.cpp
45	* @author tlin
46	* @date 11/04/2004
47	* @time 22:51am
48	* @version 1.0
49	*/
50
51	#include <algorithm>
52	#include "UseTheForce/ForceField.hpp"
53	#include "utils/simError.h"
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
64	if (tempPath == NULL) {
65	//convert a macro from compiler to a string in c++
66	STR_DEFINE(ffPath_, FRC_PATH );
67	} else {
68	ffPath_ = tempPath;
69	}
70	}
71
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	/**
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);
100
101	//try exact match first
102	BondType* bondType = bondTypeCont_.find(keys);
103	if (bondType) {
104	return bondType;
105	} else {
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	// query atom types for their chains of responsibility
117	std::vector<AtomType*> at1Chain = atype1->allYourBase();
118	std::vector<AtomType*> at2Chain = atype2->allYourBase();
119
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);
171	keys.push_back(at2);
172	keys.push_back(at3);
173
174	//try exact match first
175	BendType* bendType = bendTypeCont_.find(keys);
176	if (bendType) {
177	return bendType;
178	} else {
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,
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 {
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	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,
532	const std::string &at2,
533	const std::string &at3){
534	std::vector<std::string> keys;
535	keys.push_back(at1);
536	keys.push_back(at2);
537	keys.push_back(at3);
538	return bendTypeCont_.find(keys);
539	}
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);
548	keys.push_back(at3);
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::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);
593	}
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;
598	keys.push_back(at1);
599	keys.push_back(at2);
600	keys.push_back(at3);
601	return bendTypeCont_.add(keys, bendType);
602	}
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);
612	keys.push_back(at3);
613	keys.push_back(at4);
614	return torsionTypeCont_.add(keys, torsionType);
615	}
616
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	}
646	return rcut;
647	}
648
649
650	ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
651	std::string forceFieldFilename(filename);
652	ifstrstream* ffStream = new ifstrstream();
653
654	//try to open the force filed file in current directory first
655	ffStream->open(forceFieldFilename.c_str());
656	if(!ffStream->is_open()){
657
658	forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
659	ffStream->open( forceFieldFilename.c_str() );
660
661	//if current directory does not contain the force field file,
662	//try to open it in the path
663	if(!ffStream->is_open()){
664
665	sprintf( painCave.errMsg,
666	"Error opening the force field parameter file:\n"
667	"\t%s\n"
668	"\tHave you tried setting the FORCE_PARAM_PATH environment "
669	"variable?\n",
670	forceFieldFilename.c_str() );
671	painCave.severity = OPENMD_ERROR;
672	painCave.isFatal = 1;
673	simError();
674	}
675	}
676	return ffStream;
677	}
678
679	} //end namespace OpenMD