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.
#	User	Rev	Content
1	gezelter	1502	/*
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	gezelter	1535
46			InteractionManager* InteractionManager::_instance = NULL;
47			SimInfo* InteractionManager::info_ = NULL;
48	gezelter	1502	bool InteractionManager::initialized_ = false;
49	gezelter	1535
50			RealType InteractionManager::rCut_ = 0.0;
51			RealType InteractionManager::rSwitch_ = 0.0;
52			CutoffMethod InteractionManager::cutoffMethod_ = SHIFTED_FORCE;
53			SwitchingFunctionType InteractionManager::sft_ = cubic;
54	gezelter	1536	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	gezelter	1535
57	gezelter	1502	map<int, AtomType*> InteractionManager::typeMap_;
58			map<pair<AtomType, AtomType>, set<NonBondedInteraction*> > InteractionManager::interactions_;
59	gezelter	1528
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	gezelter	1532	MAW* InteractionManager::maw_ = new MAW();
68	gezelter	1530	SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
69	gezelter	1502
70			InteractionManager* InteractionManager::Instance() {
71			if (!_instance) {
72			_instance = new InteractionManager();
73			}
74			return _instance;
75			}
76
77			void InteractionManager::initialize() {
78
79	gezelter	1535	ForceField* forceField_ = info_->getForceField();
80
81	gezelter	1502	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	gezelter	1532	maw_->setForceField(forceField_);
89	gezelter	1502
90	gezelter	1530	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
91	gezelter	1535
92	gezelter	1530	// 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	gezelter	1536	// 0 = topologically unconnected
96	gezelter	1530	// 1 = bonded together
97			// 2 = connected via a bend
98			// 3 = connected via a torsion
99
100	gezelter	1536	vdwScale_[0] = 1.0;
101	gezelter	1530	vdwScale_[1] = fopts.getvdw12scale();
102			vdwScale_[2] = fopts.getvdw13scale();
103			vdwScale_[3] = fopts.getvdw14scale();
104
105	gezelter	1536	electrostaticScale_[0] = 1.0;
106	gezelter	1530	electrostaticScale_[1] = fopts.getelectrostatic12scale();
107			electrostaticScale_[2] = fopts.getelectrostatic13scale();
108	gezelter	1535	electrostaticScale_[3] = fopts.getelectrostatic14scale();
109	gezelter	1530
110	gezelter	1502	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	gezelter	1535
180			if (nbiType != NULL) {
181	gezelter	1502
182	gezelter	1535	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	gezelter	1502	}
194	gezelter	1535
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	gezelter	1502	}
217	gezelter	1535
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	gezelter	1502	}
230	gezelter	1535
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	gezelter	1502	}
253	gezelter	1535
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	gezelter	1502	}
277			}
278			}
279
280	gezelter	1535
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	gezelter	1502	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	gezelter	1535
304			setupCutoffs();
305			setupSwitching();
306
307	gezelter	1545	//int ljsp = cutoffMethod_ == SHIFTED_POTENTIAL ? 1 : 0;
308			//int ljsf = cutoffMethod_ == SHIFTED_FORCE ? 1 : 0;
309			//notifyFortranCutoffs(&rCut_, &rSwitch_, &ljsp, &ljsf);
310	gezelter	1535
311			initialized_ = true;
312			}
313	gezelter	1502
314	gezelter	1535	/**
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	gezelter	1536
453			switcher_->setSwitchType(sft_);
454			switcher_->setSwitch(rSwitch_, rCut_);
455	gezelter	1535	}
456
457	gezelter	1545	void InteractionManager::doPrePair(InteractionData idat){
458	gezelter	1502
459			if (!initialized_) initialize();
460	gezelter	1545
461	gezelter	1504	set<NonBondedInteraction*>::iterator it;
462
463	gezelter	1554	for (it = interactions_[ *(idat.atypes) ].begin();
464			it != interactions_[ *(idat.atypes) ].end(); ++it){
465	gezelter	1504	if ((*it)->getFamily() == METALLIC_FAMILY) {
466	gezelter	1545	dynamic_cast<MetallicInteraction>(it)->calcDensity(idat);
467	gezelter	1504	}
468			}
469	gezelter	1502
470			return;
471			}
472	gezelter	1504
473	gezelter	1545	void InteractionManager::doPreForce(SelfData sdat){
474	gezelter	1502
475	gezelter	1504	if (!initialized_) initialize();
476	gezelter	1545
477			pair<AtomType, AtomType> key = make_pair(sdat.atype, sdat.atype);
478	gezelter	1504	set<NonBondedInteraction*>::iterator it;
479	gezelter	1502
480	gezelter	1504	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
481			if ((*it)->getFamily() == METALLIC_FAMILY) {
482	gezelter	1545	dynamic_cast<MetallicInteraction>(it)->calcFunctional(sdat);
483	gezelter	1504	}
484			}
485
486	gezelter	1502	return;
487			}
488
489	gezelter	1545	void InteractionManager::doPair(InteractionData idat){
490	gezelter	1504
491	gezelter	1502	if (!initialized_) initialize();
492	gezelter	1545
493	gezelter	1502	set<NonBondedInteraction*>::iterator it;
494
495	gezelter	1554	for (it = interactions_[ *(idat.atypes) ].begin();
496			it != interactions_[ *(idat.atypes) ].end(); ++it)
497	gezelter	1502	(*it)->calcForce(idat);
498
499			return;
500			}
501
502	gezelter	1545	void InteractionManager::doSkipCorrection(InteractionData idat){
503	gezelter	1502
504	gezelter	1545	if (!initialized_) initialize();
505	gezelter	1504
506			set<NonBondedInteraction*>::iterator it;
507
508	gezelter	1554	for (it = interactions_[ *(idat.atypes) ].begin();
509			it != interactions_[ *(idat.atypes) ].end(); ++it){
510	gezelter	1504	if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
511	gezelter	1545	dynamic_cast<ElectrostaticInteraction>(it)->calcSkipCorrection(idat);
512	gezelter	1504	}
513			}
514	gezelter	1502
515	gezelter	1504	return;
516	gezelter	1502	}
517
518	gezelter	1545	void InteractionManager::doSelfCorrection(SelfData sdat){
519	gezelter	1502
520			if (!initialized_) initialize();
521	gezelter	1504
522	gezelter	1545	pair<AtomType, AtomType> key = make_pair(sdat.atype, sdat.atype);
523	gezelter	1504	set<NonBondedInteraction*>::iterator it;
524	gezelter	1502
525	gezelter	1504	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
526			if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
527	gezelter	1545	dynamic_cast<ElectrostaticInteraction>(it)->calcSelfCorrection(sdat);
528	gezelter	1504	}
529			}
530	gezelter	1536
531	gezelter	1504	return;
532	gezelter	1502	}
533
534	gezelter	1505	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	gezelter	1545	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
545	gezelter	1505	return cutoff;
546			}
547
548	gezelter	1528	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	gezelter	1545	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(key));
557	gezelter	1528	return cutoff;
558			}
559	gezelter	1502	} //end namespace OpenMD