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