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;
  ForceField* InteractionManager::forceField_ = NULL;  
  InteractionManager* InteractionManager::_instance = NULL;
  map<int, AtomType*> InteractionManager::typeMap_;
  map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_;
 
  InteractionManager* InteractionManager::Instance() {
    if (!_instance) {
      _instance = new InteractionManager();
    }
    return _instance;
  }

  void InteractionManager::initialize() {
    
    lj_ = new LJ();
    gb_ = new GB();
    sticky_ = new Sticky();
    eam_ = new EAM();
    sc_ = new SC();
    morse_ = new Morse();
    electrostatic_ = new Electrostatic();

    lj_->setForceField(forceField_);
    gb_->setForceField(forceField_);
    sticky_->setForceField(forceField_);
    eam_->setForceField(forceField_);
    sc_->setForceField(forceField_);
    morse_->setForceField(forceField_);
    electrostatic_->setForceField(forceField_);

    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, RealType *vdwMult,RealType *electroMult, 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 = *vdwMult;
    idat.electroMult = *electroMult;
    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;    
  }

} //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)

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