src/UseTheForce/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]  Vardeman & Gezelter, in progress (2009).                        
 */
 
/**
 * @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 "UseTheForce/DarkSide/atype_interface.h"
#include "UseTheForce/DarkSide/fForceOptions_interface.h"
#include "UseTheForce/DarkSide/switcheroo_interface.h"
namespace OpenMD {

  ForceField::ForceField() { 
    char* tempPath; 
    tempPath = getenv("FORCE_PARAM_PATH");

    if (tempPath == NULL) {
      ffPath_ = "ORNULL(FRC_PATH)";
    } else {
      ffPath_ = tempPath;
    }
  }


  ForceField::~ForceField() {
    deleteAtypes();
    deleteSwitch();
  }

  AtomType* ForceField::getAtomType(const std::string &at) {
    std::vector<std::string> keys;
    keys.push_back(at);
    return atomTypeCont_.find(keys);
  }

  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 {
      //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);
    return atomTypeCont_.add(keys, atomType);
  }

  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
    std::vector<std::string> keys;
    keys.push_back(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;
  }

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