src/nonbonded/InteractionManager.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
 * [4]  Vardeman & Gezelter, in progress (2009).                        
 */
 
#include "nonbonded/InteractionManager.hpp"

namespace OpenMD {
 
  bool InteractionManager::initialized_ = false;
  ForceField* InteractionManager::forceField_ = NULL;  
  InteractionManager* InteractionManager::_instance = NULL;
  map<int, AtomType*> InteractionManager::typeMap_;
  map<pair<AtomType*, AtomType*>, set<NonBondedInteraction*> > InteractionManager::interactions_;
 
  InteractionManager* InteractionManager::Instance() {
    if (!_instance) {
      _instance = new InteractionManager();
    }
    return _instance;
  }

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

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

    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
    ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
    AtomType* atype1;
    AtomType* atype2;
    pair<AtomType*, AtomType*> key;
    pair<set<NonBondedInteraction*>::iterator, bool> ret;
    
    for (atype1 = atomTypes->beginType(i1); atype1 != NULL; 
         atype1 = atomTypes->nextType(i1)) {
      
      // add it to the map:
      AtomTypeProperties atp = atype1->getATP();    
      
      pair<map<int,AtomType*>::iterator,bool> ret;    
      ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
      if (ret.second == false) {
        sprintf( painCave.errMsg,
                 "InteractionManager already had a previous entry with ident %d\n",
                 atp.ident);
        painCave.severity = OPENMD_INFO;
        painCave.isFatal = 0;
        simError();                 
      }
    }
    
    // Now, iterate over all known types and add to the interaction map:
    
    map<int, AtomType*>::iterator it1, it2;
    for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
      atype1 = (*it1).second;

      for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {        
        atype2 = (*it2).second;
        
        bool vdwExplicit = false;
        bool metExplicit = false;
        bool hbExplicit = false;
                       
        key = make_pair(atype1, atype2);
        
        if (atype1->isLennardJones() && atype2->isLennardJones()) {
          interactions_[key].insert(lj_);
        }
        if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
          interactions_[key].insert(electrostatic_);
        }
        if (atype1->isSticky() && atype2->isSticky() ) {
          interactions_[key].insert(sticky_);
        }
        if (atype1->isStickyPower() && atype2->isStickyPower() ) {
          interactions_[key].insert(sticky_);
        }
        if (atype1->isEAM() && atype2->isEAM() ) {
          interactions_[key].insert(eam_);
        }
        if (atype1->isSC() && atype2->isSC() ) {
          interactions_[key].insert(sc_);
        }
        if (atype1->isGayBerne() && atype2->isGayBerne() ) {
          interactions_[key].insert(gb_);
        }
        if ((atype1->isGayBerne() && atype2->isLennardJones())
            || (atype1->isLennardJones() && atype2->isGayBerne())) {
          interactions_[key].insert(gb_);
        } 
        
