src/nonbonded/EAM.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 <stdio.h>
#include <string.h>

#include <cmath>
#include "nonbonded/EAM.hpp"
#include "utils/simError.h"


namespace OpenMD {

  bool EAM::initialized_ = false;
  ForceField* EAM::forceField_ = NULL;
  std::map<int, AtomType*> EAM::EAMlist;
  std::map<AtomType*, EAMAtomData> EAM::EAMMap;
  std::map<std::pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
  
  EAM* EAM::_instance = NULL;

  EAM* EAM::Instance() {
    if (!_instance) {
      _instance = new EAM();
    }
    return _instance;
  }
  
  EAMParam EAM::getEAMParam(AtomType* atomType) {
    
    // Do sanity checking on the AtomType we were passed before
    // building any data structures:
    if (!atomType->isEAM()) {
      sprintf( painCave.errMsg,
               "EAM::getEAMParam was passed an atomType (%s) that does not\n"
               "\tappear to be an embedded atom method (EAM) atom.\n",
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();
    }
    
    GenericData* data = atomType->getPropertyByName("EAM");
    if (data == NULL) {
      sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
               "\tparameters for atomType %s.\n", 
               atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError(); 
    }
    
    EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
    if (eamData == NULL) {
      sprintf( painCave.errMsg,
               "EAM::getEAMParam could not convert GenericData to EAMParam for\n"
               "\tatom type %s\n", atomType->getName().c_str());
      painCave.severity = OPENMD_ERROR;
      painCave.isFatal = 1;
      simError();          
    }
    
    return eamData->getData();
  }

  CubicSpline* EAM::getZ(AtomType* atomType) {    
    EAMParam eamParam = getEAMParam(atomType);
    int nr = eamParam.nr;
    RealType dr = eamParam.dr;
    vector<RealType> rvals;
    
    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
      
    CubicSpline* cs = new CubicSpline();
    cs->addPoints(rvals, eamParam.Z);
    return cs;
  }

  CubicSpline* EAM::getRho(AtomType* atomType) {    
    EAMParam eamParam = getEAMParam(atomType);
    int nr = eamParam.nr;
    RealType dr = eamParam.dr;
    vector<RealType> rvals;
    
    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
      
    CubicSpline* cs = new CubicSpline();
    cs->addPoints(rvals, eamParam.rho);
    return cs;
  }

  CubicSpline* EAM::getF(AtomType* atomType) {    
    EAMParam eamParam = getEAMParam(atomType);
    int nrho = eamParam.nrho;
    RealType drho = eamParam.drho;
    vector<RealType> rhovals;
    vector<RealType> scaledF;
    
    for (int i = 0; i < nrho; i++) {
      rhovals.push_back(i * drho);
      scaledF.push_back( eamParam.F[i] * 23.06054 );
    }
      
    CubicSpline* cs = new CubicSpline();
    cs->addPoints(rhovals, eamParam.F);
    return cs;
  }
  
  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {    
    EAMParam eamParam1 = getEAMParam(atomType1);
    EAMParam eamParam2 = getEAMParam(atomType2);
    CubicSpline* z1 = getZ(atomType1);
    CubicSpline* z2 = getZ(atomType2);

    // make the r grid:

    // set rcut to be the smaller of the two atomic rcuts

    RealType rcut = eamParam1.rcut < eamParam2.rcut ? 
      eamParam1.rcut : eamParam2.rcut; 

    // use the smallest dr (finest grid) to build our grid:

    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
    int nr = int(rcut/dr);
    vector<RealType> rvals;
    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);

    // construct the pair potential:

    vector<RealType> phivals;
    RealType phi;
    RealType r;
    RealType zi, zj;

    phivals.push_back(0.0);

    for (int i = 1; i < rvals.size(); i++ ) {
      r = rvals[i];
      zi = z1->getValueAt(r); 
      zj = z2->getValueAt(r);

      phi = 331.999296 * (zi * zj) / r;
      phivals.push_back(phi);
    }
      
    CubicSpline* cs = new CubicSpline();
    cs->addPoints(rvals, phivals);
    return cs;
  }

