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.
#	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
46			bool InteractionManager::initialized_ = false;
47	gezelter	1529	RealType InteractionManager::rCut_ = -1.0;
48			RealType InteractionManager::rSwitch_ = -1.0;
49	gezelter	1502	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	gezelter	1528
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	gezelter	1530	SwitchingFunction* InteractionManager::switcher_ = new SwitchingFunction();
62	gezelter	1502
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	gezelter	1530	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	gezelter	1502	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	gezelter	1505	void InteractionManager::doPrePair(int atid1, int atid2, RealType rij, RealType rho_i_at_j, RealType *rho_j_at_i){
258	gezelter	1502
259			if (!initialized_) initialize();
260	gezelter	1504
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	gezelter	1502
278			return;
279			}
280	gezelter	1504
281	gezelter	1505	void InteractionManager::doPreForce(int atid, RealType rho, RealType frho, RealType dfrhodrho){
282	gezelter	1502
283	gezelter	1504	if (!initialized_) initialize();
284
285			FunctionalData fdat;
286	gezelter	1502
287	gezelter	1504	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	gezelter	1502
295	gezelter	1504	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	gezelter	1502	return;
302			}
303
304	gezelter	1530	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	gezelter	1504
306	gezelter	1502	if (!initialized_) initialize();
307	gezelter	1504
308	gezelter	1502	InteractionData idat;
309	gezelter	1504
310	gezelter	1502	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	gezelter	1530	idat.vdwMult = vdwScale_[*topoDist];
318			idat.electroMult = electrostaticScale_[*topoDist];
319	gezelter	1502	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	gezelter	1505	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	gezelter	1502
358			if (!initialized_) initialize();
359	gezelter	1504
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	gezelter	1502
378	gezelter	1504	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	gezelter	1502
387	gezelter	1504	f1[0] = skdat.f1.x();
388			f1[1] = skdat.f1.y();
389			f1[2] = skdat.f1.z();
390	gezelter	1502
391	gezelter	1504	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	gezelter	1502
399	gezelter	1504	return;
400	gezelter	1502	}
401
402	gezelter	1505	void InteractionManager::doSelfCorrection(int atid, RealType eFrame, RealType skippedCharge, RealType pot, RealType *t){
403	gezelter	1502
404			if (!initialized_) initialize();
405	gezelter	1504
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	gezelter	1502
414	gezelter	1504	pair<AtomType, AtomType> key = make_pair(scdat.atype, scdat.atype);
415			set<NonBondedInteraction*>::iterator it;
416	gezelter	1502
417	gezelter	1504	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	gezelter	1502	}
429
430	gezelter	1505
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	gezelter	1528	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	gezelter	1530
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	gezelter	1502	} //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	gezelter	1530	#define fortranSetSwitch FC_FUNC(set_switch, SET_SWITCH)
479			#define fortranGetSwitch FC_FUNC(get_switch, GET_SWITCH)
480	gezelter	1502
481			void fortranDoPrePair(int atid1, int atid2, RealType *rij,
482			RealType rho_i_at_j, RealType rho_j_at_i) {
483
484	gezelter	1505	return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij,
485			rho_i_at_j,
486			rho_j_at_i);
487	gezelter	1502	}
488	gezelter	1528	void fortranDoPreForce(int atid, RealType rho, RealType *frho,
489	gezelter	1502	RealType *dfrhodrho) {
490
491	gezelter	1505	return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho,
492			dfrhodrho);
493	gezelter	1502	}
494
495			void fortranDoPair(int atid1, int atid2, RealType d, RealType r,
496	gezelter	1530	RealType r2, RealType rcut, RealType sw, int topoDist,
497			RealType pot, RealType vpair, RealType f1, RealType eFrame1,
498	gezelter	1504	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	gezelter	1502
503	gezelter	1505	return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r,
504	gezelter	1530	r2, rcut, sw, topoDist,
505	gezelter	1505	pot, vpair, f1,
506			eFrame1, eFrame2,
507			A1, A2, t1, t2, rho1,
508			rho2, dfrho1, dfrho2,
509			fshift1, fshift2);
510	gezelter	1502	}
511	gezelter	1505
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	gezelter	1502	RealType vpair, RealType f1,
517			RealType eFrame1, RealType eFrame2,
518			RealType t1, RealType t2){
519
520	gezelter	1505	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	gezelter	1502	}
529	gezelter	1505
530	gezelter	1502	void fortranDoSelfCorrection(int atid, RealType eFrame, RealType *skippedCharge,
531			RealType pot, RealType t) {
532
533	gezelter	1505	return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid,
534			eFrame,
535			skippedCharge,
536			pot, t);
537	gezelter	1502	}
538	gezelter	1505	RealType fortranGetCutoff(int *atid) {
539			return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid);
540	gezelter	1502	}
541	gezelter	1530
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	gezelter	1502	}