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

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

    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_);
  }

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