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 {
 
  bool InteractionManager::initialized_ = false;
  RealType InteractionManager::rCut_ = -1.0;
  RealType InteractionManager::rSwitch_ = -1.0;
  ForceField* InteractionManager::forceField_ = NULL;  
  InteractionManager* InteractionManager::_instance = NULL;
  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();
  SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
 
  InteractionManager* InteractionManager::Instance() {
    if (!_instance) {
      _instance = new InteractionManager();
    }
    return _instance;
  }

  void InteractionManager::initialize() {
    
    lj_->setForceField(forceField_);
    gb_->setForceField(forceField_);
    sticky_->setForceField(forceField_);
    eam_->setForceField(forceField_);
    sc_->setForceField(forceField_);
    morse_->setForceField(forceField_);
    electrostatic_->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 = the atom itself
    // 1 = bonded together 
    // 2 = connected via a bend
    // 3 = connected via a torsion

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

    electrostaticScale_[0] = 0.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->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;
        }
      }
    }
    
    // make sure every pair of atom types has a non-bonded interaction:
    for (atype1 = atomTypes->beginType(i1); atype1 != NULL; 
         atype1 = atomTypes->nextType(i1)) {
      for (atype2 = atomTypes->beginType(i2); atype2 != NULL; 
           atype2 = atomTypes->nextType(i2)) {
        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();
        }
      }
    }
  } 
  
  void InteractionManager::doPrePair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){
    
    if (!initialized_) initialize();
          
    DensityData ddat;

    ddat.atype1 = typeMap_[*atid1];
    ddat.atype2 = typeMap_[*atid2];
    ddat.rij = *rij;
    ddat.rho_i_at_j = *rho_i_at_j;
    ddat.rho_j_at_i = *rho_j_at_i;

    pair<AtomType*, AtomType*> key = make_pair(ddat.atype1, ddat.atype2);
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
      if ((*it)->getFamily() == METALLIC_FAMILY) {
        dynamic_cast<MetallicInteraction*>(*it)->calcDensity(ddat);
      }
    }
    
    return;    
  }
  
  void InteractionManager::doPreForce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){

    if (!initialized_) initialize();
          
    FunctionalData fdat;

    fdat.atype = typeMap_[*atid];
    fdat.rho = *rho;
    fdat.frho = *frho;
    fdat.dfrhodrho = *dfrhodrho;

    pair<AtomType*, AtomType*> key = make_pair(fdat.atype, fdat.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(fdat);
      }
    }
    
    return;    
  }

  void InteractionManager::doPair(int *atid1, int *atid2, RealType *d, RealType *r, RealType *r2, RealType *rcut, RealType *sw, int *topoDist, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *A1, RealType *A2, RealType *t1, RealType *t2, RealType *rho1, RealType *rho2, RealType *dfrho1, RealType *dfrho2, RealType *fshift1, RealType *fshift2){
    
    if (!initialized_) initialize();
    
    InteractionData idat;
    
    idat.atype1 = typeMap_[*atid1];
    idat.atype2 = typeMap_[*atid2];
    idat.d = Vector3d(d);
    idat.rij = *r;
    idat.r2 = *r2;
    idat.rcut = *rcut;
    idat.sw = *sw;
    idat.vdwMult = vdwScale_[*topoDist];
    idat.electroMult = electrostaticScale_[*topoDist];
    idat.pot = *pot;
    idat.vpair = *vpair;
    idat.f1 = Vector3d(f1);
    idat.eFrame1 = Mat3x3d(eFrame1);
    idat.eFrame2 = Mat3x3d(eFrame2);
    idat.A1 = RotMat3x3d(A1);
    idat.A2 = RotMat3x3d(A2);
    idat.t1 = Vector3d(t1);
    idat.t2 = Vector3d(t2);
    idat.rho1 = *rho1;
    idat.rho2 = *rho2;
    idat.dfrho1 = *dfrho1;
    idat.dfrho2 = *dfrho2;
    idat.fshift1 = *fshift1;
    idat.fshift2 = *fshift2;

    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
      (*it)->calcForce(idat);
    
    f1[0] = idat.f1.x();
    f1[1] = idat.f1.y();
    f1[2] = idat.f1.z();
    
    t1[0] = idat.t1.x();
    t1[1] = idat.t1.y();
    t1[2] = idat.t1.z();
    
    t2[0] = idat.t2.x();
    t2[1] = idat.t2.y();
    t2[2] = idat.t2.z();

    return;    
  }

  void InteractionManager::doSkipCorrection(int *atid1, int *atid2, RealType *d, RealType *r, RealType *skippedCharge1, RealType *skippedCharge2, RealType *sw, RealType *electroMult, RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, RealType *eFrame2, RealType *t1, RealType *t2){

    if (!initialized_) initialize();
    
    SkipCorrectionData skdat;
    
    skdat.atype1 = typeMap_[*atid1];
    skdat.atype2 = typeMap_[*atid2];
    skdat.d = Vector3d(d);
    skdat.rij = *r;
    skdat.skippedCharge1 = *skippedCharge1;
    skdat.skippedCharge2 = *skippedCharge2;
    skdat.sw = *sw;
    skdat.electroMult = *electroMult;
    skdat.pot = *pot;
    skdat.vpair = *vpair;
    skdat.f1 = Vector3d(f1);
    skdat.eFrame1 = Mat3x3d(eFrame1);
    skdat.eFrame2 = Mat3x3d(eFrame2);
    skdat.t1 = Vector3d(t1);
    skdat.t2 = Vector3d(t2);

    pair<AtomType*, AtomType*> key = make_pair(skdat.atype1, skdat.atype2);
    set<NonBondedInteraction*>::iterator it;

    for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
      if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
        dynamic_cast<ElectrostaticInteraction*>(*it)->calcSkipCorrection(skdat);
      }
    }
    
    f1[0] = skdat.f1.x();
    f1[1] = skdat.f1.y();
    f1[2] = skdat.f1.z();
    
    t1[0] = skdat.t1.x();
    t1[1] = skdat.t1.y();
    t1[2] = skdat.t1.z();
    
    t2[0] = skdat.t2.x();
    t2[1] = skdat.t2.y();
    t2[2] = skdat.t2.z();

    return;    
  }

  void InteractionManager::doSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){

    if (!initialized_) initialize();
    
    SelfCorrectionData scdat;
    
    scdat.atype = typeMap_[*atid];
    scdat.eFrame = Mat3x3d(eFrame);
    scdat.skippedCharge = *skippedCharge;
    scdat.pot = *pot;
    scdat.t = Vector3d(t);

    pair<AtomType*, AtomType*> key = make_pair(scdat.atype, scdat.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(scdat);
      }
    }
        
    t[0] = scdat.t.x();
    t[1] = scdat.t.y();
    t[2] = scdat.t.z();

    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(atype, atype));   
    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(atype, atype));   
    return cutoff;    
  }


  void InteractionManager::setSwitch(RealType *rIn, RealType *rOut) {
    switcher_->setSwitch(*rIn, *rOut);    
  }

  void InteractionManager::getSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r,
                                     int *in_switching_region) {
    bool isr = switcher_->getSwitch(*r2, *sw, *dswdr, *r);    
    *in_switching_region = (int)isr;
  }

} //end namespace OpenMD

