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