src/UseTheForce/ForceField.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Vardeman & Gezelter, in progress (2009).                        
 */
 
/**
 * @file ForceField.cpp
 * @author tlin
 * @date 11/04/2004
 * @time 22:51am
 * @version 1.0
 */
  
#include <algorithm>
#include "UseTheForce/ForceField.hpp"
#include "utils/simError.h"
#include "utils/Tuple.hpp"
#include "UseTheForce/DarkSide/atype_interface.h"
#include "UseTheForce/DarkSide/fForceOptions_interface.h"
#include "UseTheForce/DarkSide/switcheroo_interface.h"
namespace OpenMD {

  ForceField::ForceField() { 

    char* tempPath; 
    tempPath = getenv("FORCE_PARAM_PATH");
    
    if (tempPath == NULL) {
      //convert a macro from compiler to a string in c++
      STR_DEFINE(ffPath_, FRC_PATH );
    } else {
      ffPath_ = tempPath;
    }
  }


  ForceField::~ForceField() {
    deleteAtypes();
    deleteSwitch();
  }

  AtomType* ForceField::getAtomType(const std::string &at) {
    std::vector<std::string> keys;
    keys.push_back(at);
    return atomTypeCont_.find(keys);
  }

  BondType* ForceField::getBondType(const std::string &at1, 
                                    const std::string &at2) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    

    //try exact match first
    BondType* bondType = bondTypeCont_.find(keys);
    if (bondType) {
      return bondType;
    } else {
      AtomType* atype1;
      AtomType* atype2;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
  
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;

      int ii = 0;
      int jj = 0;
      int bondTypeScore;

      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;

      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
        jj = 0;
        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {

          bondTypeScore = ii + jj;

          std::vector<std::string> myKeys;
          myKeys.push_back((*i)->getName());
          myKeys.push_back((*j)->getName());

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


      if (foundBonds.size() > 0) {
        // sort the foundBonds by the score:
        std::sort(foundBonds.begin(), foundBonds.end());
     
        int bestScore = foundBonds[0].first;
        std::vector<std::string> theKeys = foundBonds[0].second;
        
        BondType* bestType = bondTypeCont_.find(theKeys);
        
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
      }
    }
  }
  
  BendType* ForceField::getBendType(const std::string &at1, 
                                    const std::string &at2,
                                    const std::string &at3) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    

    //try exact match first
    BendType* bendType = bendTypeCont_.find(keys);
    if (bendType) {
      return bendType;
    } else {

      AtomType* atype1;
      AtomType* atype2;
      AtomType* atype3;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
  
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);