  void EAM::initialize() { 

    // set up the mixing method:
    ForceFieldOptions ffo = forceField_->getForceFieldOptions();
    string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());

    if (EAMMixMeth == "JOHNSON") 
      mixMeth_ = eamJohnson;    
    else if (EAMMixMeth == "DAW")
      mixMeth_ = eamDaw;
    else
      mixMeth_ = eamUnknown;
      
    // find all of the EAM atom Types:
    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
    ForceField::AtomTypeContainer::MapTypeIterator i;
    AtomType* at;

    for (at = atomTypes->beginType(i); at != NULL; 
         at = atomTypes->nextType(i)) {
      
      if (at->isEAM())
        addType(at);
    }
    
    // find all of the explicit EAM interactions (setfl):
    ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
    ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
    NonBondedInteractionType* nbt;

    for (nbt = nbiTypes->beginType(j); nbt != NULL; 
         nbt = nbiTypes->nextType(j)) {
      
      if (nbt->isEAM()) {
        
        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
        
        GenericData* data = nbt->getPropertyByName("EAM");
        if (data == NULL) {
          sprintf( painCave.errMsg, "EAM::rebuildMixingMap could not find\n"
                   "\tEAM parameters for %s - %s interaction.\n", 
                   atypes.first->getName().c_str(),
                   atypes.second->getName().c_str());
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError(); 
        }
        
        EAMMixingData* eamData = dynamic_cast<EAMMixingData*>(data);
        if (eamData == NULL) {
          sprintf( painCave.errMsg,
                   "EAM::rebuildMixingMap could not convert GenericData to\n"
                   "\tEAMMixingData for %s - %s interaction.\n", 
                   atypes.first->getName().c_str(),
                   atypes.second->getName().c_str());
          painCave.severity = OPENMD_ERROR;
          painCave.isFatal = 1;
          simError();          
        }
        
        EAMMix eamParam = eamData->getData();

        vector<RealType> phiAB = eamParam.phiAB;
        RealType dr = eamParam.dr;
        int nr = eamParam.nr;

        addExplicitInteraction(atypes.first, atypes.second, dr, nr, phiAB);
      }
    }  
    initialized_ = true;
  }
      

  void EAM::addType(AtomType* atomType){

    EAMAtomData eamAtomData;

    eamAtomData.rho = getRho(atomType);
    eamAtomData.F = getF(atomType);
    eamAtomData.Z = getZ(atomType);
    eamAtomData.rcut = getRcut(atomType);

    // add it to the map:
    AtomTypeProperties atp = atomType->getATP();    

    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
    if (ret.second == false) {
      sprintf( painCave.errMsg,
               "EAM already had a previous entry with ident %d\n",
               atp.ident);
      painCave.severity = OPENMD_INFO;
      painCave.isFatal = 0;
      simError();         
    }

    EAMMap[atomType] = eamAtomData;
    
    // Now, iterate over all known types and add to the mixing map:
    
    std::map<int, AtomType*>::iterator it;
    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
      
      AtomType* atype2 = (*it).second;

      EAMInteractionData mixer;
      mixer.phi = getPhi(atomType, atype2);
      mixer.explicitlySet = false;

      std::pair<AtomType*, AtomType*> key1, key2;
      key1 = std::make_pair(atomType, atype2);
      key2 = std::make_pair(atype2, atomType);
      
      MixingMap[key1] = mixer;
      if (key2 != key1) {
        MixingMap[key2] = mixer;
      }
    }      
    return;
  }
  
  void EAM::addExplicitInteraction(AtomType* atype1, AtomType* atype2, 
                                   RealType dr, int nr,
                                   vector<RealType> phiVals) {
    
    // in case these weren't already in the map
    addType(atype1);
    addType(atype2);

    EAMInteractionData mixer;
    CubicSpline* cs = new CubicSpline();
    vector<RealType> rvals;

    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);

    cs->addPoints(rVals, phiVals);
    mixer.phi = cs;
    mixer.explicitlySet = true;

    std::pair<AtomType*, AtomType*> key1, key2;
    key1 = std::make_pair(atype1, atype2);
    key2 = std::make_pair(atype2, atype1);
    
    MixingMap[key1] = mixer;
    if (key2 != key1) {
      MixingMap[key2] = mixer;
    }    
    return;
  }

  void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d, 
                        RealType rij, RealType r2, RealType rho_i_at_j,
                        RealType rho_j_at_i) {

    if (!initialized_) initialize();

    EAMAtomData data1 = EAMMap[at1];
    EAMAtomData data2 = EAMMap[at2];

    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
    return;
  }

  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
                           RealType dfrhodrho) {

    if (!initialized_) initialize();

    EAMAtomData data1 = EAMMap[at1];
        
    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);

    frho = result.first;
    dfrhodrho = result.second;
    return;
  }

 
  void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d, 
                      RealType rij, RealType r2, RealType sw, 
                      RealType &vpair, RealType &pot, Vector3d &f1,
                      RealType rho1, RealType rho2, RealType dfrho1,
                      RealType dfrho2, RealType fshift1, RealType fshift2) {

    if (!initialized_) initialize();
    
    pair<RealType, RealType> res;
    
    if (rij < eamRcut_) {

      EAMAtomData data1 = EAMMap[at1];
      EAMAtomData data2 = EAMMap[at2];

      // get type-specific cutoff radii

      RealType rci = data1.rcut;
      RealType rcj = data2.rcut;
      
      RealType rha, drha, rhb, drhb;
      RealType pha, dpha, phb, dphb;
      RealType phab, dvpdr;
      RealType drhoidr, drhojdr, dudr;
      
      if (rij < rci) {
        res = data1.rho->getValueAndDerivativeAt(rij);
        rha = res.first;
        drha = res.second;

        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
        pha = res.first;
        dpha = res.second;
      }

      if (rij < rcj) {
        res = data2.rho->getValueAndDerivativeAt(rij);
        rhb = res.first;
        drhb = res.second;

        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
        phb = res.first;
        dphb = res.second;
      }

      phab = 0.0;
      dvpdr = 0.0;

      switch(mixMeth_) {
      case eamJohnson:
       
        if (rij < rci) {
          phab = phab + 0.5 * (rhb / rha) * pha;
          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha + 
                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
        }

        if (rij < rcj) {
          phab = phab + 0.5 * (rha / rhb) * phb;
          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb + 
                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
        }

        break;

      case eamDaw:
                
        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
        phab = res.first;
        dvpdr = res.second;

        break;
      case eamUnknown:
      default:

        sprintf(painCave.errMsg,
                "EAM::calcForce hit a mixing method it doesn't know about!\n"
                );
        painCave.severity = OPENMD_ERROR;
        painCave.isFatal = 1;
        simError();        
          
      }
      
      drhoidr = drha;
      drhojdr = drhb;

      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr; 

      f1 = d * dudr / rij;
        
      // particle_pot is the difference between the full potential 
      // and the full potential without the presence of a particular
      // particle (atom1).
      //
      // This reduces the density at other particle locations, so
      // we need to recompute the density at atom2 assuming atom1
      // didn't contribute.  This then requires recomputing the
      // density functional for atom2 as well.
      //
      // Most of the particle_pot heavy lifting comes from the
      // pair interaction, and will be handled by vpair.
     
      fshift_i = data1.F->getValueAt( rho_i - rhb );
      fshift_j = data1.F->getValueAt( rho_j - rha );

      pot += phab;

      vpair += phab;
    }

    return;
    
  }


  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d, 
                                 RealType *rij, RealType *r2, 
                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
    if (!initialized_) initialize();

    AtomType* atype1 = EAMlist[*atid1];
    AtomType* atype2 = EAMlist[*atid2];
    
    Vector3d disp(d[0], d[1], d[2]);

    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);

    return;    
  }

  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
                                   RealType *dfrhodrho) {

    if (!initialized_) initialize();

    AtomType* atype1 = EAMlist[*atid1];   

    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
    
    return;    
  }

  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij, 
                        RealType *r2, RealType *sw, RealType *vpair, 
                        RealType *pot, RealType *f1, RealType *rho1,
                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
                        RealType *fshift1, RealType *fshift2) {

    if (!initialized_) initialize();
    
    AtomType* atype1 = EAMMap[*atid1];
    AtomType* atype2 = EAMMap[*atid2];
    
    Vector3d disp(d[0], d[1], d[2]);
    Vector3d frc(f1[0], f1[1], f1[2]);
    
    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
      
    f1[0] = frc.x();
    f1[1] = frc.y();
    f1[2] = frc.z();

    return;    
  }
  
  void EAM::setCutoffEAM(RealType *thisRcut) {
    eamRcut_ = thisRcut;
  }
}

