src/nonbonded/InteractionManager.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).                        
 */
 
#include "nonbonded/InteractionManager.hpp"
#include "UseTheForce/doForces_interface.h"

namespace OpenMD {

  InteractionManager* InteractionManager::_instance = NULL; 
  SimInfo* InteractionManager::info_ = NULL;
  bool InteractionManager::initialized_ = false;

  RealType InteractionManager::rCut_ = 0.0;
  RealType InteractionManager::rSwitch_ = 0.0;
  RealType InteractionManager::skinThickness_ = 0.0;
  RealType InteractionManager::listRadius_ = 0.0;
  CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE;
  SwitchingFunctionType InteractionManager::sft_ = cubic;
  RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0};
  RealType InteractionManager::electrostaticScale_[4] = {1.0, 0.0, 0.0, 0.0};

  map<int, AtomType*> InteractionManager::typeMap_;
  map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_;

  LJ* InteractionManager::lj_ = new LJ();
  GB* InteractionManager::gb_ = new GB();
  Sticky* InteractionManager::sticky_ = new Sticky();
  Morse* InteractionManager::morse_ = new Morse();
  EAM* InteractionManager::eam_ = new EAM();
  SC* InteractionManager::sc_ = new SC();
  Electrostatic* InteractionManager::electrostatic_ = new Electrostatic();
  MAW* InteractionManager::maw_ = new MAW();
  SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
 
  InteractionManager* InteractionManager::Instance() {
    if (!_instance) {
      _instance = new InteractionManager();
    }
    return _instance;
  }

  void InteractionManager::initialize() {
    
    ForceField* forceField_ = info_->getForceField();
    
    lj_->setForceField(forceField_);
    gb_->setForceField(forceField_);
    sticky_->setForceField(forceField_);
    eam_->setForceField(forceField_);
    sc_->setForceField(forceField_);
    morse_->setForceField(forceField_);
    electrostatic_->setForceField(forceField_);
    maw_->setForceField(forceField_);

    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();

    // Force fields can set options on how to scale van der Waals and electrostatic
    // interactions for atoms connected via bonds, bends and torsions
    // in this case the topological distance between atoms is:
    // 0 = topologically unconnected
    // 1 = bonded together 
    // 2 = connected via a bend
    // 3 = connected via a torsion

    vdwScale_[0] = 1.0;
    vdwScale_[1] = fopts.getvdw12scale();
    vdwScale_[2] = fopts.getvdw13scale();
    vdwScale_[3] = fopts.getvdw14scale();

    electrostaticScale_[0] = 1.0;
    electrostaticScale_[1] = fopts.getelectrostatic12scale();
    electrostaticScale_[2] = fopts.getelectrostatic13scale();
    electrostaticScale_[3] = fopts.getelectrostatic14scale();    

    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
    ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
    AtomType* atype1;
    AtomType* atype2;
    pair<AtomType*, AtomType*> key;
    pair<set<NonBondedInteraction*>::iterator, bool> ret;
    
    for (atype1 = atomTypes->beginType(i1); atype1 != NULL; 
         atype1 = atomTypes->nextType(i1)) {
      
      // add it to the map:
      AtomTypeProperties atp = atype1->getATP();    
      
      pair<map<int,AtomType*>::iterator,bool> ret;    
      ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
      if (ret.second == false) {
        sprintf( painCave.errMsg,
                 "InteractionManager already had a previous entry with ident %d\n",
                 atp.ident);
        painCave.severity = OPENMD_INFO;
        painCave.isFatal = 0;
        simError();                 
      }
    }
    
    // Now, iterate over all known types and add to the interaction map:
    
    map<int, AtomType*>::iterator it1, it2;
    for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
      atype1 = (*it1).second;

      for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {        
        atype2 = (*it2).second;
        
        bool vdwExplicit = false;
        bool metExplicit = false;
        bool hbExplicit = false;
                       
        key = make_pair(atype1, atype2);
        
        if (atype1->isLennardJones() && atype2->isLennardJones()) {
          interactions_[key].insert(lj_);
        }
        if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
          interactions_[key].insert(electrostatic_);
        }
        if (atype1->isSticky() && atype2->isSticky() ) {
          interactions_[key].insert(sticky_);
        }
        if (atype1->isStickyPower() && atype2->isStickyPower() ) {
          interactions_[key].insert(sticky_);
        }
        if (atype1->isEAM() && atype2->isEAM() ) {
          interactions_[key].insert(eam_);
        }
        if (atype1->isSC() && atype2->isSC() ) {
          interactions_[key].insert(sc_);
        }
        if (atype1->isGayBerne() && atype2->isGayBerne() ) {
          interactions_[key].insert(gb_);
        }
        if ((atype1->isGayBerne() && atype2->isLennardJones())
            || (atype1->isLennardJones() && atype2->isGayBerne())) {
          interactions_[key].insert(gb_);
        } 
        
        // look for an explicitly-set non-bonded interaction type using the 
        // two atom types.
        NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());
        
        if (nbiType != NULL) {

          if (nbiType->isLennardJones()) {
            // We found an explicit Lennard-Jones interaction.  
            // override all other vdw entries for this pair of atom types:
            set<NonBondedInteraction*>::iterator it;
            for (it = interactions_[key].begin(); 
                 it != interactions_[key].end(); ++it) {
              InteractionFamily ifam = (*it)->getFamily();
              if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
            }
            interactions_[key].insert(lj_);
            vdwExplicit = true;
          }
          
          if (nbiType->isMorse()) {
            if (vdwExplicit) {
              sprintf( painCave.errMsg,
                       "InteractionManager::initialize found more than one "
                       "explicit \n"
                       "\tvan der Waals interaction for atom types %s - %s\n",
                       atype1->getName().c_str(), atype2->getName().c_str());
              painCave.severity = OPENMD_ERROR;
              painCave.isFatal = 1;
              simError();
            }
            // We found an explicit Morse interaction.  
            // override all other vdw entries for this pair of atom types:
            set<NonBondedInteraction*>::iterator it;
            for (it = interactions_[key].begin(); 
                 it != interactions_[key].end(); ++it) {
              InteractionFamily ifam = (*it)->getFamily();
              if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
            }
            interactions_[key].insert(morse_);
            vdwExplicit = true;
          }
          
          if (nbiType->isEAM()) {
            // We found an explicit EAM interaction.  
            // override all other metallic entries for this pair of atom types:
            set<NonBondedInteraction*>::iterator it;
            for (it = interactions_[key].begin(); 
                 it != interactions_[key].end(); ++it) {
              InteractionFamily ifam = (*it)->getFamily();
              if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
            }
            interactions_[key].insert(eam_);
            metExplicit = true;
          }
          
          if (nbiType->isSC()) {
            if (metExplicit) {
              sprintf( painCave.errMsg,
                       "InteractionManager::initialize found more than one "
                       "explicit\n"
                       "\tmetallic interaction for atom types %s - %s\n",
                       atype1->getName().c_str(), atype2->getName().c_str());
              painCave.severity = OPENMD_ERROR;
              painCave.isFatal = 1;
              simError();
            }
            // We found an explicit Sutton-Chen interaction.  
            // override all other metallic entries for this pair of atom types:
            set<NonBondedInteraction*>::iterator it;
            for (it = interactions_[key].begin(); 
                 it != interactions_[key].end(); ++it) {
              InteractionFamily ifam = (*it)->getFamily();
              if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
            }
            interactions_[key].insert(sc_);
            metExplicit = true;
          }
          
          if (nbiType->isMAW()) {
            if (vdwExplicit) {
              sprintf( painCave.errMsg,
                       "InteractionManager::initialize found more than one "
                       "explicit\n"
                       "\tvan der Waals interaction for atom types %s - %s\n",
                       atype1->getName().c_str(), atype2->getName().c_str());
              painCave.severity = OPENMD_ERROR;
              painCave.isFatal = 1;
              simError();
            }
            // We found an explicit MAW interaction.  
            // override all other vdw entries for this pair of atom types:
            set<NonBondedInteraction*>::iterator it;
            for (it = interactions_[key].begin(); 
                 it != interactions_[key].end(); ++it) {
              InteractionFamily ifam = (*it)->getFamily();
              if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
            }
            interactions_[key].insert(maw_);
            vdwExplicit = true;
          }        
        }
      }
    }
    
    
    // make sure every pair of atom types in this simulation has a
    // non-bonded interaction:

    set<AtomType*> simTypes = info_->getSimulatedAtomTypes();
    set<AtomType*>::iterator it, jt;
    for (it = simTypes.begin(); it != simTypes.end(); ++it) {
      atype1 = (*it);
      for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) {
        atype2 = (*jt);
        key = make_pair(atype1, atype2);
        
        if (interactions_[key].size() == 0) {
          sprintf( painCave.errMsg,
                   "InteractionManager unable to find an appropriate non-bonded\n"
                   "\tinteraction for atom types %s - %s\n",
                   atype1->getName().c_str(), atype2->getName().c_str());
          painCave.severity = OPENMD_INFO;
          painCave.isFatal = 1;
          simError();
        }
      }
    }

    setupCutoffs();
    setupSwitching();
    setupNeighborlists();

    //int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0;
    //int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0;
    //notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf);
    //notifyFortranSkinThickness(&skinThickness_);

    initialized_ = true;
  }
  
  /**
   * setupCutoffs
   *
   * Sets the values of cutoffRadius and cutoffMethod
   *
   * cutoffRadius : realType
   *  If the cutoffRadius was explicitly set, use that value.
   *  If the cutoffRadius was not explicitly set:
   *      Are there electrostatic atoms?  Use 12.0 Angstroms.
   *      No electrostatic atoms?  Poll the atom types present in the
   *      simulation for suggested cutoff values (e.g. 2.5 * sigma).
   *      Use the maximum suggested value that was found.
   *
   * cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
   *      If cutoffMethod was explicitly set, use that choice.
   *      If cutoffMethod was not explicitly set, use SHIFTED_FORCE
   */
  void InteractionManager::setupCutoffs() {
    
    Globals* simParams_ = info_->getSimParams();
    
    if (simParams_->haveCutoffRadius()) {
      rCut_ = simParams_->getCutoffRadius();
    } else {      
      if (info_->usesElectrostaticAtoms()) {
        sprintf(painCave.errMsg,
                "InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n"
                "\tOpenMD will use a default value of 12.0 angstroms"
                "\tfor the cutoffRadius.\n");
        painCave.isFatal = 0;
        painCave.severity = OPENMD_INFO;
        simError();
        rCut_ = 12.0;
      } else {
        RealType thisCut;
        set<AtomType*>::iterator i;
        set<AtomType*> atomTypes;
        atomTypes = info_->getSimulatedAtomTypes();        
        for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
          thisCut = getSuggestedCutoffRadius((*i));
          rCut_ = max(thisCut, rCut_);
        }
        sprintf(painCave.errMsg,
                "InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n"
                "\tOpenMD will use %lf angstroms.\n",
                rCut_);
        painCave.isFatal = 0;
        painCave.severity = OPENMD_INFO;
        simError();
      }             
    }

    map<string, CutoffMethod> stringToCutoffMethod;
    stringToCutoffMethod["HARD"] = HARD;
    stringToCutoffMethod["SWITCHED"] = SWITCHED;
    stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;    
    stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
  
    if (simParams_->haveCutoffMethod()) {
      string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
      map<string, CutoffMethod>::iterator i;
      i = stringToCutoffMethod.find(cutMeth);
      if (i == stringToCutoffMethod.end()) {
        sprintf(painCave.errMsg,
                "InteractionManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
                "\tShould be one of: "
                "HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
                cutMeth.c_str());
        painCave.isFatal = 1;
        painCave.severity = OPENMD_ERROR;
        simError();
      } else {
        cutoffMethod_ = i->second;
      }
    } else {
      sprintf(painCave.errMsg,
              "InteractionManager::setupCutoffs: No value was set for the cutoffMethod.\n"
              "\tOpenMD will use SHIFTED_FORCE.\n");
        painCave.isFatal = 0;
        painCave.severity = OPENMD_INFO;
        simError();
        cutoffMethod_ = SHIFTED_FORCE;        
    }
  }


  /**
   * setupSwitching
   *
   * Sets the values of switchingRadius and 
   *  If the switchingRadius was explicitly set, use that value (but check it)
   *  If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
   */
  void InteractionManager::setupSwitching() {
    Globals* simParams_ = info_->getSimParams();

    if (simParams_->haveSwitchingRadius()) {
      rSwitch_ = simParams_->getSwitchingRadius();
      if (rSwitch_ > rCut_) {        
        sprintf(painCave.errMsg,
                "InteractionManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
                rSwitch_, rCut_);
        painCave.isFatal = 1;
        painCave.severity = OPENMD_ERROR;
        simError();
      }
    } else {      
      rSwitch_ = 0.85 * rCut_;
      sprintf(painCave.errMsg,
              "InteractionManager::setupSwitching: No value was set for the switchingRadius.\n"
              "\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
              "\tswitchingRadius = %f. for this simulation\n", rSwitch_);
      painCave.isFatal = 0;
      painCave.severity = OPENMD_WARNING;
      simError();
    }           
    
    if (simParams_->haveSwitchingFunctionType()) {
      string funcType = simParams_->getSwitchingFunctionType();
      toUpper(funcType);
      if (funcType == "CUBIC") {
        sft_ = cubic;
      } else {
        if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
          sft_ = fifth_order_poly;
        } else {
          // throw error        
          sprintf( painCave.errMsg,
                   "InteractionManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
                   "\tswitchingFunctionType must be one of: "
                   "\"cubic\" or \"fifth_order_polynomial\".", 
                   funcType.c_str() );
          painCave.isFatal = 1;
          painCave.severity = OPENMD_ERROR;
          simError();
        }           
      }
    }

    switcher_->setSwitchType(sft_);
    switcher_->setSwitch(rSwitch_, rCut_);
  }

  /**
   * setupNeighborlists
   *
   *  If the skinThickness was explicitly set, use that value (but check it)
   *  If the skinThickness was not explicitly set: use 1.0 angstroms
   */
  void InteractionManager::setupNeighborlists() {  

    Globals* simParams_ = info_->getSimParams();    
  
    if (simParams_->haveSkinThickness()) {
      skinThickness_ = simParams_->getSkinThickness();
    } else {      
      skinThickness_ = 1.0;
      sprintf(painCave.errMsg,
              "InteractionManager::setupNeighborlists: No value was set for the skinThickness.\n"
              "\tOpenMD will use a default value of %f Angstroms\n"
              "\tfor this simulation\n", skinThickness_);
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();
    }             

    listRadius_ = rCut_ + skinThickness_;
  }


  void InteractionManager::doPrePair(InteractionData idat){
    
    if (!initialized_) initialize();
         
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[idat.atypes].begin(); 
         it != interactions_[idat.atypes].end(); ++it){
      if ((*it)->getFamily() == METALLIC_FAMILY) {
        dynamic_cast<MetallicInteraction*>(*it)->calcDensity(idat);
      }
    }
    
    return;    
  }
  
  void InteractionManager::doPreForce(SelfData sdat){

    if (!initialized_) initialize();
    
    pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
    set<NonBondedInteraction*>::iterator it;
    
    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
      if ((*it)->getFamily() == METALLIC_FAMILY) {
        dynamic_cast<MetallicInteraction*>(*it)->calcFunctional(sdat);
      }
    }
    
    return;    
  }

  void InteractionManager::doPair(InteractionData idat){
    
    if (!initialized_) initialize();
   
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[idat.atypes].begin(); 
         it != interactions_[idat.atypes].end(); ++it)
      (*it)->calcForce(idat);
    
    return;    
  }

  void InteractionManager::doSkipCorrection(InteractionData idat){

    if (!initialized_) initialize();  
    
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[idat.atypes].begin(); 
         it != interactions_[idat.atypes].end(); ++it){
      if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
        dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(idat);
      }
    }
    
    return;    
  }

  void InteractionManager::doSelfCorrection(SelfData sdat){

    if (!initialized_) initialize();
    
    pair<AtomType*, AtomType*> key = make_pair(sdat.atype, sdat.atype);
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
      if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
        dynamic_cast<ElectrostaticInteraction*>(*it)->calcSelfCorrection(sdat);
      }
    }
      
    return;    
  }

  RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
    if (!initialized_) initialize();
    
    AtomType* atype = typeMap_[*atid];

    pair<AtomType*, AtomType*> key = make_pair(atype, atype);
    set<NonBondedInteraction*>::iterator it;
    RealType cutoff = 0.0;
    
    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
      cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));   
    return cutoff;    
  }

  RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
    if (!initialized_) initialize();
    
    pair<AtomType*, AtomType*> key = make_pair(atype, atype);
    set<NonBondedInteraction*>::iterator it;
    RealType cutoff = 0.0;
    
    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
      cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));   
    return cutoff;    
  }
} //end namespace OpenMD
Revision:	1545
Committed:	Fri Apr 8 21:25:19 2011 UTC (14 years ago) by gezelter
File size:	22181 byte(s)
Log Message:	still busted, but much progress
#	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	#include "nonbonded/InteractionManager.hpp"
43	#include "UseTheForce/doForces_interface.h"
44
45	namespace OpenMD {
46
47	InteractionManager* InteractionManager::_instance = NULL;
48	SimInfo* InteractionManager::info_ = NULL;
49	bool InteractionManager::initialized_ = false;
50
51	RealType InteractionManager::rCut_ = 0.0;
52	RealType InteractionManager::rSwitch_ = 0.0;
53	RealType InteractionManager::skinThickness_ = 0.0;
54	RealType InteractionManager::listRadius_ = 0.0;
55	CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE;
56	SwitchingFunctionType InteractionManager::sft_ = cubic;
57	RealType InteractionManager::vdwScale_[4] = {1.0, 0.0, 0.0, 0.0};
58	RealType InteractionManager::electrostaticScale_[4] = {1.0, 0.0, 0.0, 0.0};
59
60	map<int, AtomType*> InteractionManager::typeMap_;
61	map<pair<AtomType, AtomType>, set<NonBondedInteraction*> > InteractionManager::interactions_;
62
63	LJ* InteractionManager::lj_ = new LJ();
64	GB* InteractionManager::gb_ = new GB();
65	Sticky* InteractionManager::sticky_ = new Sticky();
66	Morse* InteractionManager::morse_ = new Morse();
67	EAM* InteractionManager::eam_ = new EAM();
68	SC* InteractionManager::sc_ = new SC();
69	Electrostatic* InteractionManager::electrostatic_ = new Electrostatic();
70	MAW* InteractionManager::maw_ = new MAW();
71	SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
72
73	InteractionManager* InteractionManager::Instance() {
74	if (!_instance) {
75	_instance = new InteractionManager();
76	}
77	return _instance;
78	}
79
80	void InteractionManager::initialize() {
81
82	ForceField* forceField_ = info_->getForceField();
83
84	lj_->setForceField(forceField_);
85	gb_->setForceField(forceField_);
86	sticky_->setForceField(forceField_);
87	eam_->setForceField(forceField_);
88	sc_->setForceField(forceField_);
89	morse_->setForceField(forceField_);
90	electrostatic_->setForceField(forceField_);
91	maw_->setForceField(forceField_);
92
93	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
94
95	// Force fields can set options on how to scale van der Waals and electrostatic
96	// interactions for atoms connected via bonds, bends and torsions
97	// in this case the topological distance between atoms is:
98	// 0 = topologically unconnected
99	// 1 = bonded together
100	// 2 = connected via a bend
101	// 3 = connected via a torsion
102
103	vdwScale_[0] = 1.0;
104	vdwScale_[1] = fopts.getvdw12scale();
105	vdwScale_[2] = fopts.getvdw13scale();
106	vdwScale_[3] = fopts.getvdw14scale();
107
108	electrostaticScale_[0] = 1.0;
109	electrostaticScale_[1] = fopts.getelectrostatic12scale();
110	electrostaticScale_[2] = fopts.getelectrostatic13scale();
111	electrostaticScale_[3] = fopts.getelectrostatic14scale();
112
113	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
114	ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
115	AtomType* atype1;
116	AtomType* atype2;
117	pair<AtomType, AtomType> key;
118	pair<set<NonBondedInteraction*>::iterator, bool> ret;
119
120	for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
121	atype1 = atomTypes->nextType(i1)) {
122
123	// add it to the map:
124	AtomTypeProperties atp = atype1->getATP();
125
126	pair<map<int,AtomType*>::iterator,bool> ret;
127	ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
128	if (ret.second == false) {
129	sprintf( painCave.errMsg,
130	"InteractionManager already had a previous entry with ident %d\n",
131	atp.ident);
132	painCave.severity = OPENMD_INFO;
133	painCave.isFatal = 0;
134	simError();
135	}
136	}
137
138	// Now, iterate over all known types and add to the interaction map:
139
140	map<int, AtomType*>::iterator it1, it2;
141	for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
142	atype1 = (*it1).second;
143
144	for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {
145	atype2 = (*it2).second;
146
147	bool vdwExplicit = false;
148	bool metExplicit = false;
149	bool hbExplicit = false;
150
151	key = make_pair(atype1, atype2);
152
153	if (atype1->isLennardJones() && atype2->isLennardJones()) {
154	interactions_[key].insert(lj_);
155	}
156	if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
157	interactions_[key].insert(electrostatic_);
158	}
159	if (atype1->isSticky() && atype2->isSticky() ) {
160	interactions_[key].insert(sticky_);
161	}
162	if (atype1->isStickyPower() && atype2->isStickyPower() ) {
163	interactions_[key].insert(sticky_);
164	}
165	if (atype1->isEAM() && atype2->isEAM() ) {
166	interactions_[key].insert(eam_);
167	}
168	if (atype1->isSC() && atype2->isSC() ) {
169	interactions_[key].insert(sc_);
170	}
171	if (atype1->isGayBerne() && atype2->isGayBerne() ) {
172	interactions_[key].insert(gb_);
173	}
174	if ((atype1->isGayBerne() && atype2->isLennardJones())
175	\|\| (atype1->isLennardJones() && atype2->isGayBerne())) {
176	interactions_[key].insert(gb_);
177	}
178
179	// look for an explicitly-set non-bonded interaction type using the
180	// two atom types.
181	NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());
182
183	if (nbiType != NULL) {
184
185	if (nbiType->isLennardJones()) {
186	// We found an explicit Lennard-Jones interaction.
187	// override all other vdw entries for this pair of atom types:
188	set<NonBondedInteraction*>::iterator it;
189	for (it = interactions_[key].begin();
190	it != interactions_[key].end(); ++it) {
191	InteractionFamily ifam = (*it)->getFamily();
192	if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
193	}
194	interactions_[key].insert(lj_);
195	vdwExplicit = true;
196	}
197
198	if (nbiType->isMorse()) {
199	if (vdwExplicit) {
200	sprintf( painCave.errMsg,
201	"InteractionManager::initialize found more than one "
202	"explicit \n"
203	"\tvan der Waals interaction for atom types %s - %s\n",
204	atype1->getName().c_str(), atype2->getName().c_str());
205	painCave.severity = OPENMD_ERROR;
206	painCave.isFatal = 1;
207	simError();
208	}
209	// We found an explicit Morse interaction.
210	// override all other vdw entries for this pair of atom types:
211	set<NonBondedInteraction*>::iterator it;
212	for (it = interactions_[key].begin();
213	it != interactions_[key].end(); ++it) {
214	InteractionFamily ifam = (*it)->getFamily();
215	if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
216	}
217	interactions_[key].insert(morse_);
218	vdwExplicit = true;
219	}
220
221	if (nbiType->isEAM()) {
222	// We found an explicit EAM interaction.
223	// override all other metallic entries for this pair of atom types:
224	set<NonBondedInteraction*>::iterator it;
225	for (it = interactions_[key].begin();
226	it != interactions_[key].end(); ++it) {
227	InteractionFamily ifam = (*it)->getFamily();
228	if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
229	}
230	interactions_[key].insert(eam_);
231	metExplicit = true;
232	}
233
234	if (nbiType->isSC()) {
235	if (metExplicit) {
236	sprintf( painCave.errMsg,
237	"InteractionManager::initialize found more than one "
238	"explicit\n"
239	"\tmetallic interaction for atom types %s - %s\n",
240	atype1->getName().c_str(), atype2->getName().c_str());
241	painCave.severity = OPENMD_ERROR;
242	painCave.isFatal = 1;
243	simError();
244	}
245	// We found an explicit Sutton-Chen interaction.
246	// override all other metallic entries for this pair of atom types:
247	set<NonBondedInteraction*>::iterator it;
248	for (it = interactions_[key].begin();
249	it != interactions_[key].end(); ++it) {
250	InteractionFamily ifam = (*it)->getFamily();
251	if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
252	}
253	interactions_[key].insert(sc_);
254	metExplicit = true;
255	}
256
257	if (nbiType->isMAW()) {
258	if (vdwExplicit) {
259	sprintf( painCave.errMsg,
260	"InteractionManager::initialize found more than one "
261	"explicit\n"
262	"\tvan der Waals interaction for atom types %s - %s\n",
263	atype1->getName().c_str(), atype2->getName().c_str());
264	painCave.severity = OPENMD_ERROR;
265	painCave.isFatal = 1;
266	simError();
267	}
268	// We found an explicit MAW interaction.
269	// override all other vdw entries for this pair of atom types:
270	set<NonBondedInteraction*>::iterator it;
271	for (it = interactions_[key].begin();
272	it != interactions_[key].end(); ++it) {
273	InteractionFamily ifam = (*it)->getFamily();
274	if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
275	}
276	interactions_[key].insert(maw_);
277	vdwExplicit = true;
278	}
279	}
280	}
281	}
282
283
284	// make sure every pair of atom types in this simulation has a
285	// non-bonded interaction:
286
287	set<AtomType*> simTypes = info_->getSimulatedAtomTypes();
288	set<AtomType*>::iterator it, jt;
289	for (it = simTypes.begin(); it != simTypes.end(); ++it) {
290	atype1 = (*it);
291	for (jt = simTypes.begin(); jt != simTypes.end(); ++jt) {
292	atype2 = (*jt);
293	key = make_pair(atype1, atype2);
294
295	if (interactions_[key].size() == 0) {
296	sprintf( painCave.errMsg,
297	"InteractionManager unable to find an appropriate non-bonded\n"
298	"\tinteraction for atom types %s - %s\n",
299	atype1->getName().c_str(), atype2->getName().c_str());
300	painCave.severity = OPENMD_INFO;
301	painCave.isFatal = 1;
302	simError();
303	}
304	}
305	}
306
307	setupCutoffs();
308	setupSwitching();
309	setupNeighborlists();
310
311	//int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0;
312	//int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0;
313	//notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf);
314	//notifyFortranSkinThickness(&skinThickness_);
315
316	initialized_ = true;
317	}
318
319	/**
320	* setupCutoffs
321	*
322	* Sets the values of cutoffRadius and cutoffMethod
323	*
324	* cutoffRadius : realType
325	* If the cutoffRadius was explicitly set, use that value.
326	* If the cutoffRadius was not explicitly set:
327	* Are there electrostatic atoms? Use 12.0 Angstroms.
328	* No electrostatic atoms? Poll the atom types present in the
329	* simulation for suggested cutoff values (e.g. 2.5 * sigma).
330	* Use the maximum suggested value that was found.
331	*
332	* cutoffMethod : (one of HARD, SWITCHED, SHIFTED_FORCE, SHIFTED_POTENTIAL)
333	* If cutoffMethod was explicitly set, use that choice.
334	* If cutoffMethod was not explicitly set, use SHIFTED_FORCE
335	*/
336	void InteractionManager::setupCutoffs() {
337
338	Globals* simParams_ = info_->getSimParams();
339
340	if (simParams_->haveCutoffRadius()) {
341	rCut_ = simParams_->getCutoffRadius();
342	} else {
343	if (info_->usesElectrostaticAtoms()) {
344	sprintf(painCave.errMsg,
345	"InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n"
346	"\tOpenMD will use a default value of 12.0 angstroms"
347	"\tfor the cutoffRadius.\n");
348	painCave.isFatal = 0;
349	painCave.severity = OPENMD_INFO;
350	simError();
351	rCut_ = 12.0;
352	} else {
353	RealType thisCut;
354	set<AtomType*>::iterator i;
355	set<AtomType*> atomTypes;
356	atomTypes = info_->getSimulatedAtomTypes();
357	for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
358	thisCut = getSuggestedCutoffRadius((*i));
359	rCut_ = max(thisCut, rCut_);
360	}
361	sprintf(painCave.errMsg,
362	"InteractionManager::setupCutoffs: No value was set for the cutoffRadius.\n"
363	"\tOpenMD will use %lf angstroms.\n",
364	rCut_);
365	painCave.isFatal = 0;
366	painCave.severity = OPENMD_INFO;
367	simError();
368	}
369	}
370
371	map<string, CutoffMethod> stringToCutoffMethod;
372	stringToCutoffMethod["HARD"] = HARD;
373	stringToCutoffMethod["SWITCHED"] = SWITCHED;
374	stringToCutoffMethod["SHIFTED_POTENTIAL"] = SHIFTED_POTENTIAL;
375	stringToCutoffMethod["SHIFTED_FORCE"] = SHIFTED_FORCE;
376
377	if (simParams_->haveCutoffMethod()) {
378	string cutMeth = toUpperCopy(simParams_->getCutoffMethod());
379	map<string, CutoffMethod>::iterator i;
380	i = stringToCutoffMethod.find(cutMeth);
381	if (i == stringToCutoffMethod.end()) {
382	sprintf(painCave.errMsg,
383	"InteractionManager::setupCutoffs: Could not find chosen cutoffMethod %s\n"
384	"\tShould be one of: "
385	"HARD, SWITCHED, SHIFTED_POTENTIAL, or SHIFTED_FORCE\n",
386	cutMeth.c_str());
387	painCave.isFatal = 1;
388	painCave.severity = OPENMD_ERROR;
389	simError();
390	} else {
391	cutoffMethod_ = i->second;
392	}
393	} else {
394	sprintf(painCave.errMsg,
395	"InteractionManager::setupCutoffs: No value was set for the cutoffMethod.\n"
396	"\tOpenMD will use SHIFTED_FORCE.\n");
397	painCave.isFatal = 0;
398	painCave.severity = OPENMD_INFO;
399	simError();
400	cutoffMethod_ = SHIFTED_FORCE;
401	}
402	}
403
404
405	/**
406	* setupSwitching
407	*
408	* Sets the values of switchingRadius and
409	* If the switchingRadius was explicitly set, use that value (but check it)
410	* If the switchingRadius was not explicitly set: use 0.85 * cutoffRadius_
411	*/
412	void InteractionManager::setupSwitching() {
413	Globals* simParams_ = info_->getSimParams();
414
415	if (simParams_->haveSwitchingRadius()) {
416	rSwitch_ = simParams_->getSwitchingRadius();
417	if (rSwitch_ > rCut_) {
418	sprintf(painCave.errMsg,
419	"InteractionManager::setupSwitching: switchingRadius (%f) is larger than cutoffRadius(%f)\n",
420	rSwitch_, rCut_);
421	painCave.isFatal = 1;
422	painCave.severity = OPENMD_ERROR;
423	simError();
424	}
425	} else {
426	rSwitch_ = 0.85 * rCut_;
427	sprintf(painCave.errMsg,
428	"InteractionManager::setupSwitching: No value was set for the switchingRadius.\n"
429	"\tOpenMD will use a default value of 85 percent of the cutoffRadius.\n"
430	"\tswitchingRadius = %f. for this simulation\n", rSwitch_);
431	painCave.isFatal = 0;
432	painCave.severity = OPENMD_WARNING;
433	simError();
434	}
435
436	if (simParams_->haveSwitchingFunctionType()) {
437	string funcType = simParams_->getSwitchingFunctionType();
438	toUpper(funcType);
439	if (funcType == "CUBIC") {
440	sft_ = cubic;
441	} else {
442	if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
443	sft_ = fifth_order_poly;
444	} else {
445	// throw error
446	sprintf( painCave.errMsg,
447	"InteractionManager::setupSwitching : Unknown switchingFunctionType. (Input file specified %s .)\n"
448	"\tswitchingFunctionType must be one of: "
449	"\"cubic\" or \"fifth_order_polynomial\".",
450	funcType.c_str() );
451	painCave.isFatal = 1;
452	painCave.severity = OPENMD_ERROR;
453	simError();
454	}
455	}
456	}
457
458	switcher_->setSwitchType(sft_);
459	switcher_->setSwitch(rSwitch_, rCut_);
460	}
461
462	/**
463	* setupNeighborlists
464	*
465	* If the skinThickness was explicitly set, use that value (but check it)
466	* If the skinThickness was not explicitly set: use 1.0 angstroms
467	*/
468	void InteractionManager::setupNeighborlists() {
469
470	Globals* simParams_ = info_->getSimParams();
471
472	if (simParams_->haveSkinThickness()) {
473	skinThickness_ = simParams_->getSkinThickness();
474	} else {
475	skinThickness_ = 1.0;
476	sprintf(painCave.errMsg,
477	"InteractionManager::setupNeighborlists: No value was set for the skinThickness.\n"
478	"\tOpenMD will use a default value of %f Angstroms\n"
479	"\tfor this simulation\n", skinThickness_);
480	painCave.severity = OPENMD_INFO;
481	painCave.isFatal = 0;
482	simError();
483	}
484
485	listRadius_ = rCut_ + skinThickness_;
486	}
487
488
489	void InteractionManager::doPrePair(InteractionData idat){
490
491	if (!initialized_) initialize();
492
493	set<NonBondedInteraction*>::iterator it;
494
495	for (it = interactions_[idat.atypes].begin();
496	it != interactions_[idat.atypes].end(); ++it){
497	if ((*it)->getFamily() == METALLIC_FAMILY) {
498	dynamic_cast<MetallicInteraction>(it)->calcDensity(idat);
499	}
500	}
501
502	return;
503	}
504
505	void InteractionManager::doPreForce(SelfData sdat){
506
507	if (!initialized_) initialize();
508
509	pair<AtomType, AtomType> key = make_pair(sdat.atype, sdat.atype);
510	set<NonBondedInteraction*>::iterator it;
511
512	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
513	if ((*it)->getFamily() == METALLIC_FAMILY) {
514	dynamic_cast<MetallicInteraction>(it)->calcFunctional(sdat);
515	}
516	}
517
518	return;
519	}
520
521	void InteractionManager::doPair(InteractionData idat){
522
523	if (!initialized_) initialize();
524
525	set<NonBondedInteraction*>::iterator it;
526
527	for (it = interactions_[idat.atypes].begin();
528	it != interactions_[idat.atypes].end(); ++it)
529	(*it)->calcForce(idat);
530
531	return;
532	}
533
534	void InteractionManager::doSkipCorrection(InteractionData idat){
535
536	if (!initialized_) initialize();
537
538	set<NonBondedInteraction*>::iterator it;
539
540	for (it = interactions_[idat.atypes].begin();
541	it != interactions_[idat.atypes].end(); ++it){
542	if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
543	dynamic_cast<ElectrostaticInteraction>(it)->calcSkipCorrection(idat);
544	}
545	}
546
547	return;
548	}
549
550	void InteractionManager::doSelfCorrection(SelfData sdat){
551
552	if (!initialized_) initialize();
553
554	pair<AtomType, AtomType> key = make_pair(sdat.atype, sdat.atype);
555	set<NonBondedInteraction*>::iterator it;
556
557	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
558	if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
559	dynamic_cast<ElectrostaticInteraction>(it)->calcSelfCorrection(sdat);
560	}
561	}
562
563	return;
564	}
565
566	RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
567	if (!initialized_) initialize();
568
569	AtomType* atype = typeMap_[*atid];
570
571	pair<AtomType, AtomType> key = make_pair(atype, atype);
572	set<NonBondedInteraction*>::iterator it;
573	RealType cutoff = 0.0;
574
575	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
576	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
577	return cutoff;
578	}
579
580	RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
581	if (!initialized_) initialize();
582
583	pair<AtomType, AtomType> key = make_pair(atype, atype);
584	set<NonBondedInteraction*>::iterator it;
585	RealType cutoff = 0.0;
586
587	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
588	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
589	return cutoff;
590	}
591	} //end namespace OpenMD