        // look for an explicitly-set non-bonded interaction type using the 
        // two atom types.
        NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());

        if (nbiType->isLennardJones()) {
          // We found an explicit Lennard-Jones interaction.  
          // override all other vdw entries for this pair of atom types:
          set<NonBondedInteraction*>::iterator it;
          for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
            InteractionFamily ifam = (*it)->getFamily();
            if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
          }
          interactions_[key].insert(lj_);
          vdwExplicit = true;
        }

        if (nbiType->isMorse()) {
          if (vdwExplicit) {
            sprintf( painCave.errMsg,
                     "InteractionManager::initialize found more than one explicit\n"
                     "\tvan der Waals interaction for atom types %s - %s\n",
                     atype1->getName().c_str(), atype2->getName().c_str());
            painCave.severity = OPENMD_ERROR;
            painCave.isFatal = 1;
            simError();
          }
          // We found an explicit Morse interaction.  
          // override all other vdw entries for this pair of atom types:
          set<NonBondedInteraction*>::iterator it;
          for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
            InteractionFamily ifam = (*it)->getFamily();
            if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
          }
          interactions_[key].insert(morse_);
          vdwExplicit = true;
        }

        if (nbiType->isEAM()) {
          // We found an explicit EAM interaction.  
          // override all other metallic entries for this pair of atom types:
          set<NonBondedInteraction*>::iterator it;
          for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
            InteractionFamily ifam = (*it)->getFamily();
            if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
          }
          interactions_[key].insert(eam_);
          metExplicit = true;
        }

        if (nbiType->isSC()) {
          if (metExplicit) {
            sprintf( painCave.errMsg,
                     "InteractionManager::initialize found more than one explicit\n"
                     "\tmetallic interaction for atom types %s - %s\n",
                     atype1->getName().c_str(), atype2->getName().c_str());
            painCave.severity = OPENMD_ERROR;
            painCave.isFatal = 1;
            simError();
          }
          // We found an explicit Sutton-Chen interaction.  
          // override all other metallic entries for this pair of atom types:
          set<NonBondedInteraction*>::iterator it;
          for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
            InteractionFamily ifam = (*it)->getFamily();
            if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
          }
          interactions_[key].insert(sc_);
          metExplicit = true;
        }
      }
    }
    
    // make sure every pair of atom types has a non-bonded interaction:
    for (atype1 = atomTypes->beginType(i1); atype1 != NULL; 
         atype1 = atomTypes->nextType(i1)) {
      for (atype2 = atomTypes->beginType(i2); atype2 != NULL; 
           atype2 = atomTypes->nextType(i2)) {
        key = make_pair(atype1, atype2);
        
        if (interactions_[key].size() == 0) {
          sprintf( painCave.errMsg,
                   "InteractionManager unable to find an appropriate non-bonded\n"
                   "\tinteraction for atom types %s - %s\n",
                   atype1->getName().c_str(), atype2->getName().c_str());
          painCave.severity = OPENMD_INFO;
          painCave.isFatal = 1;
          simError();
        }
      }
    }
  } 


  void InteractionManager::doPrePair(AtomType* atype1,
                                     AtomType* atype2,
                                     RealType rij,
                                     RealType &rho_i_at_j,
                                     RealType &rho_j_at_i) {
    
  }
  
  void InteractionManager::doPreForce(AtomType* atype,
                                      RealType rho,      
                                      RealType &frho,
                                      RealType &dfrhodrho) {
  }

  void InteractionManager::doSkipCorrection(AtomType* atype1,       
                                            AtomType* atype2,
                                            Vector3d d,
                                            RealType rij,
                                            RealType &skippedCharge1,
                                            RealType &skippedCharge2,
                                            RealType sw,
                                            RealType electroMult,
                                            RealType &pot,
                                            RealType &vpair,
                                            Vector3d &f1,
                                            Mat3x3d eFrame1,
                                            Mat3x3d eFrame2,
                                            Vector3d &t1,
                                            Vector3d &t2) {
  }
  
  void InteractionManager::doSelfCorrection(AtomType* atype,
                                            Mat3x3d eFrame,
                                            RealType skippedCharge,
                                            RealType &pot,
                                            Vector3d &t) {
  }

  void InteractionManager::do_prepair(int *atid1, int *atid2, RealType *rij, RealType *rho_i_at_j, RealType *rho_j_at_i){
    
    if (!initialized_) initialize();
    AtomType* atype1 = typeMap_[*atid1];
    AtomType* atype2 = typeMap_[*atid2];
    
    doPrePair(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
    
    return;    
  }

  void InteractionManager::do_preforce(int *atid, RealType *rho, RealType *frho, RealType *dfrhodrho){

    if (!initialized_) initialize();    
    AtomType* atype = typeMap_[*atid];

    doPreForce(atype, *rho, *frho, *dfrhodrho);
    
    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();
           
    AtomType* atype1 = typeMap_[*atid1];
    AtomType* atype2 = typeMap_[*atid2];
    Vector3d disp(d);
    Vector3d frc(f1);
    Vector3d trq1(t1);
    Vector3d trq2(t2);
    Mat3x3d eFi(eFrame1);
    Mat3x3d eFj(eFrame2);
       
    doSkipCorrection(atype1, atype2, disp, *r, *skippedCharge1, *skippedCharge2, *sw, 
                     *electroMult, *pot,  *vpair, frc, eFi, eFj, trq1, trq2);

    f1[0] = frc.x();
    f1[1] = frc.y();
    f1[2] = frc.z();
    
    t1[0] = trq1.x();
    t1[1] = trq1.y();
    t1[2] = trq1.z();
    
    t2[0] = trq2.x();
    t2[1] = trq2.y();
    t2[2] = trq2.z();

    return;
  }

  void InteractionManager::do_self_correction(int *atid, RealType *eFrame, RealType *skippedCharge, RealType *pot, RealType *t){

    if (!initialized_) initialize();
           
    AtomType* atype = typeMap_[*atid];
    Mat3x3d eFi(eFrame);
    Vector3d trq1(t);
       
    doSelfCorrection(atype, eFi, *skippedCharge, *pot, trq1);

    t[0] = trq1.x();
    t[1] = trq1.y();
    t[2] = trq1.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:	1502
Committed:	Sat Oct 2 19:53:32 2010 UTC (14 years, 7 months ago) by gezelter
File size:	18529 byte(s)
Log Message:	Many changes, and we're not quite done with this phase yet.
#	User	Rev	Content
1	gezelter	1502	/*
2			* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Redistributions of source code must retain the above copyright
10			* notice, this list of conditions and the following disclaimer.
11			*
12			* 2. Redistributions in binary form must reproduce the above copyright
13			* notice, this list of conditions and the following disclaimer in the
14			* documentation and/or other materials provided with the
15			* distribution.
16			*
17			* This software is provided "AS IS," without a warranty of any
18			* kind. All express or implied conditions, representations and
19			* warranties, including any implied warranty of merchantability,
20			* fitness for a particular purpose or non-infringement, are hereby
21			* excluded. The University of Notre Dame and its licensors shall not
22			* be liable for any damages suffered by licensee as a result of
23			* using, modifying or distributing the software or its
24			* derivatives. In no event will the University of Notre Dame or its
25			* licensors be liable for any lost revenue, profit or data, or for
26			* direct, indirect, special, consequential, incidental or punitive
27			* damages, however caused and regardless of the theory of liability,
28			* arising out of the use of or inability to use software, even if the
29			* University of Notre Dame has been advised of the possibility of
30			* such damages.
31			*
32			* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33			* research, please cite the appropriate papers when you publish your
34			* work. Good starting points are:
35			*
36			* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37			* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38			* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39			* [4] Vardeman & Gezelter, in progress (2009).
40			*/
41
42			#include "nonbonded/InteractionManager.hpp"
43
44			namespace OpenMD {
45
46			bool InteractionManager::initialized_ = false;
47			ForceField* InteractionManager::forceField_ = NULL;
48			InteractionManager* InteractionManager::_instance = NULL;
49			map<int, AtomType*> InteractionManager::typeMap_;
50			map<pair<AtomType, AtomType>, set<NonBondedInteraction*> > InteractionManager::interactions_;
51
52			InteractionManager* InteractionManager::Instance() {
53			if (!_instance) {
54			_instance = new InteractionManager();
55			}
56			return _instance;
57			}
58
59			void InteractionManager::initialize() {
60
61			lj_ = new LJ();
62			gb_ = new GB();
63			sticky_ = new Sticky();
64			eam_ = new EAM();
65			sc_ = new SC();
66			morse_ = new Morse();
67			electrostatic_ = new Electrostatic();
68
69			lj_->setForceField(forceField_);
70			gb_->setForceField(forceField_);
71			sticky_->setForceField(forceField_);
72			eam_->setForceField(forceField_);
73			sc_->setForceField(forceField_);
74			morse_->setForceField(forceField_);
75			electrostatic_->setForceField(forceField_);
76
77			ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
78			ForceField::AtomTypeContainer::MapTypeIterator i1, i2;
79			AtomType* atype1;
80			AtomType* atype2;
81			pair<AtomType, AtomType> key;
82			pair<set<NonBondedInteraction*>::iterator, bool> ret;
83
84			for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
85			atype1 = atomTypes->nextType(i1)) {
86
87			// add it to the map:
88			AtomTypeProperties atp = atype1->getATP();
89
90			pair<map<int,AtomType*>::iterator,bool> ret;
91			ret = typeMap_.insert( pair<int, AtomType*>(atp.ident, atype1) );
92			if (ret.second == false) {
93			sprintf( painCave.errMsg,
94			"InteractionManager already had a previous entry with ident %d\n",
95			atp.ident);
96			painCave.severity = OPENMD_INFO;
97			painCave.isFatal = 0;
98			simError();
99			}
100			}
101
102			// Now, iterate over all known types and add to the interaction map:
103
104			map<int, AtomType*>::iterator it1, it2;
105			for (it1 = typeMap_.begin(); it1 != typeMap_.end(); ++it1) {
106			atype1 = (*it1).second;
107
108			for( it2 = typeMap_.begin(); it2 != typeMap_.end(); ++it2) {
109			atype2 = (*it2).second;
110
111			bool vdwExplicit = false;
112			bool metExplicit = false;
113			bool hbExplicit = false;
114
115			key = make_pair(atype1, atype2);
116
117			if (atype1->isLennardJones() && atype2->isLennardJones()) {
118			interactions_[key].insert(lj_);
119			}
120			if (atype1->isElectrostatic() && atype2->isElectrostatic() ) {
121			interactions_[key].insert(electrostatic_);
122			}
123			if (atype1->isSticky() && atype2->isSticky() ) {
124			interactions_[key].insert(sticky_);
125			}
126			if (atype1->isStickyPower() && atype2->isStickyPower() ) {
127			interactions_[key].insert(sticky_);
128			}
129			if (atype1->isEAM() && atype2->isEAM() ) {
130			interactions_[key].insert(eam_);
131			}
132			if (atype1->isSC() && atype2->isSC() ) {
133			interactions_[key].insert(sc_);
134			}
135			if (atype1->isGayBerne() && atype2->isGayBerne() ) {
136			interactions_[key].insert(gb_);
137			}
138			if ((atype1->isGayBerne() && atype2->isLennardJones())
139			\|\| (atype1->isLennardJones() && atype2->isGayBerne())) {
140			interactions_[key].insert(gb_);
141			}
142
143			// look for an explicitly-set non-bonded interaction type using the
144			// two atom types.
145			NonBondedInteractionType* nbiType = forceField_->getNonBondedInteractionType(atype1->getName(), atype2->getName());
146
147			if (nbiType->isLennardJones()) {
148			// We found an explicit Lennard-Jones interaction.
149			// override all other vdw entries for this pair of atom types:
150			set<NonBondedInteraction*>::iterator it;
151			for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
152			InteractionFamily ifam = (*it)->getFamily();
153			if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
154			}
155			interactions_[key].insert(lj_);
156			vdwExplicit = true;
157			}
158
159			if (nbiType->isMorse()) {
160			if (vdwExplicit) {
161			sprintf( painCave.errMsg,
162			"InteractionManager::initialize found more than one explicit\n"
163			"\tvan der Waals interaction for atom types %s - %s\n",
164			atype1->getName().c_str(), atype2->getName().c_str());
165			painCave.severity = OPENMD_ERROR;
166			painCave.isFatal = 1;
167			simError();
168			}
169			// We found an explicit Morse interaction.
170			// override all other vdw entries for this pair of atom types:
171			set<NonBondedInteraction*>::iterator it;
172			for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
173			InteractionFamily ifam = (*it)->getFamily();
174			if (ifam == VANDERWAALS_FAMILY) interactions_[key].erase(*it);
175			}
176			interactions_[key].insert(morse_);
177			vdwExplicit = true;
178			}
179
180			if (nbiType->isEAM()) {
181			// We found an explicit EAM interaction.
182			// override all other metallic entries for this pair of atom types:
183			set<NonBondedInteraction*>::iterator it;
184			for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
185			InteractionFamily ifam = (*it)->getFamily();
186			if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
187			}
188			interactions_[key].insert(eam_);
189			metExplicit = true;
190			}
191
192			if (nbiType->isSC()) {
193			if (metExplicit) {
194			sprintf( painCave.errMsg,
195			"InteractionManager::initialize found more than one explicit\n"
196			"\tmetallic interaction for atom types %s - %s\n",
197			atype1->getName().c_str(), atype2->getName().c_str());
198			painCave.severity = OPENMD_ERROR;
199			painCave.isFatal = 1;
200			simError();
201			}
202			// We found an explicit Sutton-Chen interaction.
203			// override all other metallic entries for this pair of atom types:
204			set<NonBondedInteraction*>::iterator it;
205			for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it) {
206			InteractionFamily ifam = (*it)->getFamily();
207			if (ifam == METALLIC_FAMILY) interactions_[key].erase(*it);
208			}
209			interactions_[key].insert(sc_);
210			metExplicit = true;
211			}
212			}
213			}
214
215			// make sure every pair of atom types has a non-bonded interaction:
216			for (atype1 = atomTypes->beginType(i1); atype1 != NULL;
217			atype1 = atomTypes->nextType(i1)) {
218			for (atype2 = atomTypes->beginType(i2); atype2 != NULL;
219			atype2 = atomTypes->nextType(i2)) {
220			key = make_pair(atype1, atype2);
221
222			if (interactions_[key].size() == 0) {
223			sprintf( painCave.errMsg,
224			"InteractionManager unable to find an appropriate non-bonded\n"
225			"\tinteraction for atom types %s - %s\n",
226			atype1->getName().c_str(), atype2->getName().c_str());
227			painCave.severity = OPENMD_INFO;
228			painCave.isFatal = 1;
229			simError();
230			}
231			}
232			}
233			}
234
235
236			void InteractionManager::doPrePair(AtomType* atype1,
237			AtomType* atype2,
238			RealType rij,
239			RealType &rho_i_at_j,
240			RealType &rho_j_at_i) {
241
242			}
243
244			void InteractionManager::doPreForce(AtomType* atype,
245			RealType rho,
246			RealType &frho,
247			RealType &dfrhodrho) {
248			}
249
250			void InteractionManager::doSkipCorrection(AtomType* atype1,
251			AtomType* atype2,
252			Vector3d d,
253			RealType rij,
254			RealType &skippedCharge1,
255			RealType &skippedCharge2,
256			RealType sw,
257			RealType electroMult,
258			RealType &pot,
259			RealType &vpair,
260			Vector3d &f1,
261			Mat3x3d eFrame1,
262			Mat3x3d eFrame2,
263			Vector3d &t1,
264			Vector3d &t2) {
265			}
266
267			void InteractionManager::doSelfCorrection(AtomType* atype,
268			Mat3x3d eFrame,
269			RealType skippedCharge,
270			RealType &pot,
271			Vector3d &t) {
272			}
273
274			void InteractionManager::do_prepair(int atid1, int atid2, RealType rij, RealType rho_i_at_j, RealType *rho_j_at_i){
275
276			if (!initialized_) initialize();
277			AtomType* atype1 = typeMap_[*atid1];
278			AtomType* atype2 = typeMap_[*atid2];
279
280			doPrePair(atype1, atype2, rij, rho_i_at_j, *rho_j_at_i);
281
282			return;
283			}
284
285			void InteractionManager::do_preforce(int atid, RealType rho, RealType frho, RealType dfrhodrho){
286
287			if (!initialized_) initialize();
288			AtomType* atype = typeMap_[*atid];
289
290			doPreForce(atype, rho, frho, *dfrhodrho);
291
292			return;
293			}
294
295			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){
296
297			if (!initialized_) initialize();
298
299			InteractionData idat;
300
301			idat.atype1 = typeMap_[*atid1];
302			idat.atype2 = typeMap_[*atid2];
303			idat.d = Vector3d(d);
304			idat.rij = *r;
305			idat.r2 = *r2;
306			idat.rcut = *rcut;
307			idat.sw = *sw;
308			idat.vdwMult = *vdwMult;
309			idat.electroMult = *electroMult;
310			idat.pot = *pot;
311			idat.vpair = *vpair;
312			idat.f1 = Vector3d(f1);
313			idat.eFrame1 = Mat3x3d(eFrame1);
314			idat.eFrame2 = Mat3x3d(eFrame2);
315			idat.A1 = RotMat3x3d(A1);
316			idat.A2 = RotMat3x3d(A2);
317			idat.t1 = Vector3d(t1);
318			idat.t2 = Vector3d(t2);
319			idat.rho1 = *rho1;
320			idat.rho2 = *rho2;
321			idat.dfrho1 = *dfrho1;
322			idat.dfrho2 = *dfrho2;
323			idat.fshift1 = *fshift1;
324			idat.fshift2 = *fshift2;
325
326			pair<AtomType, AtomType> key = make_pair(idat.atype1, idat.atype2);
327			set<NonBondedInteraction*>::iterator it;
328
329			for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
330			(*it)->calcForce(idat);
331
332			f1[0] = idat.f1.x();
333			f1[1] = idat.f1.y();
334			f1[2] = idat.f1.z();
335
336			t1[0] = idat.t1.x();
337			t1[1] = idat.t1.y();
338			t1[2] = idat.t1.z();
339
340			t2[0] = idat.t2.x();
341			t2[1] = idat.t2.y();
342			t2[2] = idat.t2.z();
343
344			return;
345			}
346
347			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){
348
349			if (!initialized_) initialize();
350
351			AtomType* atype1 = typeMap_[*atid1];
352			AtomType* atype2 = typeMap_[*atid2];
353			Vector3d disp(d);
354			Vector3d frc(f1);
355			Vector3d trq1(t1);
356			Vector3d trq2(t2);
357			Mat3x3d eFi(eFrame1);
358			Mat3x3d eFj(eFrame2);
359
360			doSkipCorrection(atype1, atype2, disp, r, skippedCharge1, skippedCharge2, sw,
361			electroMult, pot, *vpair, frc, eFi, eFj, trq1, trq2);
362
363			f1[0] = frc.x();
364			f1[1] = frc.y();
365			f1[2] = frc.z();
366
367			t1[0] = trq1.x();
368			t1[1] = trq1.y();
369			t1[2] = trq1.z();
370
371			t2[0] = trq2.x();
372			t2[1] = trq2.y();
373			t2[2] = trq2.z();
374
375			return;
376			}
377
378			void InteractionManager::do_self_correction(int atid, RealType eFrame, RealType skippedCharge, RealType pot, RealType *t){
379
380			if (!initialized_) initialize();
381
382			AtomType* atype = typeMap_[*atid];
383			Mat3x3d eFi(eFrame);
384			Vector3d trq1(t);
385
386			doSelfCorrection(atype, eFi, skippedCharge, pot, trq1);
387
388			t[0] = trq1.x();
389			t[1] = trq1.y();
390			t[2] = trq1.z();
391
392			return;
393			}
394
395			} //end namespace OpenMD
396
397			extern "C" {
398
399			#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR)
400			#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE)
401			#define fortranDoPair FC_FUNC(do_pair, DO_PAIR)
402			#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION)
403			#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION)
404			#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF)
405
406			void fortranDoPrePair(int atid1, int atid2, RealType *rij,
407			RealType rho_i_at_j, RealType rho_j_at_i) {
408
409			return OpenMD::InteractionManager::Instance()->do_prepair(atid1, atid2, rij,
410			rho_i_at_j,
411			rho_j_at_i);
412			}
413			void fortranDoPreforce(int atid, RealType rho, RealType *frho,
414			RealType *dfrhodrho) {
415
416			return OpenMD::InteractionManager::Instance()->do_preforce(atid, rho, frho,
417			dfrhodrho);
418			}
419
420			void fortranDoPair(int atid1, int atid2, RealType d, RealType r,
421			RealType r2, RealType rcut, RealType sw, RealType vdwMult,
422			RealType electroMult, RealType pot, RealType vpair, RealType f1,
423			RealType eFrame1, RealType eFrame2, RealType A1, RealType A2,
424			RealType t1, RealType t2,
425			RealType rho1, RealType rho2, RealType dfrho1, RealType dfrho2,
426			RealType fshift1, RealType fshift2){
427
428			return OpenMD::InteractionManager::Instance()->do_pair(atid1, atid2, d, r, r2, rcut,
429			sw, vdwMult, electroMult, pot,
430			vpair, f1, eFrame1, eFrame2,
431			A1, A2, t1, t2, rho1, rho2,
432			dfrho1, dfrho2, fshift1, fshift2);
433			}
434
435			void fortranDoSkipCorrection(int atid1, int atid2, RealType d, RealType r,
436			RealType skippedCharge1, RealType skippedCharge2,
437			RealType sw, RealType electroMult, RealType *pot,
438			RealType vpair, RealType f1,
439			RealType eFrame1, RealType eFrame2,
440			RealType t1, RealType t2){
441
442			return OpenMD::InteractionManager::Instance()->do_skip_correction(atid1, atid2, d, r,
443			skippedCharge1,
444			skippedCharge2,
445			sw, electroMult, pot,
446			vpair, f1, eFrame1,
447			eFrame2, t1, t2);
448			}
449
450			void fortranDoSelfCorrection(int atid, RealType eFrame, RealType *skippedCharge,
451			RealType pot, RealType t) {
452
453			return OpenMD::InteractionManager::Instance()->do_self_correction(atid, eFrame,
454			skippedCharge,
455			pot, t);
456			}
457			RealType fortranGetCutoff() {
458			return OpenMD::InteractionManager::Instance()->getCutoff();
459			}
460			}