extern "C" {
  
#define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
#define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
#define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
#define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
  
  RealType fortranCalcDensity(int *atid1, int *atid2, RealType *d, 
                              RealType *rij, RealType *r2, 
                              RealType *rho_i_at_j, RealType *rho_j_at_i) {

    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d, 
                                                         *rij, *r2, 
                                                         *rho_i_at_j,  
                                                         *rho_j_at_i);
  }
  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
                                 RealType *dfrhodrho) {  

    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1, 
                                                           *rho, 
                                                           *frho,
                                                           *dfrhodrho);

  }
  void fortranSetEAMCutoff(RealType *rcut) {
    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
  }
  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij, 
                        RealType *r2, RealType *sw, RealType *vpair, 
                        RealType *pot, RealType *f1, RealType *rho1,
                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
                        RealType *fshift1, RealType *fshift2){
    
    return OpenMD::EAM::Instance()->do_eam_pair(*atid1, *atid2, *d, *rij, 
                                                *r2, *sw,  *vpair, 
                                                *pot, *f1,  *rho1,
                                                *rho2,  *dfrho1,  *dfrho2,
                                                *fshift1,  *fshift2);
  }
}
Revision:	1478
Committed:	Fri Jul 23 20:45:40 2010 UTC (14 years, 9 months ago) by gezelter
File size:	18547 byte(s)
Log Message:	busticated EAM
#	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 <stdio.h>
43	#include <string.h>
44
45	#include <cmath>
46	#include "nonbonded/EAM.hpp"
47	#include "utils/simError.h"
48
49
50	namespace OpenMD {
51
52	bool EAM::initialized_ = false;
53	ForceField* EAM::forceField_ = NULL;
54	std::map<int, AtomType*> EAM::EAMlist;
55	std::map<AtomType*, EAMAtomData> EAM::EAMMap;
56	std::map<std::pair<AtomType, AtomType>, EAMInteractionData> EAM::MixingMap;
57
58	EAM* EAM::_instance = NULL;
59
60	EAM* EAM::Instance() {
61	if (!_instance) {
62	_instance = new EAM();
63	}
64	return _instance;
65	}
66
67	EAMParam EAM::getEAMParam(AtomType* atomType) {
68
69	// Do sanity checking on the AtomType we were passed before
70	// building any data structures:
71	if (!atomType->isEAM()) {
72	sprintf( painCave.errMsg,
73	"EAM::getEAMParam was passed an atomType (%s) that does not\n"
74	"\tappear to be an embedded atom method (EAM) atom.\n",
75	atomType->getName().c_str());
76	painCave.severity = OPENMD_ERROR;
77	painCave.isFatal = 1;
78	simError();
79	}
80
81	GenericData* data = atomType->getPropertyByName("EAM");
82	if (data == NULL) {
83	sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
84	"\tparameters for atomType %s.\n",
85	atomType->getName().c_str());
86	painCave.severity = OPENMD_ERROR;
87	painCave.isFatal = 1;
88	simError();
89	}
90
91	EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
92	if (eamData == NULL) {
93	sprintf( painCave.errMsg,
94	"EAM::getEAMParam could not convert GenericData to EAMParam for\n"
95	"\tatom type %s\n", atomType->getName().c_str());
96	painCave.severity = OPENMD_ERROR;
97	painCave.isFatal = 1;
98	simError();
99	}
100
101	return eamData->getData();
102	}
103
104	CubicSpline* EAM::getZ(AtomType* atomType) {
105	EAMParam eamParam = getEAMParam(atomType);
106	int nr = eamParam.nr;
107	RealType dr = eamParam.dr;
108	vector<RealType> rvals;
109
110	for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
111
112	CubicSpline* cs = new CubicSpline();
113	cs->addPoints(rvals, eamParam.Z);
114	return cs;
115	}
116
117	CubicSpline* EAM::getRho(AtomType* atomType) {
118	EAMParam eamParam = getEAMParam(atomType);
119	int nr = eamParam.nr;
120	RealType dr = eamParam.dr;
121	vector<RealType> rvals;
122
123	for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
124
125	CubicSpline* cs = new CubicSpline();
126	cs->addPoints(rvals, eamParam.rho);
127	return cs;
128	}
129
130	CubicSpline* EAM::getF(AtomType* atomType) {
131	EAMParam eamParam = getEAMParam(atomType);
132	int nrho = eamParam.nrho;
133	RealType drho = eamParam.drho;
134	vector<RealType> rhovals;
135	vector<RealType> scaledF;
136
137	for (int i = 0; i < nrho; i++) {
138	rhovals.push_back(i * drho);
139	scaledF.push_back( eamParam.F[i] * 23.06054 );
140	}
141
142	CubicSpline* cs = new CubicSpline();
143	cs->addPoints(rhovals, eamParam.F);
144	return cs;
145	}
146
147	CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {
148	EAMParam eamParam1 = getEAMParam(atomType1);
149	EAMParam eamParam2 = getEAMParam(atomType2);
150	CubicSpline* z1 = getZ(atomType1);
151	CubicSpline* z2 = getZ(atomType2);
152
153	// make the r grid:
154
155	// set rcut to be the smaller of the two atomic rcuts
156
157	RealType rcut = eamParam1.rcut < eamParam2.rcut ?
158	eamParam1.rcut : eamParam2.rcut;
159
160	// use the smallest dr (finest grid) to build our grid:
161
162	RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
163	int nr = int(rcut/dr);
164	vector<RealType> rvals;
165	for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
166
167	// construct the pair potential:
168
169	vector<RealType> phivals;
170	RealType phi;
171	RealType r;
172	RealType zi, zj;
173
174	phivals.push_back(0.0);
175
176	for (int i = 1; i < rvals.size(); i++ ) {
177	r = rvals[i];
178	zi = z1->getValueAt(r);
179	zj = z2->getValueAt(r);
180
181	phi = 331.999296 * (zi * zj) / r;
182	phivals.push_back(phi);
183	}
184
185	CubicSpline* cs = new CubicSpline();
186	cs->addPoints(rvals, phivals);
187	return cs;
188	}
189
190	void EAM::initialize() {
191
192	// set up the mixing method:
193	ForceFieldOptions ffo = forceField_->getForceFieldOptions();
194	string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());
195
196	if (EAMMixMeth == "JOHNSON")
197	mixMeth_ = eamJohnson;
198	else if (EAMMixMeth == "DAW")
199	mixMeth_ = eamDaw;
200	else
201	mixMeth_ = eamUnknown;
202
203	// find all of the EAM atom Types:
204	ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
205	ForceField::AtomTypeContainer::MapTypeIterator i;
206	AtomType* at;
207
208	for (at = atomTypes->beginType(i); at != NULL;
209	at = atomTypes->nextType(i)) {
210
211	if (at->isEAM())
212	addType(at);
213	}
214
215	// find all of the explicit EAM interactions (setfl):
216	ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
217	ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
218	NonBondedInteractionType* nbt;
219
220	for (nbt = nbiTypes->beginType(j); nbt != NULL;
221	nbt = nbiTypes->nextType(j)) {
222
223	if (nbt->isEAM()) {
224
225	std::pair<AtomType, AtomType> atypes = nbt->getAtomTypes();
226
227	GenericData* data = nbt->getPropertyByName("EAM");
228	if (data == NULL) {
229	sprintf( painCave.errMsg, "EAM::rebuildMixingMap could not find\n"
230	"\tEAM parameters for %s - %s interaction.\n",
231	atypes.first->getName().c_str(),
232	atypes.second->getName().c_str());
233	painCave.severity = OPENMD_ERROR;
234	painCave.isFatal = 1;
235	simError();
236	}
237
238	EAMMixingData* eamData = dynamic_cast<EAMMixingData*>(data);
239	if (eamData == NULL) {
240	sprintf( painCave.errMsg,
241	"EAM::rebuildMixingMap could not convert GenericData to\n"
242	"\tEAMMixingData for %s - %s interaction.\n",
243	atypes.first->getName().c_str(),
244	atypes.second->getName().c_str());
245	painCave.severity = OPENMD_ERROR;
246	painCave.isFatal = 1;
247	simError();
248	}
249
250	EAMMix eamParam = eamData->getData();
251
252	vector<RealType> phiAB = eamParam.phiAB;
253	RealType dr = eamParam.dr;
254	int nr = eamParam.nr;
255
256	addExplicitInteraction(atypes.first, atypes.second, dr, nr, phiAB);
257	}
258	}
259	initialized_ = true;
260	}
261
262
263
264	void EAM::addType(AtomType* atomType){
265
266	EAMAtomData eamAtomData;
267
268	eamAtomData.rho = getRho(atomType);
269	eamAtomData.F = getF(atomType);
270	eamAtomData.Z = getZ(atomType);
271	eamAtomData.rcut = getRcut(atomType);
272
273	// add it to the map:
274	AtomTypeProperties atp = atomType->getATP();
275
276	std::pair<std::map<int,AtomType*>::iterator,bool> ret;
277	ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, atomType) );
278	if (ret.second == false) {
279	sprintf( painCave.errMsg,
280	"EAM already had a previous entry with ident %d\n",
281	atp.ident);
282	painCave.severity = OPENMD_INFO;
283	painCave.isFatal = 0;
284	simError();
285	}
286
287	EAMMap[atomType] = eamAtomData;
288
289	// Now, iterate over all known types and add to the mixing map:
290
291	std::map<int, AtomType*>::iterator it;
292	for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
293
294	AtomType* atype2 = (*it).second;
295
296	EAMInteractionData mixer;
297	mixer.phi = getPhi(atomType, atype2);
298	mixer.explicitlySet = false;
299
300	std::pair<AtomType, AtomType> key1, key2;
301	key1 = std::make_pair(atomType, atype2);
302	key2 = std::make_pair(atype2, atomType);
303
304	MixingMap[key1] = mixer;
305	if (key2 != key1) {
306	MixingMap[key2] = mixer;
307	}
308	}
309	return;
310	}
311
312	void EAM::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
313	RealType dr, int nr,
314	vector<RealType> phiVals) {
315
316	// in case these weren't already in the map
317	addType(atype1);
318	addType(atype2);
319
320	EAMInteractionData mixer;
321	CubicSpline* cs = new CubicSpline();
322	vector<RealType> rvals;
323
324	for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
325
326	cs->addPoints(rVals, phiVals);
327	mixer.phi = cs;
328	mixer.explicitlySet = true;
329
330	std::pair<AtomType, AtomType> key1, key2;
331	key1 = std::make_pair(atype1, atype2);
332	key2 = std::make_pair(atype2, atype1);
333
334	MixingMap[key1] = mixer;
335	if (key2 != key1) {
336	MixingMap[key2] = mixer;
337	}
338	return;
339	}
340
341	void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d,
342	RealType rij, RealType r2, RealType rho_i_at_j,
343	RealType rho_j_at_i) {
344
345	if (!initialized_) initialize();
346
347	EAMAtomData data1 = EAMMap[at1];
348	EAMAtomData data2 = EAMMap[at2];
349
350	if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
351	if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
352	return;
353	}
354
355	void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
356	RealType dfrhodrho) {
357
358	if (!initialized_) initialize();
359
360	EAMAtomData data1 = EAMMap[at1];
361
362	pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
363
364	frho = result.first;
365	dfrhodrho = result.second;
366	return;
367	}
368
369
370	void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
371	RealType rij, RealType r2, RealType sw,
372	RealType &vpair, RealType &pot, Vector3d &f1,
373	RealType rho1, RealType rho2, RealType dfrho1,
374	RealType dfrho2, RealType fshift1, RealType fshift2) {
375
376	if (!initialized_) initialize();
377
378	pair<RealType, RealType> res;
379
380	if (rij < eamRcut_) {
381
382	EAMAtomData data1 = EAMMap[at1];
383	EAMAtomData data2 = EAMMap[at2];
384
385	// get type-specific cutoff radii
386
387	RealType rci = data1.rcut;
388	RealType rcj = data2.rcut;
389
390	RealType rha, drha, rhb, drhb;
391	RealType pha, dpha, phb, dphb;
392	RealType phab, dvpdr;
393	RealType drhoidr, drhojdr, dudr;
394
395	if (rij < rci) {
396	res = data1.rho->getValueAndDerivativeAt(rij);
397	rha = res.first;
398	drha = res.second;
399
400	res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
401	pha = res.first;
402	dpha = res.second;
403	}
404
405	if (rij < rcj) {
406	res = data2.rho->getValueAndDerivativeAt(rij);
407	rhb = res.first;
408	drhb = res.second;
409
410	res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
411	phb = res.first;
412	dphb = res.second;
413	}
414
415	phab = 0.0;
416	dvpdr = 0.0;
417
418	switch(mixMeth_) {
419	case eamJohnson:
420
421	if (rij < rci) {
422	phab = phab + 0.5 * (rhb / rha) * pha;
423	dvpdr = dvpdr + 0.5((rhb/rha)dpha +
424	pha((drhb/rha) - (rhbdrha/rha/rha)));
425	}
426
427	if (rij < rcj) {
428	phab = phab + 0.5 * (rha / rhb) * phb;
429	dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
430	phb((drha/rhb) - (rhadrhb/rhb/rhb)));
431	}
432
433	break;
434
435	case eamDaw:
436
437	res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
438	phab = res.first;
439	dvpdr = res.second;
440
441	break;
442	case eamUnknown:
443	default:
444
445	sprintf(painCave.errMsg,
446	"EAM::calcForce hit a mixing method it doesn't know about!\n"
447	);
448	painCave.severity = OPENMD_ERROR;
449	painCave.isFatal = 1;
450	simError();
451
452	}
453
454	drhoidr = drha;
455	drhojdr = drhb;
456
457	dudr = drhojdrdfrhodrho_i + drhoidrdfrhodrho_j + dvpdr;
458
459	f1 = d * dudr / rij;
460
461	// particle_pot is the difference between the full potential
462	// and the full potential without the presence of a particular
463	// particle (atom1).
464	//
465	// This reduces the density at other particle locations, so
466	// we need to recompute the density at atom2 assuming atom1
467	// didn't contribute. This then requires recomputing the
468	// density functional for atom2 as well.
469	//
470	// Most of the particle_pot heavy lifting comes from the
471	// pair interaction, and will be handled by vpair.
472
473	fshift_i = data1.F->getValueAt( rho_i - rhb );
474	fshift_j = data1.F->getValueAt( rho_j - rha );
475
476	pot += phab;
477
478	vpair += phab;
479	}
480
481	return;
482
483	}
484
485
486	void EAM::calc_eam_prepair_rho(int atid1, int atid2, RealType *d,
487	RealType rij, RealType r2,
488	RealType* rho_i_at_j, RealType* rho_j_at_i){
489	if (!initialized_) initialize();
490
491	AtomType* atype1 = EAMlist[*atid1];
492	AtomType* atype2 = EAMlist[*atid2];
493
494	Vector3d disp(d[0], d[1], d[2]);
495
496	calcDensity(atype1, atype2, disp, rij, r2, rho_i_at_j, rho_j_at_i);
497
498	return;
499	}
500
501	void EAM::calc_eam_preforce_Frho(int atid1, RealType rho, RealType *frho,
502	RealType *dfrhodrho) {
503
504	if (!initialized_) initialize();
505
506	AtomType* atype1 = EAMlist[*atid1];
507
508	calcFunctional(atype1, rho, frho, *dfrhodrho);
509
510	return;
511	}
512
513	void EAM::do_eam_pair(int atid1, int atid2, RealType d, RealType rij,
514	RealType r2, RealType sw, RealType *vpair,
515	RealType pot, RealType f1, RealType *rho1,
516	RealType rho2, RealType dfrho1, RealType *dfrho2,
517	RealType fshift1, RealType fshift2) {
518
519	if (!initialized_) initialize();
520
521	AtomType* atype1 = EAMMap[*atid1];
522	AtomType* atype2 = EAMMap[*atid2];
523
524	Vector3d disp(d[0], d[1], d[2]);
525	Vector3d frc(f1[0], f1[1], f1[2]);
526
527	calcForce(atype1, atype2, disp, rij, r2, sw, vpair, *pot, frc,
528	rho1, rho2, dfrho1, dfrho2, fshift1, fshift2);
529
530	f1[0] = frc.x();
531	f1[1] = frc.y();
532	f1[2] = frc.z();
533
534	return;
535	}
536
537	void EAM::setCutoffEAM(RealType *thisRcut) {
538	eamRcut_ = thisRcut;
539	}
540	}
541
542	extern "C" {
543
544	#define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
545	#define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
546	#define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
547	#define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
548
549	RealType fortranCalcDensity(int atid1, int atid2, RealType *d,
550	RealType rij, RealType r2,
551	RealType rho_i_at_j, RealType rho_j_at_i) {
552
553	return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, *d,
554	rij, r2,
555	*rho_i_at_j,
556	*rho_j_at_i);
557	}
558	RealType fortranCalcFunctional(int atid1, RealType rho, RealType *frho,
559	RealType *dfrhodrho) {
560
561	return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
562	*rho,
563	*frho,
564	*dfrhodrho);
565
566	}
567	void fortranSetEAMCutoff(RealType *rcut) {
568	return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
569	}
570	void fortranDoEAMPair(int atid1, int atid2, RealType d, RealType rij,
571	RealType r2, RealType sw, RealType *vpair,
572	RealType pot, RealType f1, RealType *rho1,
573	RealType rho2, RealType dfrho1, RealType *dfrho2,
574	RealType fshift1, RealType fshift2){
575
576	return OpenMD::EAM::Instance()->do_eam_pair(atid1, atid2, d, rij,
577	r2, sw, *vpair,
578	pot, f1, *rho1,
579	rho2, dfrho1, *dfrho2,
580	fshift1, fshift2);
581	}
582	}