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"
#include "types/NonBondedInteractionType.hpp"


namespace OpenMD {

  bool EAM::initialized_ = false;
  RealType EAM::eamRcut_ = 0.0;
  EAMMixingMethod EAM::mixMeth_ = eamJohnson;
  ForceField* EAM::forceField_ = NULL;
  map<int, AtomType*> EAM::EAMlist;
  map<AtomType*, EAMAtomData> EAM::EAMMap;
  map<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;
  }

  RealType EAM::getRcut(AtomType* atomType) {    
    EAMParam eamParam = getEAMParam(atomType);
    return eamParam.rcut;
  }

  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:


    // we need phi out to the largest value we'll encounter in the radial space;
    
    RealType rmax = 0.0;
    rmax = max(rmax, eamParam1.rcut);
    rmax = max(rmax, eamParam1.nr * eamParam1.dr);

    rmax = max(rmax, eamParam2.rcut);
    rmax = max(rmax, eamParam2.nr * eamParam2.dr);

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

    RealType dr = min(eamParam1.dr, eamParam2.dr); 

    int nr = int(rmax/dr + 0.5);

    vector<RealType> rvals;
    for (int i = 0; i < nr; i++) rvals.push_back(RealType(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];

      // only use z(r) if we're inside this atoms cutoff radius, otherwise, we'll use zero for the charge.
      // This effectively means that our phi grid goes out beyond the cutoff of the pair potential

      zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
      zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;

      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& fopts = forceField_->getForceFieldOptions();
    string EAMMixMeth = fopts.getEAMMixingMethod();
    toUpper(EAMMixMeth);
   
    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()) {
        
        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();          
        }
        
        EAMMixingParam eamParam = eamData->getData();

        vector<RealType> phiAB = eamParam.phi;
        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();    

    pair<map<int,AtomType*>::iterator,bool> ret;    
    ret = EAMlist.insert( 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:
    
    map<AtomType*, EAMAtomData>::iterator it;
    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
      
      AtomType* atype2 = (*it).first;

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

      pair<AtomType*, AtomType*> key1, key2;
      key1 = make_pair(atomType, atype2);
      key2 = 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;

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

  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij, 
                        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 rho_i, RealType rho_j, 
                      RealType dfrhodrho_i, RealType dfrhodrho_j, 
                      RealType &fshift_i, RealType &fshift_j) {

    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 *rij,
                                 RealType* rho_i_at_j, RealType* rho_j_at_i){

    if (!initialized_) initialize();
    
    AtomType* atype1 = EAMlist[*atid1];
    AtomType* atype2 = EAMlist[*atid2];
    
    calcDensity(atype1, atype2, *rij, *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;    
  }
  RealType EAM::getEAMcut(int *atid1) {

    if (!initialized_) initialize();
    
    AtomType* atype1 = EAMlist[*atid1];   
       
    return getRcut(atype1);
  }

  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 = EAMlist[*atid1];
    AtomType* atype2 = EAMlist[*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)
#define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)

  
  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij, 
                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
    
    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij, 
                                                         rho_i_at_j,  
                                                         rho_j_at_i);
  }
  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
                             RealType *dfrhodrho) {  
    
    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
                                                           dfrhodrho);
    
  }
  void fortranSetCutoffEAM(RealType *rcut) {
    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
  }
  void fortranCalcForce(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);
  }
  RealType fortranGetEAMcut(int* atid) {
    return OpenMD::EAM::Instance()->getEAMcut(atid);
  }

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