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
#	Content
1	/*
2	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	*
4	* The University of Notre Dame grants you ("Licensee") a
5	* non-exclusive, royalty free, license to use, modify and
6	* redistribute this software in source and binary code form, provided
7	* that the following conditions are met:
8	*
9	* 1. Redistributions of source code must retain the above copyright
10	* notice, this list of conditions and the following disclaimer.
11	*
12	* 2. Redistributions in binary form must reproduce the above copyright
13	* notice, this list of conditions and the following disclaimer in the
14	* documentation and/or other materials provided with the
15	* distribution.
16	*
17	* This software is provided "AS IS," without a warranty of any
18	* kind. All express or implied conditions, representations and
19	* warranties, including any implied warranty of merchantability,
20	* fitness for a particular purpose or non-infringement, are hereby
21	* excluded. The University of Notre Dame and its licensors shall not
22	* be liable for any damages suffered by licensee as a result of
23	* using, modifying or distributing the software or its
24	* derivatives. In no event will the University of Notre Dame or its
25	* licensors be liable for any lost revenue, profit or data, or for
26	* direct, indirect, special, consequential, incidental or punitive
27	* damages, however caused and regardless of the theory of liability,
28	* arising out of the use of or inability to use software, even if the
29	* University of Notre Dame has been advised of the possibility of
30	* such damages.
31	*
32	* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33	* research, please cite the appropriate papers when you publish your
34	* work. Good starting points are:
35	*
36	* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39	* [4] Vardeman & Gezelter, in progress (2009).
40	*/
41
42	#include "nonbonded/InteractionManager.hpp"
43
44	namespace OpenMD {
45
46	bool InteractionManager::initialized_ = false;
47	RealType InteractionManager::rCut_ = -1.0;
48	RealType InteractionManager::rSwitch_ = -1.0;
49	ForceField* InteractionManager::forceField_ = NULL;
50	InteractionManager* InteractionManager::_instance = NULL;
51	map<int, AtomType*> InteractionManager::typeMap_;
52	map<pair<AtomType, AtomType>, set<NonBondedInteraction*> > InteractionManager::interactions_;
53
54	LJ* InteractionManager::lj_ = new LJ();
55	GB* InteractionManager::gb_ = new GB();
56	Sticky* InteractionManager::sticky_ = new Sticky();
57	Morse* InteractionManager::morse_ = new Morse();
58	EAM* InteractionManager::eam_ = new EAM();
59	SC* InteractionManager::sc_ = new SC();
60	Electrostatic* InteractionManager::electrostatic_ = new Electrostatic();
61
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	void InteractionManager::doPrePair(int atid1, int atid2, RealType rij, RealType rho_i_at_j, RealType *rho_j_at_i){
238
239	if (!initialized_) initialize();
240
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
258	return;
259	}
260
261	void InteractionManager::doPreForce(int atid, RealType rho, RealType frho, RealType dfrhodrho){
262
263	if (!initialized_) initialize();
264
265	FunctionalData fdat;
266
267	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
275	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	return;
282	}
283
284	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
286	if (!initialized_) initialize();
287
288	InteractionData idat;
289
290	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	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
338	if (!initialized_) initialize();
339
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
358	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
367	f1[0] = skdat.f1.x();
368	f1[1] = skdat.f1.y();
369	f1[2] = skdat.f1.z();
370
371	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
379	return;
380	}
381
382	void InteractionManager::doSelfCorrection(int atid, RealType eFrame, RealType skippedCharge, RealType pot, RealType *t){
383
384	if (!initialized_) initialize();
385
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
394	pair<AtomType, AtomType> key = make_pair(scdat.atype, scdat.atype);
395	set<NonBondedInteraction*>::iterator it;
396
397	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	}
409
410
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	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	} //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	return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij,
452	rho_i_at_j,
453	rho_j_at_i);
454	}
455	void fortranDoPreForce(int atid, RealType rho, RealType *frho,
456	RealType *dfrhodrho) {
457
458	return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho,
459	dfrhodrho);
460	}
461
462	void fortranDoPair(int atid1, int atid2, RealType d, RealType r,
463	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
471	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	}
480
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	RealType vpair, RealType f1,
486	RealType eFrame1, RealType eFrame2,
487	RealType t1, RealType t2){
488
489	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	}
498
499	void fortranDoSelfCorrection(int atid, RealType eFrame, RealType *skippedCharge,
500	RealType pot, RealType t) {
501
502	return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid,
503	eFrame,
504	skippedCharge,
505	pot, t);
506	}
507	RealType fortranGetCutoff(int *atid) {
508	return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid);
509	}
510	}