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();
 
  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_);

    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;    
  }

  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;    
  }

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