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

  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:	1528
Committed:	Fri Dec 17 20:11:05 2010 UTC (14 years, 4 months ago) by gezelter
File size:	19998 byte(s)
Log Message:	Doesn't build yet, but getting much closer...
#	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	LJ* InteractionManager::lj_ = new LJ();
53	GB* InteractionManager::gb_ = new GB();
54	Sticky* InteractionManager::sticky_ = new Sticky();
55	Morse* InteractionManager::morse_ = new Morse();
56	EAM* InteractionManager::eam_ = new EAM();
57	SC* InteractionManager::sc_ = new SC();
58	Electrostatic* InteractionManager::electrostatic_ = new Electrostatic();
59
60	InteractionManager* InteractionManager::Instance() {
61	if (!_instance) {
62	_instance = new InteractionManager();
63	}
64	return _instance;
65	}
66
67	void InteractionManager::initialize() {
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::doPrePair(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::doPreForce(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::doPair(int atid1, int atid2, RealType d, RealType r, RealType r2, RealType rcut, RealType sw, RealType vdwMult,RealType electroMult, RealType pot, RealType vpair, RealType f1, RealType eFrame1, RealType eFrame2, RealType A1, RealType A2, RealType t1, RealType t2, RealType rho1, RealType rho2, RealType dfrho1, RealType dfrho2, RealType fshift1, RealType fshift2){
283
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::doSkipCorrection(int atid1, int atid2, RealType d, RealType r, RealType skippedCharge1, RealType skippedCharge2, RealType sw, RealType electroMult, RealType pot, RealType vpair, RealType f1, RealType eFrame1, RealType eFrame2, RealType t1, RealType *t2){
335
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::doSelfCorrection(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
409	RealType InteractionManager::getSuggestedCutoffRadius(int *atid) {
410	if (!initialized_) initialize();
411
412	AtomType* atype = typeMap_[*atid];
413
414	pair<AtomType, AtomType> key = make_pair(atype, atype);
415	set<NonBondedInteraction*>::iterator it;
416	RealType cutoff = 0.0;
417
418	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
419	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype));
420	return cutoff;
421	}
422
423	RealType InteractionManager::getSuggestedCutoffRadius(AtomType* atype) {
424	if (!initialized_) initialize();
425
426	pair<AtomType, AtomType> key = make_pair(atype, atype);
427	set<NonBondedInteraction*>::iterator it;
428	RealType cutoff = 0.0;
429
430	for (it = interactions_[key].begin(); it != interactions_[key].end(); ++it)
431	cutoff = max(cutoff, (*it)->getSuggestedCutoffRadius(atype, atype));
432	return cutoff;
433	}
434
435	} //end namespace OpenMD
436
437	extern "C" {
438
439	#define fortranDoPrePair FC_FUNC(do_prepair, DO_PREPAIR)
440	#define fortranDoPreForce FC_FUNC(do_preforce, DO_PREFORCE)
441	#define fortranDoPair FC_FUNC(do_pair, DO_PAIR)
442	#define fortranDoSkipCorrection FC_FUNC(do_skip_correction, DO_SKIP_CORRECTION)
443	#define fortranDoSelfCorrection FC_FUNC(do_self_correction, DO_SELF_CORRECTION)
444	#define fortranGetCutoff FC_FUNC(get_cutoff, GET_CUTOFF)
445
446	void fortranDoPrePair(int atid1, int atid2, RealType *rij,
447	RealType rho_i_at_j, RealType rho_j_at_i) {
448
449	return OpenMD::InteractionManager::Instance()->doPrePair(atid1, atid2, rij,
450	rho_i_at_j,
451	rho_j_at_i);
452	}
453	void fortranDoPreForce(int atid, RealType rho, RealType *frho,
454	RealType *dfrhodrho) {
455
456	return OpenMD::InteractionManager::Instance()->doPreForce(atid, rho, frho,
457	dfrhodrho);
458	}
459
460	void fortranDoPair(int atid1, int atid2, RealType d, RealType r,
461	RealType r2, RealType rcut, RealType *sw,
462	RealType vdwMult, RealType electroMult, RealType *pot,
463	RealType vpair, RealType f1, RealType *eFrame1,
464	RealType eFrame2, RealType A1, RealType *A2,
465	RealType t1, RealType t2, RealType rho1, RealType rho2,
466	RealType dfrho1, RealType dfrho2, RealType *fshift1,
467	RealType *fshift2){
468
469	return OpenMD::InteractionManager::Instance()->doPair(atid1, atid2, d, r,
470	r2, rcut, sw,
471	vdwMult, electroMult,
472	pot, vpair, f1,
473	eFrame1, eFrame2,
474	A1, A2, t1, t2, rho1,
475	rho2, dfrho1, dfrho2,
476	fshift1, fshift2);
477	}
478
479	void fortranDoSkipCorrection(int atid1, int atid2, RealType *d,
480	RealType r, RealType skippedCharge1,
481	RealType skippedCharge2, RealType sw,
482	RealType electroMult, RealType pot,
483	RealType vpair, RealType f1,
484	RealType eFrame1, RealType eFrame2,
485	RealType t1, RealType t2){
486
487	return OpenMD::InteractionManager::Instance()->doSkipCorrection(atid1,
488	atid2, d,
489	r,
490	skippedCharge1,
491	skippedCharge2,
492	sw, electroMult, pot,
493	vpair, f1, eFrame1,
494	eFrame2, t1, t2);
495	}
496
497	void fortranDoSelfCorrection(int atid, RealType eFrame, RealType *skippedCharge,
498	RealType pot, RealType t) {
499
500	return OpenMD::InteractionManager::Instance()->doSelfCorrection(atid,
501	eFrame,
502	skippedCharge,
503	pot, t);
504	}
505	RealType fortranGetCutoff(int *atid) {
506	return OpenMD::InteractionManager::Instance()->getSuggestedCutoffRadius(atid);
507	}
508	}