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

          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 = 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:	1303
Committed:	Mon Oct 13 21:35:22 2008 UTC (16 years, 6 months ago) by cli2
File size:	19562 byte(s)
Log Message:	Fixes for Inversions for use with Amber
#	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	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
137	if (foundBonds.size() > 0) {
138	// sort the foundBonds by the score:
139	std::sort(foundBonds.begin(), foundBonds.end());
140
141	int bestScore = foundBonds[0].first;
142	std::vector<std::string> theKeys = foundBonds[0].second;
143
144	BondType* bestType = bondTypeCont_.find(theKeys);
145
146	return bestType;
147	} else {
148	//if no exact match found, try wild card match
149	return bondTypeCont_.find(keys, wildCardAtomTypeName_);
150	}
151	}
152	}
153
154	BendType* ForceField::getBendType(const std::string &at1,
155	const std::string &at2,
156	const std::string &at3) {
157	std::vector<std::string> keys;
158	keys.push_back(at1);
159	keys.push_back(at2);
160	keys.push_back(at3);
161
162	//try exact match first
163	BendType* bendType = bendTypeCont_.find(keys);
164	if (bendType) {
165	return bendType;
166	} else {
167
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	if (foundBends.size() > 0) {
224	std::sort(foundBends.begin(), foundBends.end());
225	int jscore = foundBends[0].first;
226	int ikscore = foundBends[0].second;
227	std::vector<std::string> theKeys = foundBends[0].third;
228
229	BendType* bestType = bendTypeCont_.find(theKeys);
230	return bestType;
231	} else {
232	//if no exact match found, try wild card match
233	return bendTypeCont_.find(keys, wildCardAtomTypeName_);
234	}
235	}
236	}
237
238	TorsionType* ForceField::getTorsionType(const std::string &at1,
239	const std::string &at2,
240	const std::string &at3,
241	const std::string &at4) {
242	std::vector<std::string> keys;
243	keys.push_back(at1);
244	keys.push_back(at2);
245	keys.push_back(at3);
246	keys.push_back(at4);
247
248
249	//try exact match first
250	TorsionType* torsionType = torsionTypeCont_.find(keys);
251	if (torsionType) {
252	return torsionType;
253	} else {
254
255	AtomType* atype1;
256	AtomType* atype2;
257	AtomType* atype3;
258	AtomType* atype4;
259	std::vector<std::string> at1key;
260	at1key.push_back(at1);
261	atype1 = atomTypeCont_.find(at1key);
262
263	std::vector<std::string> at2key;
264	at2key.push_back(at2);
265	atype2 = atomTypeCont_.find(at2key);
266
267	std::vector<std::string> at3key;
268	at3key.push_back(at3);
269	atype3 = atomTypeCont_.find(at3key);
270
271	std::vector<std::string> at4key;
272	at4key.push_back(at4);
273	atype4 = atomTypeCont_.find(at4key);
274
275	// query atom types for their chains of responsibility
276	std::vector<AtomType*> at1Chain = atype1->allYourBase();
277	std::vector<AtomType*> at2Chain = atype2->allYourBase();
278	std::vector<AtomType*> at3Chain = atype3->allYourBase();
279	std::vector<AtomType*> at4Chain = atype4->allYourBase();
280
281	std::vector<AtomType*>::iterator i;
282	std::vector<AtomType*>::iterator j;
283	std::vector<AtomType*>::iterator k;
284	std::vector<AtomType*>::iterator l;
285
286	int ii = 0;
287	int jj = 0;
288	int kk = 0;
289	int ll = 0;
290	int ILscore;
291	int JKscore;
292
293	std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
294
295	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
296	kk = 0;
297	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
298	ii = 0;
299	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
300	ll = 0;
301	for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
302
303	ILscore = ii + ll;
304	JKscore = jj + kk;
305
306	std::vector<std::string> myKeys;
307	myKeys.push_back((*i)->getName());
308	myKeys.push_back((*j)->getName());
309	myKeys.push_back((*k)->getName());
310	myKeys.push_back((*l)->getName());
311
312	TorsionType* torsionType = torsionTypeCont_.find(myKeys);
313	if (torsionType) {
314	foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
315	}
316	ll++;
317	}
318	ii++;
319	}
320	kk++;
321	}
322	jj++;
323	}
324
325	if (foundTorsions.size() > 0) {
326	std::sort(foundTorsions.begin(), foundTorsions.end());
327	int jkscore = foundTorsions[0].first;
328	int ilscore = foundTorsions[0].second;
329	std::vector<std::string> theKeys = foundTorsions[0].third;
330
331	TorsionType* bestType = torsionTypeCont_.find(theKeys);
332	return bestType;
333	} else {
334	//if no exact match found, try wild card match
335	return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
336	}
337	}
338	}
339
340	InversionType* ForceField::getInversionType(const std::string &at1,
341	const std::string &at2,
342	const std::string &at3,
343	const std::string &at4) {
344	std::vector<std::string> keys;
345	keys.push_back(at1);
346	keys.push_back(at2);
347	keys.push_back(at3);
348	keys.push_back(at4);
349
350	//try exact match first
351	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
352	if (inversionType) {
353	return inversionType;
354	} else {
355
356	AtomType* atype1;
357	AtomType* atype2;
358	AtomType* atype3;
359	AtomType* atype4;
360	std::vector<std::string> at1key;
361	at1key.push_back(at1);
362	atype1 = atomTypeCont_.find(at1key);
363
364	std::vector<std::string> at2key;
365	at2key.push_back(at2);
366	atype2 = atomTypeCont_.find(at2key);
367
368	std::vector<std::string> at3key;
369	at3key.push_back(at3);
370	atype3 = atomTypeCont_.find(at3key);
371
372	std::vector<std::string> at4key;
373	at4key.push_back(at4);
374	atype4 = atomTypeCont_.find(at4key);
375
376	// query atom types for their chains of responsibility
377	std::vector<AtomType*> at1Chain = atype1->allYourBase();
378	std::vector<AtomType*> at2Chain = atype2->allYourBase();
379	std::vector<AtomType*> at3Chain = atype3->allYourBase();
380	std::vector<AtomType*> at4Chain = atype4->allYourBase();
381
382	std::vector<AtomType*>::iterator i;
383	std::vector<AtomType*>::iterator j;
384	std::vector<AtomType*>::iterator k;
385	std::vector<AtomType*>::iterator l;
386
387	int ii = 0;
388	int jj = 0;
389	int kk = 0;
390	int ll = 0;
391	int Iscore;
392	int JKLscore;
393
394	std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
395
396	for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
397	kk = 0;
398	for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
399	ii = 0;
400	for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
401	ll = 0;
402	for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
403
404	Iscore = ii;
405	JKLscore = jj + kk + ll;
406
407	std::vector<std::string> myKeys;
408	myKeys.push_back((*i)->getName());
409	myKeys.push_back((*j)->getName());
410	myKeys.push_back((*k)->getName());
411	myKeys.push_back((*l)->getName());
412
413	InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
414	if (inversionType) {
415	foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
416	}
417	ll++;
418	}
419	ii++;
420	}
421	kk++;
422	}
423	jj++;
424	}
425
426	if (foundInversions.size() > 0) {
427	std::sort(foundInversions.begin(), foundInversions.end());
428	int iscore = foundInversions[0].first;
429	int jklscore = foundInversions[0].second;
430	std::vector<std::string> theKeys = foundInversions[0].third;
431
432	InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
433	return bestType;
434	} else {
435	//if no exact match found, try wild card match
436	return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
437	}
438	}
439	}
440
441	NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
442	std::vector<std::string> keys;
443	keys.push_back(at1);
444	keys.push_back(at2);
445
446	//try exact match first
447	NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
448	if (nbiType) {
449	return nbiType;
450	} else {
451	//if no exact match found, try wild card match
452	return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
453	}
454	}
455
456	BondType* ForceField::getExactBondType(const std::string &at1,
457	const std::string &at2){
458	std::vector<std::string> keys;
459	keys.push_back(at1);
460	keys.push_back(at2);
461	return bondTypeCont_.find(keys);
462	}
463
464	BendType* ForceField::getExactBendType(const std::string &at1,
465	const std::string &at2,
466	const std::string &at3){
467	std::vector<std::string> keys;
468	keys.push_back(at1);
469	keys.push_back(at2);
470	keys.push_back(at3);
471	return bendTypeCont_.find(keys);
472	}
473
474	TorsionType* ForceField::getExactTorsionType(const std::string &at1,
475	const std::string &at2,
476	const std::string &at3,
477	const std::string &at4){
478	std::vector<std::string> keys;
479	keys.push_back(at1);
480	keys.push_back(at2);
481	keys.push_back(at3);
482	keys.push_back(at4);
483	return torsionTypeCont_.find(keys);
484	}
485
486	InversionType* ForceField::getExactInversionType(const std::string &at1,
487	const std::string &at2,
488	const std::string &at3,
489	const std::string &at4){
490	std::vector<std::string> keys;
491	keys.push_back(at1);
492	keys.push_back(at2);
493	keys.push_back(at3);
494	keys.push_back(at4);
495	return inversionTypeCont_.find(keys);
496	}
497
498	NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
499	std::vector<std::string> keys;
500	keys.push_back(at1);
501	keys.push_back(at2);
502	return nonBondedInteractionTypeCont_.find(keys);
503	}
504
505
506	bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
507	std::vector<std::string> keys;
508	keys.push_back(at);
509	return atomTypeCont_.add(keys, atomType);
510	}
511
512	bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
513	std::vector<std::string> keys;
514	keys.push_back(at);
515	return atomTypeCont_.replace(keys, atomType);
516	}
517
518	bool ForceField::addBondType(const std::string &at1, const std::string &at2,
519	BondType* bondType) {
520	std::vector<std::string> keys;
521	keys.push_back(at1);
522	keys.push_back(at2);
523	return bondTypeCont_.add(keys, bondType);
524	}
525
526	bool ForceField::addBendType(const std::string &at1, const std::string &at2,
527	const std::string &at3, BendType* bendType) {
528	std::vector<std::string> keys;
529	keys.push_back(at1);
530	keys.push_back(at2);
531	keys.push_back(at3);
532	return bendTypeCont_.add(keys, bendType);
533	}
534
535	bool ForceField::addTorsionType(const std::string &at1,
536	const std::string &at2,
537	const std::string &at3,
538	const std::string &at4,
539	TorsionType* torsionType) {
540	std::vector<std::string> keys;
541	keys.push_back(at1);
542	keys.push_back(at2);
543	keys.push_back(at3);
544	keys.push_back(at4);
545	return torsionTypeCont_.add(keys, torsionType);
546	}
547
548	bool ForceField::addInversionType(const std::string &at1,
549	const std::string &at2,
550	const std::string &at3,
551	const std::string &at4,
552	InversionType* inversionType) {
553	std::vector<std::string> keys;
554	keys.push_back(at1);
555	keys.push_back(at2);
556	keys.push_back(at3);
557	keys.push_back(at4);
558	return inversionTypeCont_.add(keys, inversionType);
559	}
560
561	bool ForceField::addNonBondedInteractionType(const std::string &at1,
562	const std::string &at2,
563	NonBondedInteractionType* nbiType) {
564	std::vector<std::string> keys;
565	keys.push_back(at1);
566	keys.push_back(at2);
567	return nonBondedInteractionTypeCont_.add(keys, nbiType);
568	}
569
570	RealType ForceField::getRcutFromAtomType(AtomType* at) {
571	/*@todo /
572	GenericData* data;
573	RealType rcut = 0.0;
574
575	if (at->isLennardJones()) {
576	data = at->getPropertyByName("LennardJones");
577	if (data != NULL) {
578	LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
579
580	if (ljData != NULL) {
581	LJParam ljParam = ljData->getData();
582
583	//by default use 2.5*sigma as cutoff radius
584	rcut = 2.5 * ljParam.sigma;
585
586	} else {
587	sprintf( painCave.errMsg,
588	"Can not cast GenericData to LJParam\n");
589	painCave.severity = OOPSE_ERROR;
590	painCave.isFatal = 1;
591	simError();
592	}
593	} else {
594	sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
595	painCave.severity = OOPSE_ERROR;
596	painCave.isFatal = 1;
597	simError();
598	}
599	}
600	return rcut;
601	}
602
603
604	ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
605	std::string forceFieldFilename(filename);
606	ifstrstream* ffStream = new ifstrstream();
607
608	//try to open the force filed file in current directory first
609	ffStream->open(forceFieldFilename.c_str());
610	if(!ffStream->is_open()){
611
612	forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
613	ffStream->open( forceFieldFilename.c_str() );
614
615	//if current directory does not contain the force field file,
616	//try to open it in the path
617	if(!ffStream->is_open()){
618
619	sprintf( painCave.errMsg,
620	"Error opening the force field parameter file:\n"
621	"\t%s\n"
622	"\tHave you tried setting the FORCE_PARAM_PATH environment "
623	"variable?\n",
624	forceFieldFilename.c_str() );
625	painCave.severity = OOPSE_ERROR;
626	painCave.isFatal = 1;
627	simError();
628	}
629	}
630	return ffStream;
631	}
632
633	void ForceField::setFortranForceOptions(){
634	ForceOptions theseFortranOptions;
635	forceFieldOptions_.makeFortranOptions(theseFortranOptions);
636	setfForceOptions(&theseFortranOptions);
637	}
638	} //end namespace oopse