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"

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