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. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. 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.
 */
 
/**
 * @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 oopse {

  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());

          std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
          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;
        
        std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
        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;
        
        std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
        
        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;
        
        std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
                
        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_.find(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_.find(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;
        
        std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
                
        InversionType* bestType = inversionTypeCont_.find(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 = OOPSE_ERROR;
          painCave.isFatal = 1;
          simError();          
        }            
      } else {
        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
        painCave.severity = OOPSE_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 = OOPSE_ERROR;
        painCave.isFatal = 1;
        simError();
      }
    }  
    return ffStream;
  }

  void ForceField::setFortranForceOptions(){
    ForceOptions theseFortranOptions;
    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
    setfForceOptions(&theseFortranOptions);
  }
} //end namespace oopse
Revision:	1282
Committed:	Wed Jul 30 18:11:19 2008 UTC (16 years, 9 months ago) by gezelter
File size:	20258 byte(s)
Log Message:	Many fixes
#	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. Acknowledgement of the program authors must be made in any
10	* publication of scientific results based in part on use of the
11	* program. An acceptable form of acknowledgement is citation of
12	* the article in which the program was described (Matthew
13	* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14	* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15	* Parallel Simulation Engine for Molecular Dynamics,"
16	* J. Comput. Chem. 26, pp. 252-271 (2005))
17	*
18	* 2. Redistributions of source code must retain the above copyright
19	* notice, this list of conditions and the following disclaimer.
20	*
21	* 3. Redistributions in binary form must reproduce the above copyright
22	* notice, this list of conditions and the following disclaimer in the
23	* documentation and/or other materials provided with the
24	* distribution.
25	*
26	* This software is provided "AS IS," without a warranty of any
27	* kind. All express or implied conditions, representations and
28	* warranties, including any implied warranty of merchantability,
29	* fitness for a particular purpose or non-infringement, are hereby
30	* excluded. The University of Notre Dame and its licensors shall not
31	* be liable for any damages suffered by licensee as a result of
32	* using, modifying or distributing the software or its
33	* derivatives. In no event will the University of Notre Dame or its
34	* licensors be liable for any lost revenue, profit or data, or for
35	* direct, indirect, special, consequential, incidental or punitive
36	* damages, however caused and regardless of the theory of liability,
37	* arising out of the use of or inability to use software, even if the
38	* University of Notre Dame has been advised of the possibility of
39	* such damages.
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 oopse {
58
59	ForceField::ForceField() {
60	char* tempPath;
61	tempPath = getenv("FORCE_PARAM_PATH");
62
63	if (tempPath == NULL) {
64	//convert a macro from compiler to a string in c++
65	STR_DEFINE(ffPath_, FRC_PATH );
66	} else {
67	ffPath_ = tempPath;
68	}
69	}
70
71
72	ForceField::~ForceField() {
73	deleteAtypes();
74	deleteSwitch();
75	}
76
77	AtomType* ForceField::getAtomType(const std::string &at) {
78	std::vector<std::string> keys;
79	keys.push_back(at);
80	return atomTypeCont_.find(keys);
81	}
82
83	BondType* ForceField::getBondType(const std::string &at1,
84	const std::string &at2) {
85	std::vector<std::string> keys;
86	keys.push_back(at1);
87	keys.push_back(at2);
88
89	//try exact match first
90	BondType* bondType = bondTypeCont_.find(keys);
91	if (bondType) {
92	return bondType;
93	} else {
94	AtomType* atype1;
95	AtomType* atype2;
96	std::vector<std::string> at1key;
97	at1key.push_back(at1);
98	atype1 = atomTypeCont_.find(at1key);
99
100	std::vector<std::string> at2key;
101	at2key.push_back(at2);
102	atype2 = atomTypeCont_.find(at2key);
103
104	// query atom types for their chains of responsibility
105	std::vector<AtomType*> at1Chain = atype1->allYourBase();
106	std::vector<AtomType*> at2Chain = atype2->allYourBase();
107
108	std::vector<AtomType*>::iterator i;
109	std::vector<AtomType*>::iterator j;
110
111	int ii = 0;
112	int jj = 0;
113	int bondTypeScore;
114
115	std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
116
117	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
118	jj = 0;
119	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
120
121	bondTypeScore = ii + jj;
122
123	std::vector<std::string> myKeys;
124	myKeys.push_back((*i)->getName());
125	myKeys.push_back((*j)->getName());
126
127	std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
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	std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1] << "\t(score = "<< bestScore << ")\n";
146	BondType* bestType = bondTypeCont_.find(theKeys);
147
148	return bestType;
149	} else {
150	//if no exact match found, try wild card match
151	return bondTypeCont_.find(keys, wildCardAtomTypeName_);
152	}
153	}
154	}
155
156	BendType* ForceField::getBendType(const std::string &at1,
157	const std::string &at2,
158	const std::string &at3) {
159	std::vector<std::string> keys;
160	keys.push_back(at1);
161	keys.push_back(at2);
162	keys.push_back(at3);
163
164	//try exact match first
165	BendType* bendType = bendTypeCont_.find(keys);
166	if (bendType) {
167	return bendType;
168	} else {
169
170	AtomType* atype1;
171	AtomType* atype2;
172	AtomType* atype3;
173	std::vector<std::string> at1key;
174	at1key.push_back(at1);
175	atype1 = atomTypeCont_.find(at1key);
176
177	std::vector<std::string> at2key;
178	at2key.push_back(at2);
179	atype2 = atomTypeCont_.find(at2key);
180
181	std::vector<std::string> at3key;
182	at3key.push_back(at3);
183	atype3 = atomTypeCont_.find(at3key);
184
185	// query atom types for their chains of responsibility
186	std::vector<AtomType*> at1Chain = atype1->allYourBase();
187	std::vector<AtomType*> at2Chain = atype2->allYourBase();
188	std::vector<AtomType*> at3Chain = atype3->allYourBase();
189
190	std::vector<AtomType*>::iterator i;
191	std::vector<AtomType*>::iterator j;
192	std::vector<AtomType*>::iterator k;
193
194	int ii = 0;
195	int jj = 0;
196	int kk = 0;
197	int IKscore;
198
199	std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
200
201	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
202	ii = 0;
203	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
204	kk = 0;
205	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
206
207	IKscore = ii + kk;
208
209	std::vector<std::string> myKeys;
210	myKeys.push_back((*i)->getName());
211	myKeys.push_back((*j)->getName());
212	myKeys.push_back((*k)->getName());
213
214	BendType* bendType = bendTypeCont_.find(myKeys);
215	if (bendType) {
216	foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
217	}
218	kk++;
219	}
220	ii++;
221	}
222	jj++;
223	}
224
225	if (foundBends.size() > 0) {
226	std::sort(foundBends.begin(), foundBends.end());
227	int jscore = foundBends[0].first;
228	int ikscore = foundBends[0].second;
229	std::vector<std::string> theKeys = foundBends[0].third;
230
231	std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";
232
233	BendType* bestType = bendTypeCont_.find(theKeys);
234	return bestType;
235	} else {
236	//if no exact match found, try wild card match
237	return bendTypeCont_.find(keys, wildCardAtomTypeName_);
238	}
239	}
240	}
241
242	TorsionType* ForceField::getTorsionType(const std::string &at1,
243	const std::string &at2,
244	const std::string &at3,
245	const std::string &at4) {
246	std::vector<std::string> keys;
247	keys.push_back(at1);
248	keys.push_back(at2);
249	keys.push_back(at3);
250	keys.push_back(at4);
251
252
253	//try exact match first
254	TorsionType* torsionType = torsionTypeCont_.find(keys);
255	if (torsionType) {
256	return torsionType;
257	} else {
258
259	AtomType* atype1;
260	AtomType* atype2;
261	AtomType* atype3;
262	AtomType* atype4;
263	std::vector<std::string> at1key;
264	at1key.push_back(at1);
265	atype1 = atomTypeCont_.find(at1key);
266
267	std::vector<std::string> at2key;
268	at2key.push_back(at2);
269	atype2 = atomTypeCont_.find(at2key);
270
271	std::vector<std::string> at3key;
272	at3key.push_back(at3);
273	atype3 = atomTypeCont_.find(at3key);
274
275	std::vector<std::string> at4key;
276	at4key.push_back(at4);
277	atype4 = atomTypeCont_.find(at4key);
278
279	// query atom types for their chains of responsibility
280	std::vector<AtomType*> at1Chain = atype1->allYourBase();
281	std::vector<AtomType*> at2Chain = atype2->allYourBase();
282	std::vector<AtomType*> at3Chain = atype3->allYourBase();
283	std::vector<AtomType*> at4Chain = atype4->allYourBase();
284
285	std::vector<AtomType*>::iterator i;
286	std::vector<AtomType*>::iterator j;
287	std::vector<AtomType*>::iterator k;
288	std::vector<AtomType*>::iterator l;
289
290	int ii = 0;
291	int jj = 0;
292	int kk = 0;
293	int ll = 0;
294	int ILscore;
295	int JKscore;
296
297	std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
298
299	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
300	kk = 0;
301	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
302	ii = 0;
303	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
304	ll = 0;
305	for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
306
307	ILscore = ii + ll;
308	JKscore = jj + kk;
309
310	std::vector<std::string> myKeys;
311	myKeys.push_back((*i)->getName());
312	myKeys.push_back((*j)->getName());
313	myKeys.push_back((*k)->getName());
314	myKeys.push_back((*l)->getName());
315
316	TorsionType* torsionType = torsionTypeCont_.find(myKeys);
317	if (torsionType) {
318	foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
319	}
320	ll++;
321	}
322	ii++;
323	}
324	kk++;
325	}
326	jj++;
327	}
328
329	if (foundTorsions.size() > 0) {
330	std::sort(foundTorsions.begin(), foundTorsions.end());
331	int jkscore = foundTorsions[0].first;
332	int ilscore = foundTorsions[0].second;
333	std::vector<std::string> theKeys = foundTorsions[0].third;
334
335	std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
336
337	TorsionType* bestType = torsionTypeCont_.find(theKeys);
338	return bestType;
339	} else {
340	//if no exact match found, try wild card match
341	return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
342	}
343	}
344	}
345
346	InversionType* ForceField::getInversionType(const std::string &at1,
347	const std::string &at2,
348	const std::string &at3,
349	const std::string &at4) {
350	std::vector<std::string> keys;
351	keys.push_back(at1);
352	keys.push_back(at2);
353	keys.push_back(at3);
354	keys.push_back(at4);
355
356	//try exact match first
357	InversionType* inversionType = inversionTypeCont_.find(keys);
358	if (inversionType) {
359	return inversionType;
360	} else {
361
362	AtomType* atype1;
363	AtomType* atype2;
364	AtomType* atype3;
365	AtomType* atype4;
366	std::vector<std::string> at1key;
367	at1key.push_back(at1);
368	atype1 = atomTypeCont_.find(at1key);
369
370	std::vector<std::string> at2key;
371	at2key.push_back(at2);
372	atype2 = atomTypeCont_.find(at2key);
373
374	std::vector<std::string> at3key;
375	at3key.push_back(at3);
376	atype3 = atomTypeCont_.find(at3key);
377
378	std::vector<std::string> at4key;
379	at4key.push_back(at4);
380	atype4 = atomTypeCont_.find(at4key);
381
382	// query atom types for their chains of responsibility
383	std::vector<AtomType*> at1Chain = atype1->allYourBase();
384	std::vector<AtomType*> at2Chain = atype2->allYourBase();
385	std::vector<AtomType*> at3Chain = atype3->allYourBase();
386	std::vector<AtomType*> at4Chain = atype4->allYourBase();
387
388	std::vector<AtomType*>::iterator i;
389	std::vector<AtomType*>::iterator j;
390	std::vector<AtomType*>::iterator k;
391	std::vector<AtomType*>::iterator l;
392
393	int ii = 0;
394	int jj = 0;
395	int kk = 0;
396	int ll = 0;
397	int Iscore;
398	int JKLscore;
399
400	std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
401
402	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
403	kk = 0;
404	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
405	ii = 0;
406	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
407	ll = 0;
408	for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
409
410	Iscore = ii;
411	JKLscore = jj + kk + ll;
412
413	std::vector<std::string> myKeys;
414	myKeys.push_back((*i)->getName());
415	myKeys.push_back((*j)->getName());
416	myKeys.push_back((*k)->getName());
417	myKeys.push_back((*l)->getName());
418
419	InversionType* inversionType = inversionTypeCont_.find(myKeys);
420	if (inversionType) {
421	foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
422	}
423	ll++;
424	}
425	ii++;
426	}
427	kk++;
428	}
429	jj++;
430	}
431
432	if (foundInversions.size() > 0) {
433	std::sort(foundInversions.begin(), foundInversions.end());
434	int iscore = foundInversions[0].first;
435	int jklscore = foundInversions[0].second;
436	std::vector<std::string> theKeys = foundInversions[0].third;
437
438	std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
439
440	InversionType* bestType = inversionTypeCont_.find(theKeys);
441	return bestType;
442	} else {
443	//if no exact match found, try wild card match
444	return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
445	}
446	}
447	}
448
449	NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
450	std::vector<std::string> keys;
451	keys.push_back(at1);
452	keys.push_back(at2);
453
454	//try exact match first
455	NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
456	if (nbiType) {
457	return nbiType;
458	} else {
459	//if no exact match found, try wild card match
460	return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
461	}
462	}
463
464	BondType* ForceField::getExactBondType(const std::string &at1,
465	const std::string &at2){
466	std::vector<std::string> keys;
467	keys.push_back(at1);
468	keys.push_back(at2);
469	return bondTypeCont_.find(keys);
470	}
471
472	BendType* ForceField::getExactBendType(const std::string &at1,
473	const std::string &at2,
474	const std::string &at3){
475	std::vector<std::string> keys;
476	keys.push_back(at1);
477	keys.push_back(at2);
478	keys.push_back(at3);
479	return bendTypeCont_.find(keys);
480	}
481
482	TorsionType* ForceField::getExactTorsionType(const std::string &at1,
483	const std::string &at2,
484	const std::string &at3,
485	const std::string &at4){
486	std::vector<std::string> keys;
487	keys.push_back(at1);
488	keys.push_back(at2);
489	keys.push_back(at3);
490	keys.push_back(at4);
491	return torsionTypeCont_.find(keys);
492	}
493
494	InversionType* ForceField::getExactInversionType(const std::string &at1,
495	const std::string &at2,
496	const std::string &at3,
497	const std::string &at4){
498	std::vector<std::string> keys;
499	keys.push_back(at1);
500	keys.push_back(at2);
501	keys.push_back(at3);
502	keys.push_back(at4);
503	return inversionTypeCont_.find(keys);
504	}
505
506	NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
507	std::vector<std::string> keys;
508	keys.push_back(at1);
509	keys.push_back(at2);
510	return nonBondedInteractionTypeCont_.find(keys);
511	}
512
513
514	bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
515	std::vector<std::string> keys;
516	keys.push_back(at);
517	return atomTypeCont_.add(keys, atomType);
518	}
519
520	bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
521	std::vector<std::string> keys;
522	keys.push_back(at);
523	return atomTypeCont_.replace(keys, atomType);
524	}
525
526	bool ForceField::addBondType(const std::string &at1, const std::string &at2,
527	BondType* bondType) {
528	std::vector<std::string> keys;
529	keys.push_back(at1);
530	keys.push_back(at2);
531	return bondTypeCont_.add(keys, bondType);
532	}
533
534	bool ForceField::addBendType(const std::string &at1, const std::string &at2,
535	const std::string &at3, BendType* bendType) {
536	std::vector<std::string> keys;
537	keys.push_back(at1);
538	keys.push_back(at2);
539	keys.push_back(at3);
540	return bendTypeCont_.add(keys, bendType);
541	}
542
543	bool ForceField::addTorsionType(const std::string &at1,
544	const std::string &at2,
545	const std::string &at3,
546	const std::string &at4,
547	TorsionType* torsionType) {
548	std::vector<std::string> keys;
549	keys.push_back(at1);
550	keys.push_back(at2);
551	keys.push_back(at3);
552	keys.push_back(at4);
553	return torsionTypeCont_.add(keys, torsionType);
554	}
555
556	bool ForceField::addInversionType(const std::string &at1,
557	const std::string &at2,
558	const std::string &at3,
559	const std::string &at4,
560	InversionType* inversionType) {
561	std::vector<std::string> keys;
562	keys.push_back(at1);
563	keys.push_back(at2);
564	keys.push_back(at3);
565	keys.push_back(at4);
566	return inversionTypeCont_.add(keys, inversionType);
567	}
568
569	bool ForceField::addNonBondedInteractionType(const std::string &at1,
570	const std::string &at2,
571	NonBondedInteractionType* nbiType) {
572	std::vector<std::string> keys;
573	keys.push_back(at1);
574	keys.push_back(at2);
575	return nonBondedInteractionTypeCont_.add(keys, nbiType);
576	}
577
578	RealType ForceField::getRcutFromAtomType(AtomType* at) {
579	/*@todo /
580	GenericData* data;
581	RealType rcut = 0.0;
582
583	if (at->isLennardJones()) {
584	data = at->getPropertyByName("LennardJones");
585	if (data != NULL) {
586	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
587
588	if (ljData != NULL) {
589	LJParam ljParam = ljData->getData();
590
591	//by default use 2.5*sigma as cutoff radius
592	rcut = 2.5 * ljParam.sigma;
593
594	} else {
595	sprintf( painCave.errMsg,
596	"Can not cast GenericData to LJParam\n");
597	painCave.severity = OOPSE_ERROR;
598	painCave.isFatal = 1;
599	simError();
600	}
601	} else {
602	sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
603	painCave.severity = OOPSE_ERROR;
604	painCave.isFatal = 1;
605	simError();
606	}
607	}
608	return rcut;
609	}
610
611
612	ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
613	std::string forceFieldFilename(filename);
614	ifstrstream* ffStream = new ifstrstream();
615
616	//try to open the force filed file in current directory first
617	ffStream->open(forceFieldFilename.c_str());
618	if(!ffStream->is_open()){
619
620	forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
621	ffStream->open( forceFieldFilename.c_str() );
622
623	//if current directory does not contain the force field file,
624	//try to open it in the path
625	if(!ffStream->is_open()){
626
627	sprintf( painCave.errMsg,
628	"Error opening the force field parameter file:\n"
629	"\t%s\n"
630	"\tHave you tried setting the FORCE_PARAM_PATH environment "
631	"variable?\n",
632	forceFieldFilename.c_str() );
633	painCave.severity = OOPSE_ERROR;
634	painCave.isFatal = 1;
635	simError();
636	}
637	}
638	return ffStream;
639	}
640
641	void ForceField::setFortranForceOptions(){
642	ForceOptions theseFortranOptions;
643	forceFieldOptions_.makeFortranOptions(theseFortranOptions);
644	setfForceOptions(&theseFortranOptions);
645	}
646	} //end namespace oopse