extern "C" {
  
#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR)
#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE)
#define fortranDoPair FC_FUNC(do_pair, DO_PAIR)
#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION)
#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION)
#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF)
#define fortranSetSwitch FC_FUNC(set_switch, SET_SWITCH)
#define fortranGetSwitch FC_FUNC(get_switch, GET_SWITCH)

  void fortranDoPrePair(int *atid1, int *atid2, RealType *rij, 
                        RealType *rho_i_at_j, RealType *rho_j_at_i) {
            
    return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij,
                                                             rho_i_at_j,  
                                                             rho_j_at_i);
  }
  void fortranDoPreForce(int *atid, RealType *rho, RealType *frho,
                         RealType *dfrhodrho) {  
    
    return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho,
                                                              dfrhodrho);    
  }
  
  void fortranDoPair(int *atid1, int *atid2, RealType *d, RealType *r, 
                     RealType *r2, RealType *rcut, RealType *sw, int *topoDist,
                     RealType *pot, RealType *vpair, RealType *f1, RealType *eFrame1, 
                     RealType *eFrame2, RealType *A1, RealType *A2, 
                     RealType *t1, RealType *t2, RealType *rho1, RealType *rho2,
                     RealType *dfrho1, RealType *dfrho2, RealType *fshift1, 
                     RealType *fshift2){
    
    return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r, 
                                                          r2, rcut, sw, topoDist,
                                                          pot, vpair, f1, 
                                                          eFrame1, eFrame2, 
                                                          A1, A2, t1, t2, rho1,
                                                          rho2, dfrho1, dfrho2,
                                                          fshift1, fshift2);
  }
  
  void fortranDoSkipCorrection(int *atid1, int *atid2, RealType *d, 
                               RealType *r, RealType *skippedCharge1, 
                               RealType *skippedCharge2, RealType *sw, 
                               RealType *electroMult, RealType *pot, 
                               RealType *vpair, RealType *f1,
                               RealType *eFrame1, RealType *eFrame2,
                               RealType *t1, RealType *t2){
    
    return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1, 
                                                                    atid2, d, 
                                                                    r, 
                                                                    skippedCharge1,
                                                                    skippedCharge2,
                                                                    sw, electroMult, pot,
                                                                    vpair, f1, eFrame1, 
                                                                    eFrame2, t1, t2);
  }
  
  void fortranDoSelfCorrection(int *atid, RealType *eFrame, RealType *skippedCharge, 
                               RealType *pot, RealType *t) {
    
    return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid, 
                                                                    eFrame, 
                                                                    skippedCharge,
                                                                    pot, t);
  }
  RealType fortranGetCutoff(int *atid) {    
    return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid);
  }

  void fortranGetSwitch(RealType *r2, RealType *sw, RealType *dswdr, RealType *r,
                        int *in_switching_region) {
    
    return OpenMD::InteractionManager::Instance()->getSwitch(r2, sw, dswdr, r, 
                                                             in_switching_region);
  }

  void fortranSetSwitch(RealType *rIn, RealType *rOut) {    
    return OpenMD::InteractionManager::Instance()->setSwitch(rIn, rOut);
  }
  
}
Revision:	1530
Committed:	Tue Dec 28 21:47:55 2010 UTC (14 years, 4 months ago) by gezelter
File size:	21779 byte(s)
Log Message:	Moved switching functions and force options over to the C++ side, and removed them from Fortran.
#	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	bool InteractionManager::initialized_ = false;
47	RealType InteractionManager::rCut_ = -1.0;
48	RealType InteractionManager::rSwitch_ = -1.0;
49	ForceField* InteractionManager::forceField_ = NULL;
50	InteractionManager* InteractionManager::_instance = NULL;
51	map<int, AtomType*> InteractionManager::typeMap_;
52	map<pair<AtomType, AtomType>, set<NonBondedInteraction*> > InteractionManager::interactions_;
53
54	LJ* InteractionManager::lj_ = new LJ();
55	GB* InteractionManager::gb_ = new GB();
56	Sticky* InteractionManager::sticky_ = new Sticky();
57	Morse* InteractionManager::morse_ = new Morse();
58	EAM* InteractionManager::eam_ = new EAM();
59	SC* InteractionManager::sc_ = new SC();
60	Electrostatic* InteractionManager::electrostatic_ = new Electrostatic();
61	SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
62
63	InteractionManager* InteractionManager::Instance() {
64	if (!_instance) {
65	_instance = new InteractionManager();
66	}
67	return _instance;
68	}
69
70	void InteractionManager::initialize() {
71
72	lj_->setForceField(forceField_);
73	gb_->setForceField(forceField_);
74	sticky_->setForceField(forceField_);
75	eam_->setForceField(forceField_);
76	sc_->setForceField(forceField_);
77	morse_->setForceField(forceField_);
78	electrostatic_->setForceField(forceField_);
79
80	ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
81	// Force fields can set options on how to scale van der Waals and electrostatic
82	// interactions for atoms connected via bonds, bends and torsions
83	// in this case the topological distance between atoms is:
84	// 0 = the atom itself
85	// 1 = bonded together
86	// 2 = connected via a bend
87	// 3 = connected via a torsion
88
89	vdwScale_[0] = 0.0;
90	vdwScale_[1] = fopts.getvdw12scale();
91	vdwScale_[2] = fopts.getvdw13scale();
92	vdwScale_[3] = fopts.getvdw14scale();
93
94	electrostaticScale_[0] = 0.0;
95	electrostaticScale_[1] = fopts.getelectrostatic12scale();
96	electrostaticScale_[2] = fopts.getelectrostatic13scale();
97	electrostaticScale_[3] = fopts.getelectrostatic14scale();
98
99	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
100	ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
101	AtomType* atype1;
102	AtomType* atype2;
103	pair<AtomType, AtomType> key;
104	pair<set<NonBondedInteraction*>::iterator, bool> ret;
105
106	for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
107	atype1 = atomTypes->nextType(i1)) {
108
109	// add it to the map:
110	AtomTypeProperties atp = atype1->getATP();
111
112	pair<map<int,AtomType*>::iterator,bool> ret;
113	ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
114	if (ret.second == false) {
115	sprintf( painCave.errMsg,
116	"InteractionManager already had a previous entry with ident %d\n",
117	atp.ident);
118	painCave.severity = OPENMD_INFO;
119	painCave.isFatal = 0;
120	simError();
121	}
122	}
123
124	// Now, iterate over all known types and add to the interaction map:
125
126	map<int, AtomType*>::iterator it1, it2;
127	for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
128	atype1 = (*it1).second;
129
130	for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {
131	atype2 = (*it2).second;
132
133	bool vdwExplicit = false;
134	bool metExplicit = false;
135	bool hbExplicit = false;
136
137	key = make_pair(atype1, atype2);
138
139	if (atype1->isLennardJones() && atype2->isLennardJones()) {
140	interactions_[key].insert(lj_);
141	}
142	if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
143	interactions_[key].insert(electrostatic_);
144	}
145	if (atype1->isSticky() && atype2->isSticky() ) {
146	interactions_[key].insert(sticky_);
147	}
148	if (atype1->isStickyPower() && atype2->isStickyPower() ) {
149	interactions_[key].insert(sticky_);
150	}
151	if (atype1->isEAM() && atype2->isEAM() ) {
152	interactions_[key].insert(eam_);
153	}
154	if (atype1->isSC() && atype2->isSC() ) {
155	interactions_[key].insert(sc_);
156	}
157	if (atype1->isGayBerne() && atype2->isGayBerne() ) {
158	interactions_[key].insert(gb_);
159	}
160	if ((atype1->isGayBerne() && atype2->isLennardJones())
161	\|\| (atype1->isLennardJones() && atype2->isGayBerne())) {
162	interactions_[key].insert(gb_);
163	}
164
165	// look for an explicitly-set non-bonded interaction type using the
166	// two atom types.
167	NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());
168
169	if (nbiType->isLennardJones()) {
170	// We found an explicit Lennard-Jones interaction.
171	// override all other vdw entries for this pair of atom types:
172	set<NonBondedInteraction*>::iterator it;
173	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
174	InteractionFamily ifam = (*it)->getFamily();
175	if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
176	}
177	interactions_[key].insert(lj_);
178	vdwExplicit = true;
179	}
180
181	if (nbiType->isMorse()) {
182	if (vdwExplicit) {
183	sprintf( painCave.errMsg,
184	"InteractionManager::initialize found more than one explicit\n"
185	"\tvan der Waals interaction for atom types %s - %s\n",
186	atype1->getName().c_str(), atype2->getName().c_str());
187	painCave.severity = OPENMD_ERROR;
188	painCave.isFatal = 1;
189	simError();
190	}
191	// We found an explicit Morse interaction.
192	// override all other vdw entries for this pair of atom types:
193	set<NonBondedInteraction*>::iterator it;
194	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
195	InteractionFamily ifam = (*it)->getFamily();
196	if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
197	}
198	interactions_[key].insert(morse_);
199	vdwExplicit = true;
200	}
201
202	if (nbiType->isEAM()) {
203	// We found an explicit EAM interaction.
204	// override all other metallic entries for this pair of atom types:
205	set<NonBondedInteraction*>::iterator it;
206	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
207	InteractionFamily ifam = (*it)->getFamily();
208	if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
209	}
210	interactions_[key].insert(eam_);
211	metExplicit = true;
212	}
213
214	if (nbiType->isSC()) {
215	if (metExplicit) {
216	sprintf( painCave.errMsg,
217	"InteractionManager::initialize found more than one explicit\n"
218	"\tmetallic interaction for atom types %s - %s\n",
219	atype1->getName().c_str(), atype2->getName().c_str());
220	painCave.severity = OPENMD_ERROR;
221	painCave.isFatal = 1;
222	simError();
223	}
224	// We found an explicit Sutton-Chen interaction.
225	// override all other metallic entries for this pair of atom types:
226	set<NonBondedInteraction*>::iterator it;
227	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
228	InteractionFamily ifam = (*it)->getFamily();
229	if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
230	}
231	interactions_[key].insert(sc_);
232	metExplicit = true;
233	}
234	}
235	}
236
237	// make sure every pair of atom types has a non-bonded interaction:
238	for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
239	atype1 = atomTypes->nextType(i1)) {
240	for (atype2 = atomTypes->beginType(i2); atype2 != NULL;
241	atype2 = atomTypes->nextType(i2)) {
242	key = make_pair(atype1, atype2);
243
244	if (interactions_[key].size() == 0) {
245	sprintf( painCave.errMsg,
246	"InteractionManager unable to find an appropriate non-bonded\n"
247	"\tinteraction for atom types %s - %s\n",
248	atype1->getName().c_str(), atype2->getName().c_str());
249	painCave.severity = OPENMD_INFO;
250	painCave.isFatal = 1;
251	simError();
252	}
253	}
254	}
255	}
256
257	void InteractionManager::doPrePair(int atid1, int atid2, RealType rij, RealType rho_i_at_j, RealType *rho_j_at_i){
258
259	if (!initialized_) initialize();
260
261	DensityData ddat;
262
263	ddat.atype1 = typeMap_[*atid1];
264	ddat.atype2 = typeMap_[*atid2];
265	ddat.rij = *rij;
266	ddat.rho_i_at_j = *rho_i_at_j;
267	ddat.rho_j_at_i = *rho_j_at_i;
268
269	pair<AtomType, AtomType> key = make_pair(ddat.atype1, ddat.atype2);
270	set<NonBondedInteraction*>::iterator it;
271
272	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
273	if ((*it)->getFamily() == METALLIC_FAMILY) {
274	dynamic_cast<MetallicInteraction>(it)->calcDensity(ddat);
275	}
276	}
277
278	return;
279	}
280
281	void InteractionManager::doPreForce(int atid, RealType rho, RealType frho, RealType dfrhodrho){
282
283	if (!initialized_) initialize();
284
285	FunctionalData fdat;
286
287	fdat.atype = typeMap_[*atid];
288	fdat.rho = *rho;
289	fdat.frho = *frho;
290	fdat.dfrhodrho = *dfrhodrho;
291
292	pair<AtomType, AtomType> key = make_pair(fdat.atype, fdat.atype);
293	set<NonBondedInteraction*>::iterator it;
294
295	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
296	if ((*it)->getFamily() == METALLIC_FAMILY) {
297	dynamic_cast<MetallicInteraction>(it)->calcFunctional(fdat);
298	}
299	}
300
301	return;
302	}
303
304	void InteractionManager::doPair(int atid1, int atid2, RealType d, RealType r, RealType r2, RealType rcut, RealType sw, int topoDist, RealType pot, RealType vpair, RealType f1, RealType eFrame1, RealType eFrame2, RealType A1, RealType A2, RealType t1, RealType t2, RealType rho1, RealType rho2, RealType dfrho1, RealType dfrho2, RealType fshift1, RealType *fshift2){
305
306	if (!initialized_) initialize();
307
308	InteractionData idat;
309
310	idat.atype1 = typeMap_[*atid1];
311	idat.atype2 = typeMap_[*atid2];
312	idat.d = Vector3d(d);
313	idat.rij = *r;
314	idat.r2 = *r2;
315	idat.rcut = *rcut;
316	idat.sw = *sw;
317	idat.vdwMult = vdwScale_[*topoDist];
318	idat.electroMult = electrostaticScale_[*topoDist];
319	idat.pot = *pot;
320	idat.vpair = *vpair;
321	idat.f1 = Vector3d(f1);
322	idat.eFrame1 = Mat3x3d(eFrame1);
323	idat.eFrame2 = Mat3x3d(eFrame2);
324	idat.A1 = RotMat3x3d(A1);
325	idat.A2 = RotMat3x3d(A2);
326	idat.t1 = Vector3d(t1);
327	idat.t2 = Vector3d(t2);
328	idat.rho1 = *rho1;
329	idat.rho2 = *rho2;
330	idat.dfrho1 = *dfrho1;
331	idat.dfrho2 = *dfrho2;
332	idat.fshift1 = *fshift1;
333	idat.fshift2 = *fshift2;
334
335	pair<AtomType, AtomType> key = make_pair(idat.atype1, idat.atype2);
336	set<NonBondedInteraction*>::iterator it;
337
338	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
339	(*it)->calcForce(idat);
340
341	f1[0] = idat.f1.x();
342	f1[1] = idat.f1.y();
343	f1[2] = idat.f1.z();
344
345	t1[0] = idat.t1.x();
346	t1[1] = idat.t1.y();
347	t1[2] = idat.t1.z();
348
349	t2[0] = idat.t2.x();
350	t2[1] = idat.t2.y();
351	t2[2] = idat.t2.z();
352
353	return;
354	}
355
356	void InteractionManager::doSkipCorrection(int atid1, int atid2, RealType d, RealType r, RealType skippedCharge1, RealType skippedCharge2, RealType sw, RealType electroMult, RealType pot, RealType vpair, RealType f1, RealType eFrame1, RealType eFrame2, RealType t1, RealType *t2){
357
358	if (!initialized_) initialize();
359
360	SkipCorrectionData skdat;
361
362	skdat.atype1 = typeMap_[*atid1];
363	skdat.atype2 = typeMap_[*atid2];
364	skdat.d = Vector3d(d);
365	skdat.rij = *r;
366	skdat.skippedCharge1 = *skippedCharge1;
367	skdat.skippedCharge2 = *skippedCharge2;
368	skdat.sw = *sw;
369	skdat.electroMult = *electroMult;
370	skdat.pot = *pot;
371	skdat.vpair = *vpair;
372	skdat.f1 = Vector3d(f1);
373	skdat.eFrame1 = Mat3x3d(eFrame1);
374	skdat.eFrame2 = Mat3x3d(eFrame2);
375	skdat.t1 = Vector3d(t1);
376	skdat.t2 = Vector3d(t2);
377
378	pair<AtomType, AtomType> key = make_pair(skdat.atype1, skdat.atype2);
379	set<NonBondedInteraction*>::iterator it;
380
381	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
382	if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
383	dynamic_cast<ElectrostaticInteraction>(it)->calcSkipCorrection(skdat);
384	}
385	}
386
387	f1[0] = skdat.f1.x();
388	f1[1] = skdat.f1.y();
389	f1[2] = skdat.f1.z();
390
391	t1[0] = skdat.t1.x();
392	t1[1] = skdat.t1.y();
393	t1[2] = skdat.t1.z();
394
395	t2[0] = skdat.t2.x();
396	t2[1] = skdat.t2.y();
397	t2[2] = skdat.t2.z();
398
399	return;
400	}
401
402	void InteractionManager::doSelfCorrection(int atid, RealType eFrame, RealType skippedCharge, RealType pot, RealType *t){
403
404	if (!initialized_) initialize();
405
406	SelfCorrectionData scdat;
407
408	scdat.atype = typeMap_[*atid];
409	scdat.eFrame = Mat3x3d(eFrame);
410	scdat.skippedCharge = *skippedCharge;
411	scdat.pot = *pot;
412	scdat.t = Vector3d(t);
413
414	pair<AtomType, AtomType> key = make_pair(scdat.atype, scdat.atype);
415	set<NonBondedInteraction*>::iterator it;
416
417	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it){
418	if ((*it)->getFamily() == ELECTROSTATIC_FAMILY) {
419	dynamic_cast<ElectrostaticInteraction>(it)->calcSelfCorrection(scdat);
420	}
421	}
422
423	t[0] = scdat.t.x();
424	t[1] = scdat.t.y();
425	t[2] = scdat.t.z();
426
427	return;
428	}
429
430
431	RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
432	if (!initialized_) initialize();
433
434	AtomType* atype = typeMap_[*atid];
435
436	pair<AtomType, AtomType> key = make_pair(atype, atype);
437	set<NonBondedInteraction*>::iterator it;
438	RealType cutoff = 0.0;
439
440	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
441	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype));
442	return cutoff;
443	}
444
445	RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
446	if (!initialized_) initialize();
447
448	pair<AtomType, AtomType> key = make_pair(atype, atype);
449	set<NonBondedInteraction*>::iterator it;
450	RealType cutoff = 0.0;
451
452	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
453	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype));
454	return cutoff;
455	}
456
457
458	void InteractionManager::setSwitch(RealType rIn, RealType rOut) {
459	switcher_->setSwitch(rIn, rOut);
460	}
461
462	void InteractionManager::getSwitch(RealType r2, RealType sw, RealType dswdr, RealType r,
463	int *in_switching_region) {
464	bool isr = switcher_->getSwitch(r2, sw, dswdr, r);
465	*in_switching_region = (int)isr;
466	}
467
468	} //end namespace OpenMD
469
470	extern "C" {
471
472	#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR)
473	#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE)
474	#define fortranDoPair FC_FUNC(do_pair, DO_PAIR)
475	#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION)
476	#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION)
477	#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF)
478	#define fortranSetSwitch FC_FUNC(set_switch, SET_SWITCH)
479	#define fortranGetSwitch FC_FUNC(get_switch, GET_SWITCH)
480
481	void fortranDoPrePair(int atid1, int atid2, RealType *rij,
482	RealType rho_i_at_j, RealType rho_j_at_i) {
483
484	return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij,
485	rho_i_at_j,
486	rho_j_at_i);
487	}
488	void fortranDoPreForce(int atid, RealType rho, RealType *frho,
489	RealType *dfrhodrho) {
490
491	return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho,
492	dfrhodrho);
493	}
494
495	void fortranDoPair(int atid1, int atid2, RealType d, RealType r,
496	RealType r2, RealType rcut, RealType sw, int topoDist,
497	RealType pot, RealType vpair, RealType f1, RealType eFrame1,
498	RealType eFrame2, RealType A1, RealType *A2,
499	RealType t1, RealType t2, RealType rho1, RealType rho2,
500	RealType dfrho1, RealType dfrho2, RealType *fshift1,
501	RealType *fshift2){
502
503	return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r,
504	r2, rcut, sw, topoDist,
505	pot, vpair, f1,
506	eFrame1, eFrame2,
507	A1, A2, t1, t2, rho1,
508	rho2, dfrho1, dfrho2,
509	fshift1, fshift2);
510	}
511
512	void fortranDoSkipCorrection(int atid1, int atid2, RealType *d,
513	RealType r, RealType skippedCharge1,
514	RealType skippedCharge2, RealType sw,
515	RealType electroMult, RealType pot,
516	RealType vpair, RealType f1,
517	RealType eFrame1, RealType eFrame2,
518	RealType t1, RealType t2){
519
520	return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1,
521	atid2, d,
522	r,
523	skippedCharge1,
524	skippedCharge2,
525	sw, electroMult, pot,
526	vpair, f1, eFrame1,
527	eFrame2, t1, t2);
528	}
529
530	void fortranDoSelfCorrection(int atid, RealType eFrame, RealType *skippedCharge,
531	RealType pot, RealType t) {
532
533	return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid,
534	eFrame,
535	skippedCharge,
536	pot, t);
537	}
538	RealType fortranGetCutoff(int *atid) {
539	return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid);
540	}
541
542	void fortranGetSwitch(RealType r2, RealType sw, RealType dswdr, RealType r,
543	int *in_switching_region) {
544
545	return OpenMD::InteractionManager::Instance()->getSwitch(r2, sw, dswdr, r,
546	in_switching_region);
547	}
548
549	void fortranSetSwitch(RealType rIn, RealType rOut) {
550	return OpenMD::InteractionManager::Instance()->setSwitch(rIn, rOut);
551	}
552
553	}