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

          BondType* bondType = bondTypeCont_.find(myKeys);
          if (bondType) {
            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
          }
          jj++;
        }
        ii++;
      }

      // 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);
      if (bestType) 
        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++;
      }
      
      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);  
      if (bestType) 
        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++;
      }
      
      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);
      if (bestType) {
        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++;
      }
      
      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);
      if (bestType) {
        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::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:	1275
Committed:	Fri Jul 4 20:54:29 2008 UTC (16 years, 9 months ago) by cli2
Original Path:	*trunk/src/UseTheForce/ForceField.cpp*
File size:	19780 byte(s)
Log Message:	Changes required for Inversions and Base Atom types. This will break OOPSE badly for a few days or so...
#	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			* 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	gezelter	507	/**
43			* @file ForceField.cpp
44			* @author tlin
45			* @date 11/04/2004
46			* @time 22:51am
47			* @version 1.0
48			*/
49	gezelter	246
50	gezelter	1269	#include <algorithm>
51	gezelter	206	#include "UseTheForce/ForceField.hpp"
52	gezelter	246	#include "utils/simError.h"
53	gezelter	1269	#include "utils/Tuple.hpp"
54	tim	475	#include "UseTheForce/DarkSide/atype_interface.h"
55	chuckv	821	#include "UseTheForce/DarkSide/fForceOptions_interface.h"
56	gezelter	939	#include "UseTheForce/DarkSide/switcheroo_interface.h"
57	gezelter	246	namespace oopse {
58	gezelter	206
59	gezelter	507	ForceField::ForceField() {
60	gezelter	246	char* tempPath;
61			tempPath = getenv("FORCE_PARAM_PATH");
62	gezelter	206
63	gezelter	246	if (tempPath == NULL) {
64	gezelter	507	//convert a macro from compiler to a string in c++
65			STR_DEFINE(ffPath_, FRC_PATH );
66	gezelter	246	} else {
67	gezelter	507	ffPath_ = tempPath;
68	gezelter	246	}
69	gezelter	507	}
70	gezelter	206
71	tim	475
72	gezelter	507	ForceField::~ForceField() {
73	tim	475	deleteAtypes();
74	gezelter	939	deleteSwitch();
75	gezelter	507	}
76	tim	475
77	gezelter	507	AtomType* ForceField::getAtomType(const std::string &at) {
78	gezelter	246	std::vector<std::string> keys;
79			keys.push_back(at);
80			return atomTypeCont_.find(keys);
81	gezelter	507	}
82	gezelter	206
83	cpuglis	1195	BondType* ForceField::getBondType(const std::string &at1,
84			const std::string &at2) {
85	gezelter	246	std::vector<std::string> keys;
86			keys.push_back(at1);
87			keys.push_back(at2);
88	gezelter	206
89	gezelter	246	//try exact match first
90			BondType* bondType = bondTypeCont_.find(keys);
91			if (bondType) {
92	gezelter	507	return bondType;
93	gezelter	246	} else {
94	gezelter	1269	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			BondType* bondType = bondTypeCont_.find(myKeys);
128			if (bondType) {
129			foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
130			}
131			jj++;
132			}
133			ii++;
134			}
135
136			// sort the foundBonds by the score:
137
138			std::sort(foundBonds.begin(), foundBonds.end());
139
140			int bestScore = foundBonds[0].first;
141			std::vector<std::string> theKeys = foundBonds[0].second;
142
143			std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1] << "\t(score = "<< bestScore << ")\n";
144			BondType* bestType = bondTypeCont_.find(theKeys);
145			if (bestType)
146			return bestType;
147			else {
148			//if no exact match found, try wild card match
149			return bondTypeCont_.find(keys, wildCardAtomTypeName_);
150			}
151	gezelter	246	}
152	gezelter	507	}
153	gezelter	1269
154	cpuglis	1195	BendType* ForceField::getBendType(const std::string &at1,
155			const std::string &at2,
156	gezelter	507	const std::string &at3) {
157	gezelter	246	std::vector<std::string> keys;
158			keys.push_back(at1);
159			keys.push_back(at2);
160			keys.push_back(at3);
161	gezelter	206
162	gezelter	246	//try exact match first
163			BendType* bendType = bendTypeCont_.find(keys);
164			if (bendType) {
165	gezelter	507	return bendType;
166	gezelter	246	} else {
167	gezelter	1269
168			AtomType* atype1;
169			AtomType* atype2;
170			AtomType* atype3;
171			std::vector<std::string> at1key;
172			at1key.push_back(at1);
173			atype1 = atomTypeCont_.find(at1key);
174
175			std::vector<std::string> at2key;
176			at2key.push_back(at2);
177			atype2 = atomTypeCont_.find(at2key);
178
179			std::vector<std::string> at3key;
180			at3key.push_back(at3);
181			atype3 = atomTypeCont_.find(at3key);
182
183			// query atom types for their chains of responsibility
184			std::vector<AtomType*> at1Chain = atype1->allYourBase();
185			std::vector<AtomType*> at2Chain = atype2->allYourBase();
186			std::vector<AtomType*> at3Chain = atype3->allYourBase();
187
188			std::vector<AtomType*>::iterator i;
189			std::vector<AtomType*>::iterator j;
190			std::vector<AtomType*>::iterator k;
191
192			int ii = 0;
193			int jj = 0;
194			int kk = 0;
195			int IKscore;
196
197			std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
198
199			for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
200			ii = 0;
201			for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
202			kk = 0;
203			for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
204
205			IKscore = ii + kk;
206
207			std::vector<std::string> myKeys;
208			myKeys.push_back((*i)->getName());
209			myKeys.push_back((*j)->getName());
210			myKeys.push_back((*k)->getName());
211
212			BendType* bendType = bendTypeCont_.find(myKeys);
213			if (bendType) {
214			foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
215			}
216			kk++;
217			}
218			ii++;
219			}
220			jj++;
221			}
222
223			std::sort(foundBends.begin(), foundBends.end());
224
225			int jscore = foundBends[0].first;
226			int ikscore = foundBends[0].second;
227			std::vector<std::string> theKeys = foundBends[0].third;
228
229			std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";
230
231			BendType* bestType = bendTypeCont_.find(theKeys);
232			if (bestType)
233			return bestType;
234			else {
235
236			//if no exact match found, try wild card match
237			return bendTypeCont_.find(keys, wildCardAtomTypeName_);
238			}
239	gezelter	246	}
240	gezelter	507	}
241	gezelter	206
242	gezelter	1269
243	cpuglis	1195	TorsionType* ForceField::getTorsionType(const std::string &at1,
244			const std::string &at2,
245			const std::string &at3,
246			const std::string &at4) {
247	gezelter	246	std::vector<std::string> keys;
248			keys.push_back(at1);
249			keys.push_back(at2);
250			keys.push_back(at3);
251			keys.push_back(at4);
252	gezelter	206
253	gezelter	1269
254			//try exact match first
255	gezelter	246	TorsionType* torsionType = torsionTypeCont_.find(keys);
256			if (torsionType) {
257	gezelter	507	return torsionType;
258	gezelter	246	} else {
259	gezelter	1269
260			AtomType* atype1;
261			AtomType* atype2;
262			AtomType* atype3;
263			AtomType* atype4;
264			std::vector<std::string> at1key;
265			at1key.push_back(at1);
266			atype1 = atomTypeCont_.find(at1key);
267
268			std::vector<std::string> at2key;
269			at2key.push_back(at2);
270			atype2 = atomTypeCont_.find(at2key);
271
272			std::vector<std::string> at3key;
273			at3key.push_back(at3);
274			atype3 = atomTypeCont_.find(at3key);
275
276			std::vector<std::string> at4key;
277			at4key.push_back(at4);
278			atype4 = atomTypeCont_.find(at4key);
279
280			// query atom types for their chains of responsibility
281			std::vector<AtomType*> at1Chain = atype1->allYourBase();
282			std::vector<AtomType*> at2Chain = atype2->allYourBase();
283			std::vector<AtomType*> at3Chain = atype3->allYourBase();
284			std::vector<AtomType*> at4Chain = atype4->allYourBase();
285
286			std::vector<AtomType*>::iterator i;
287			std::vector<AtomType*>::iterator j;
288			std::vector<AtomType*>::iterator k;
289			std::vector<AtomType*>::iterator l;
290
291			int ii = 0;
292			int jj = 0;
293			int kk = 0;
294			int ll = 0;
295			int ILscore;
296			int JKscore;
297
298			std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
299
300			for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
301			kk = 0;
302			for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
303			ii = 0;
304			for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
305			ll = 0;
306			for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
307
308			ILscore = ii + ll;
309			JKscore = jj + kk;
310
311			std::vector<std::string> myKeys;
312			myKeys.push_back((*i)->getName());
313			myKeys.push_back((*j)->getName());
314			myKeys.push_back((*k)->getName());
315			myKeys.push_back((*l)->getName());
316
317			TorsionType* torsionType = torsionTypeCont_.find(myKeys);
318			if (torsionType) {
319			foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
320			}
321			ll++;
322			}
323			ii++;
324			}
325			kk++;
326			}
327			jj++;
328			}
329
330			std::sort(foundTorsions.begin(), foundTorsions.end());
331
332			int jkscore = foundTorsions[0].first;
333			int ilscore = foundTorsions[0].second;
334			std::vector<std::string> theKeys = foundTorsions[0].third;
335
336			std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
337
338
339			TorsionType* bestType = torsionTypeCont_.find(theKeys);
340			if (bestType) {
341			return bestType;
342			} 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			InversionType* inversionType = inversionTypeCont_.find(keys);
361			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			InversionType* inversionType = inversionTypeCont_.find(myKeys);
423			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
435			std::sort(foundInversions.begin(), foundInversions.end());
436
437			int iscore = foundInversions[0].first;
438			int jklscore = foundInversions[0].second;
439			std::vector<std::string> theKeys = foundInversions[0].third;
440
441			std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1] << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
442
443
444			InversionType* bestType = inversionTypeCont_.find(theKeys);
445			if (bestType) {
446			return bestType;
447			} else {
448			//if no exact match found, try wild card match
449			return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
450			}
451			}
452			}
453
454	chuckv	1151	NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
455			std::vector<std::string> keys;
456			keys.push_back(at1);
457			keys.push_back(at2);
458	cpuglis	1195
459	chuckv	1151	//try exact match first
460			NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
461			if (nbiType) {
462			return nbiType;
463			} else {
464			//if no exact match found, try wild card match
465			return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
466	cpuglis	1195	}
467	chuckv	1151	}
468	cpuglis	1195
469			BondType* ForceField::getExactBondType(const std::string &at1,
470			const std::string &at2){
471	gezelter	246	std::vector<std::string> keys;
472			keys.push_back(at1);
473			keys.push_back(at2);
474			return bondTypeCont_.find(keys);
475	gezelter	507	}
476	cpuglis	1195
477			BendType* ForceField::getExactBendType(const std::string &at1,
478			const std::string &at2,
479	gezelter	507	const std::string &at3){
480	gezelter	246	std::vector<std::string> keys;
481			keys.push_back(at1);
482			keys.push_back(at2);
483			keys.push_back(at3);
484			return bendTypeCont_.find(keys);
485	gezelter	507	}
486	cpuglis	1195
487			TorsionType* ForceField::getExactTorsionType(const std::string &at1,
488			const std::string &at2,
489			const std::string &at3,
490			const std::string &at4){
491	gezelter	246	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 torsionTypeCont_.find(keys);
497	gezelter	507	}
498	cli2	1275
499			InversionType* ForceField::getExactInversionType(const std::string &at1,
500			const std::string &at2,
501			const std::string &at3,
502			const std::string &at4){
503			std::vector<std::string> keys;
504			keys.push_back(at1);
505			keys.push_back(at2);
506			keys.push_back(at3);
507			keys.push_back(at4);
508			return inversionTypeCont_.find(keys);
509			}
510
511	chuckv	1151	NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
512			std::vector<std::string> keys;
513			keys.push_back(at1);
514			keys.push_back(at2);
515			return nonBondedInteractionTypeCont_.find(keys);
516			}
517	cli2	1275
518	chuckv	1151
519	gezelter	507	bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
520	gezelter	246	std::vector<std::string> keys;
521			keys.push_back(at);
522			return atomTypeCont_.add(keys, atomType);
523	gezelter	507	}
524	gezelter	206
525	cpuglis	1195	bool ForceField::addBondType(const std::string &at1, const std::string &at2,
526			BondType* bondType) {
527	gezelter	246	std::vector<std::string> keys;
528			keys.push_back(at1);
529			keys.push_back(at2);
530	cpuglis	1195	return bondTypeCont_.add(keys, bondType);
531	gezelter	507	}
532	cpuglis	1195
533	gezelter	507	bool ForceField::addBendType(const std::string &at1, const std::string &at2,
534			const std::string &at3, BendType* bendType) {
535	gezelter	246	std::vector<std::string> keys;
536			keys.push_back(at1);
537			keys.push_back(at2);
538			keys.push_back(at3);
539			return bendTypeCont_.add(keys, bendType);
540	gezelter	507	}
541	cpuglis	1195
542			bool ForceField::addTorsionType(const std::string &at1,
543			const std::string &at2,
544			const std::string &at3,
545			const std::string &at4,
546			TorsionType* torsionType) {
547	gezelter	246	std::vector<std::string> keys;
548			keys.push_back(at1);
549			keys.push_back(at2);
550			keys.push_back(at3);
551			keys.push_back(at4);
552			return torsionTypeCont_.add(keys, torsionType);
553	gezelter	507	}
554	gezelter	206
555	cli2	1275	bool ForceField::addInversionType(const std::string &at1,
556			const std::string &at2,
557			const std::string &at3,
558			const std::string &at4,
559			InversionType* inversionType) {
560			std::vector<std::string> keys;
561			keys.push_back(at1);
562			keys.push_back(at2);
563			keys.push_back(at3);
564			keys.push_back(at4);
565			return inversionTypeCont_.add(keys, inversionType);
566			}
567
568	cpuglis	1195	bool ForceField::addNonBondedInteractionType(const std::string &at1,
569			const std::string &at2,
570			NonBondedInteractionType* nbiType) {
571	chuckv	1151	std::vector<std::string> keys;
572			keys.push_back(at1);
573			keys.push_back(at2);
574			return nonBondedInteractionTypeCont_.add(keys, nbiType);
575			}
576	cpuglis	1195
577	tim	963	RealType ForceField::getRcutFromAtomType(AtomType* at) {
578	gezelter	246	/*@todo /
579			GenericData* data;
580	tim	963	RealType rcut = 0.0;
581	cpuglis	1195
582	gezelter	246	if (at->isLennardJones()) {
583	gezelter	507	data = at->getPropertyByName("LennardJones");
584			if (data != NULL) {
585			LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
586	cpuglis	1195
587	gezelter	507	if (ljData != NULL) {
588			LJParam ljParam = ljData->getData();
589	cpuglis	1195
590	gezelter	507	//by default use 2.5*sigma as cutoff radius
591			rcut = 2.5 * ljParam.sigma;
592	cpuglis	1195
593	gezelter	507	} else {
594			sprintf( painCave.errMsg,
595			"Can not cast GenericData to LJParam\n");
596			painCave.severity = OOPSE_ERROR;
597			painCave.isFatal = 1;
598			simError();
599			}
600			} else {
601			sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
602			painCave.severity = OOPSE_ERROR;
603			painCave.isFatal = 1;
604			simError();
605			}
606	gezelter	246	}
607			return rcut;
608	gezelter	507	}
609	cpuglis	1195
610	gezelter	206
611	gezelter	507	ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
612	gezelter	246	std::string forceFieldFilename(filename);
613			ifstrstream* ffStream = new ifstrstream();
614
615			//try to open the force filed file in current directory first
616			ffStream->open(forceFieldFilename.c_str());
617			if(!ffStream->is_open()){
618
619	gezelter	507	forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
620			ffStream->open( forceFieldFilename.c_str() );
621	gezelter	246
622	gezelter	507	//if current directory does not contain the force field file,
623			//try to open it in the path
624			if(!ffStream->is_open()){
625	gezelter	246
626	gezelter	507	sprintf( painCave.errMsg,
627			"Error opening the force field parameter file:\n"
628			"\t%s\n"
629			"\tHave you tried setting the FORCE_PARAM_PATH environment "
630			"variable?\n",
631			forceFieldFilename.c_str() );
632			painCave.severity = OOPSE_ERROR;
633			painCave.isFatal = 1;
634			simError();
635			}
636	gezelter	246	}
637			return ffStream;
638	gezelter	507	}
639	gezelter	246
640	chuckv	821	void ForceField::setFortranForceOptions(){
641			ForceOptions theseFortranOptions;
642			forceFieldOptions_.makeFortranOptions(theseFortranOptions);
643			setfForceOptions(&theseFortranOptions);
644			}
645	gezelter	246	} //end namespace oopse