      std::vector<std::string> at3key;
      at3key.push_back(at3);
      atype3 = atomTypeCont_.find(at3key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      std::vector<AtomType*> at3Chain = atype3->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      std::vector<AtomType*>::iterator k;

      int ii = 0;
      int jj = 0;
      int kk = 0;
      int IKscore;

      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;

      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
        ii = 0;
        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
          kk = 0;
          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
          
            IKscore = ii + kk;

            std::vector<std::string> myKeys;
            myKeys.push_back((*i)->getName());
            myKeys.push_back((*j)->getName());
            myKeys.push_back((*k)->getName());

            BendType* bendType = bendTypeCont_.find(myKeys);
            if (bendType) { 
              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
            }
            kk++;
          }
          ii++;
        }
        jj++;
      }
      
      if (foundBends.size() > 0) {
        std::sort(foundBends.begin(), foundBends.end());
        int jscore = foundBends[0].first;
        int ikscore = foundBends[0].second;
        std::vector<std::string> theKeys = foundBends[0].third;       
        
        BendType* bestType = bendTypeCont_.find(theKeys);  
        return bestType;
      } else {        
        //if no exact match found, try wild card match
        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
      }
    }
  }

  TorsionType* ForceField::getTorsionType(const std::string &at1, 
                                          const std::string &at2,
                                          const std::string &at3, 
                                          const std::string &at4) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    


    //try exact match first
    TorsionType* torsionType = torsionTypeCont_.find(keys);
    if (torsionType) {
      return torsionType;
    } else {

      AtomType* atype1;
      AtomType* atype2;
      AtomType* atype3;
      AtomType* atype4;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
  
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);

      std::vector<std::string> at3key;
      at3key.push_back(at3);
      atype3 = atomTypeCont_.find(at3key);

      std::vector<std::string> at4key;
      at4key.push_back(at4);
      atype4 = atomTypeCont_.find(at4key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      std::vector<AtomType*> at3Chain = atype3->allYourBase();
      std::vector<AtomType*> at4Chain = atype4->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      std::vector<AtomType*>::iterator k;
      std::vector<AtomType*>::iterator l;

      int ii = 0;
      int jj = 0;
      int kk = 0;
      int ll = 0;
      int ILscore;
      int JKscore;

      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;

      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
        kk = 0;
        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
          ii = 0;       
          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
            ll = 0;
            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
          
              ILscore = ii + ll;
              JKscore = jj + kk;

              std::vector<std::string> myKeys;
              myKeys.push_back((*i)->getName());
              myKeys.push_back((*j)->getName());
              myKeys.push_back((*k)->getName());
              myKeys.push_back((*l)->getName());

              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
              if (torsionType) { 
                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
              }
              ll++;
            }
            ii++;
          }
          kk++;
        }
        jj++;
      }
      
      if (foundTorsions.size() > 0) {
        std::sort(foundTorsions.begin(), foundTorsions.end());
        int jkscore = foundTorsions[0].first;
        int ilscore = foundTorsions[0].second;
        std::vector<std::string> theKeys = foundTorsions[0].third;
        
        TorsionType* bestType = torsionTypeCont_.find(theKeys);
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
      }
    }
  }

  InversionType* ForceField::getInversionType(const std::string &at1, 
                                              const std::string &at2,
                                              const std::string &at3, 
                                              const std::string &at4) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    

    //try exact match first
    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
    if (inversionType) {
      return inversionType;
    } else {
      
      AtomType* atype1;
      AtomType* atype2;
      AtomType* atype3;
      AtomType* atype4;
      std::vector<std::string> at1key;
      at1key.push_back(at1);
      atype1 = atomTypeCont_.find(at1key);
      
      std::vector<std::string> at2key;
      at2key.push_back(at2);
      atype2 = atomTypeCont_.find(at2key);
      
      std::vector<std::string> at3key;
      at3key.push_back(at3);
      atype3 = atomTypeCont_.find(at3key);
      
      std::vector<std::string> at4key;
      at4key.push_back(at4);
      atype4 = atomTypeCont_.find(at4key);

      // query atom types for their chains of responsibility
      std::vector<AtomType*> at1Chain = atype1->allYourBase();
      std::vector<AtomType*> at2Chain = atype2->allYourBase();
      std::vector<AtomType*> at3Chain = atype3->allYourBase();
      std::vector<AtomType*> at4Chain = atype4->allYourBase();

      std::vector<AtomType*>::iterator i;
      std::vector<AtomType*>::iterator j;
      std::vector<AtomType*>::iterator k;
      std::vector<AtomType*>::iterator l;

      int ii = 0;
      int jj = 0;
      int kk = 0;
      int ll = 0;
      int Iscore;
      int JKLscore;
      
      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
      
      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
        kk = 0;
        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
          ii = 0;       
          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
            ll = 0;
            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
              
              Iscore = ii;
              JKLscore = jj + kk + ll;
              
              std::vector<std::string> myKeys;
              myKeys.push_back((*i)->getName());
              myKeys.push_back((*j)->getName());
              myKeys.push_back((*k)->getName());
              myKeys.push_back((*l)->getName());
              
              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
              if (inversionType) { 
                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
              }
              ll++;
            }
            ii++;
          }
          kk++;
        }
        jj++;
      }
         
      if (foundInversions.size() > 0) {
        std::sort(foundInversions.begin(), foundInversions.end());
        int iscore = foundInversions[0].first;
        int jklscore = foundInversions[0].second;
        std::vector<std::string> theKeys = foundInversions[0].third;
        
        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
      }
    }
  }
  
  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    
    //try exact match first
    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
    if (nbiType) {
      return nbiType;
    } else {
      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 nbiTypeScore;

      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;

      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
        jj = 0;
        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {

          nbiTypeScore = ii + jj;

          std::vector<std::string> myKeys;
          myKeys.push_back((*i)->getName());
          myKeys.push_back((*j)->getName());

          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
          if (nbiType) {
            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
          }
          jj++;
        }
        ii++;
      }


      if (foundNBI.size() > 0) {
        // sort the foundNBI by the score:
        std::sort(foundNBI.begin(), foundNBI.end());
     
        int bestScore = foundNBI[0].first;
        std::vector<std::string> theKeys = foundNBI[0].second;
        
        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
        return bestType;
      } else {
        //if no exact match found, try wild card match
        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
      }
    }
  }
  
  BondType* ForceField::getExactBondType(const std::string &at1, 
                                         const std::string &at2){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return bondTypeCont_.find(keys);
  }
  
  BendType* ForceField::getExactBendType(const std::string &at1, 
                                         const std::string &at2,
                                         const std::string &at3){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    return bendTypeCont_.find(keys);
  }
  
  TorsionType* ForceField::getExactTorsionType(const std::string &at1, 
                                               const std::string &at2,
                                               const std::string &at3, 
                                               const std::string &at4){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);   
    return torsionTypeCont_.find(keys);
  }
  
  InversionType* ForceField::getExactInversionType(const std::string &at1, 
                                                   const std::string &at2,
                                                   const std::string &at3, 
                                                   const std::string &at4){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);   
    return inversionTypeCont_.find(keys);
  }
  
  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ 
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return nonBondedInteractionTypeCont_.find(keys);
  }
  

  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
    std::vector<std::string> keys;
    keys.push_back(at);
    return atomTypeCont_.add(keys, atomType);
  }

  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
    std::vector<std::string> keys;
    keys.push_back(at);
    return atomTypeCont_.replace(keys, atomType);
  }

  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
                               BondType* bondType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return bondTypeCont_.add(keys, bondType);    
  }
  
  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
                               const std::string &at3, BendType* bendType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    return bendTypeCont_.add(keys, bendType);
  }
  
  bool ForceField::addTorsionType(const std::string &at1, 
                                  const std::string &at2,
                                  const std::string &at3, 
                                  const std::string &at4, 
                                  TorsionType* torsionType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    
    return torsionTypeCont_.add(keys, torsionType);
  }

  bool ForceField::addInversionType(const std::string &at1, 
                                    const std::string &at2,
                                    const std::string &at3, 
                                    const std::string &at4, 
                                    InversionType* inversionType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    keys.push_back(at3);    
    keys.push_back(at4);    
    return inversionTypeCont_.add(keys, inversionType);
  }
  
  bool ForceField::addNonBondedInteractionType(const std::string &at1, 
                                               const std::string &at2, 
                                               NonBondedInteractionType* nbiType) {
    std::vector<std::string> keys;
    keys.push_back(at1);
    keys.push_back(at2);    
    return nonBondedInteractionTypeCont_.add(keys, nbiType);
  }
  
  RealType ForceField::getRcutFromAtomType(AtomType* at) {
    /**@todo */
    GenericData* data;
    RealType rcut = 0.0;
    
    if (at->isLennardJones()) {
      data = at->getPropertyByName("LennardJones");
      if (data != NULL) {
        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
        
        if (ljData != NULL) {
          LJParam ljParam = ljData->getData();
          
          //by default use 2.5*sigma as cutoff radius
          rcut = 2.5 * ljParam.sigma;
          
        } else {
          sprintf( painCave.errMsg,
                   "Can not cast GenericData to LJParam\n");
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();          
        }            
      } else {
        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();          
      }
    }
    return rcut;    
  }
  

  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
    std::string forceFieldFilename(filename);
    ifstrstream* ffStream = new ifstrstream();
    
    //try to open the force filed file in current directory first    
    ffStream->open(forceFieldFilename.c_str());
    if(!ffStream->is_open()){

      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
      ffStream->open( forceFieldFilename.c_str() );

      //if current directory does not contain the force field file,
      //try to open it in the path        
      if(!ffStream->is_open()){

        sprintf( painCave.errMsg,
                 "Error opening the force field parameter file:\n"
                 "\t%s\n"
                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
                 "variable?\n",
                 forceFieldFilename.c_str() );
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();
      }
    }  
    return ffStream;
  }

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