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::do_prepair(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::do_preforce(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::do_pair(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::do_skip_correction(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::do_self_correction(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;    
  }

} //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()->do_prepair(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()->do_preforce(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()->do_pair(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()->do_skip_correction(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()->do_self_correction(atid, eFrame, 
                                                                      skippedCharge,
                                                                      pot, t);
  }
  RealType fortranGetCutoff() {
    return OpenMD::InteractionManager::Instance()->getCutoff();
  }
}
Revision:	1504
Committed:	Sat Oct 2 20:41:53 2010 UTC (14 years, 7 months ago) by gezelter
File size:	18491 byte(s)
Log Message:	The C++ side now compiles. Moving on to doForces
#	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	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	void InteractionManager::do_prepair(int atid1, int atid2, RealType rij, RealType rho_i_at_j, RealType *rho_j_at_i){
236
237	if (!initialized_) initialize();
238
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
256	return;
257	}
258
259	void InteractionManager::do_preforce(int atid, RealType rho, RealType frho, RealType dfrhodrho){
260
261	if (!initialized_) initialize();
262
263	FunctionalData fdat;
264
265	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
273	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	return;
280	}
281
282	void InteractionManager::do_pair(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
284	if (!initialized_) initialize();
285
286	InteractionData idat;
287
288	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	void InteractionManager::do_skip_correction(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
336	if (!initialized_) initialize();
337
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
356	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
365	f1[0] = skdat.f1.x();
366	f1[1] = skdat.f1.y();
367	f1[2] = skdat.f1.z();
368
369	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
377	return;
378	}
379
380	void InteractionManager::do_self_correction(int atid, RealType eFrame, RealType skippedCharge, RealType pot, RealType *t){
381
382	if (!initialized_) initialize();
383
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
392	pair<AtomType, AtomType> key = make_pair(scdat.atype, scdat.atype);
393	set<NonBondedInteraction*>::iterator it;
394
395	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	}
407
408	} //end namespace OpenMD
409
410	extern "C" {
411
412	#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR)
413	#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE)
414	#define fortranDoPair FC_FUNC(do_pair, DO_PAIR)
415	#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION)
416	#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION)
417	#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF)
418
419	void fortranDoPrePair(int atid1, int atid2, RealType *rij,
420	RealType rho_i_at_j, RealType rho_j_at_i) {
421
422	return OpenMD::InteractionManager::Instance()->do_prepair(atid1, atid2, rij,
423	rho_i_at_j,
424	rho_j_at_i);
425	}
426	void fortranDoPreforce(int atid, RealType rho, RealType *frho,
427	RealType *dfrhodrho) {
428
429	return OpenMD::InteractionManager::Instance()->do_preforce(atid, rho, frho,
430	dfrhodrho);
431	}
432
433	void fortranDoPair(int atid1, int atid2, RealType d, RealType r,
434	RealType r2, RealType rcut, RealType *sw,
435	RealType vdwMult, RealType electroMult, RealType *pot,
436	RealType vpair, RealType f1, RealType *eFrame1,
437	RealType eFrame2, RealType A1, RealType *A2,
438	RealType t1, RealType t2, RealType rho1, RealType rho2,
439	RealType dfrho1, RealType dfrho2, RealType *fshift1,
440	RealType *fshift2){
441
442	return OpenMD::InteractionManager::Instance()->do_pair(atid1, atid2, d, r, r2, rcut,
443	sw, vdwMult, electroMult, pot,
444	vpair, f1, eFrame1, eFrame2,
445	A1, A2, t1, t2, rho1, rho2,
446	dfrho1, dfrho2, fshift1, fshift2);
447	}
448
449	void fortranDoSkipCorrection(int atid1, int atid2, RealType d, RealType r,
450	RealType skippedCharge1, RealType skippedCharge2,
451	RealType sw, RealType electroMult, RealType *pot,
452	RealType vpair, RealType f1,
453	RealType eFrame1, RealType eFrame2,
454	RealType t1, RealType t2){
455
456	return OpenMD::InteractionManager::Instance()->do_skip_correction(atid1, atid2, d, r,
457	skippedCharge1,
458	skippedCharge2,
459	sw, electroMult, pot,
460	vpair, f1, eFrame1,
461	eFrame2, t1, t2);
462	}
463
464	void fortranDoSelfCorrection(int atid, RealType eFrame, RealType *skippedCharge,
465	RealType pot, RealType t) {
466
467	return OpenMD::InteractionManager::Instance()->do_self_correction(atid, eFrame,
468	skippedCharge,
469	pot, t);
470	}
471	RealType fortranGetCutoff() {
472	return OpenMD::InteractionManager::Instance()->getCutoff();
473	}
474	}