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();
  MAW* InteractionManager::maw_ = new MAW();
  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_);
    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 = 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;
        }

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