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