UseTheForce/DarkSide/eam.F90

!!
!! Copyright (c) 2005, 2009 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).
!!
module eam
  use definitions
  use simulation
  use force_globals
  use status
  use atype_module
  use vector_class
  use interpolation
  implicit none
  PRIVATE
#define __FORTRAN90
#include "UseTheForce/DarkSide/fInteractionMap.h"

  logical, save :: EAM_FF_initialized = .false.
  integer, save :: EAM_Mixing_Policy
  real(kind = dp), save :: EAM_rcut
  logical, save :: haveRcut = .false.

  character(len = statusMsgSize) :: errMesg
  integer :: eam_err

  character(len = 200) :: errMsg
  character(len=*), parameter :: RoutineName =  "EAM MODULE"
  !! Logical that determines if eam arrays should be zeroed
  logical :: cleanme = .true.

  type, private :: EAMtype
     integer           :: eam_atype       
     real( kind = DP ) :: eam_lattice        
     real( kind = DP ) :: eam_rcut     
     integer           :: eam_atype_map
     type(cubicSpline) :: rho 
     type(cubicSpline) :: Z
     type(cubicSpline) :: F
     type(cubicSpline) :: phi
  end type EAMtype


  type, private :: EAMTypeList
     integer           :: n_eam_types = 0
     integer           :: currentAddition = 0

     type (EAMtype), pointer  :: EAMParams(:) => null()
     integer, pointer         :: atidToEAMType(:) => null()
  end type EAMTypeList

  type (eamTypeList), save :: EAMList

  !! standard eam stuff  


  public :: setCutoffEAM
  public :: do_eam_pair
  public :: newEAMtype
  public :: calc_eam_prepair_rho
  public :: calc_eam_preforce_Frho
  public :: destroyEAMTypes
  public :: getEAMCut
  public :: lookupEAMSpline
  public :: lookupEAMSpline1d

contains

  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho, c_ident, status)
    real (kind = dp )                      :: lattice_constant
    integer                                :: eam_nrho
    real (kind = dp )                      :: eam_drho
    integer                                :: eam_nr
    real (kind = dp )                      :: eam_dr
    real (kind = dp )                      :: rcut
    real (kind = dp ), dimension(eam_nr)   :: eam_Z_r, rvals
    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r, eam_phi_r
    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho, rhovals
    integer                                :: c_ident
    integer                                :: status

    integer                                :: nAtypes,nEAMTypes,myATID
    integer                                :: maxVals
    integer                                :: alloc_stat
    integer                                :: current, j
    integer,pointer                        :: Matchlist(:) => null()

    status = 0

    !! Assume that atypes has already been set and get the total number of types in atypes
    !! Also assume that every member of atypes is a EAM model.

    ! check to see if this is the first time into 
    if (.not.associated(EAMList%EAMParams)) then
       call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
       EAMList%n_eam_types = nEAMtypes
       allocate(EAMList%EAMParams(nEAMTypes))
       nAtypes = getSize(atypes)
       allocate(EAMList%atidToEAMType(nAtypes))
       EAMList%atidToEAMType = -1
    end if

    EAMList%currentAddition = EAMList%currentAddition + 1
    current = EAMList%currentAddition

    myATID =  getFirstMatchingElement(atypes, "c_ident", c_ident)
    EAMList%atidToEAMType(myATID) = current

    EAMList%EAMParams(current)%eam_atype    = c_ident
    EAMList%EAMParams(current)%eam_lattice  = lattice_constant
    EAMList%EAMParams(current)%eam_rcut     = rcut

    ! Build array of r values
    do j = 1, eam_nr
       rvals(j) = real(j-1,kind=dp) * eam_dr
    end do
    
    ! Build array of rho values
    do j = 1, eam_nrho
       rhovals(j) = real(j-1,kind=dp) * eam_drho
    end do

    ! convert from eV to kcal / mol:
    do j = 1, eam_nrho
       eam_F_rho(j) = eam_F_rho(j) * 23.06054E0_DP
    end do
    
    ! precompute the pair potential and get it into kcal / mol:
    eam_phi_r(1) = 0.0E0_DP
    do j = 2, eam_nr
       eam_phi_r(j) = 331.999296E0_DP * (eam_Z_r(j)**2) / rvals(j)
    enddo

    call newSpline(EAMList%EAMParams(current)%rho, rvals, eam_rho_r, .true.)
    call newSpline(EAMList%EAMParams(current)%Z,   rvals, eam_Z_r, .true.)
    call newSpline(EAMList%EAMParams(current)%F, rhovals, eam_F_rho, .true.)
    call newSpline(EAMList%EAMParams(current)%phi, rvals, eam_phi_r, .true.)
  end subroutine newEAMtype


  ! kills all eam types entered and sets simulation to uninitalized
  subroutine destroyEAMtypes()
    integer :: i
    type(EAMType), pointer :: tempEAMType=>null()

    do i = 1, EAMList%n_eam_types
       tempEAMType => eamList%EAMParams(i)
       call deallocate_EAMType(tempEAMType)
    end do
    if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
    eamList%EAMParams => null()

    eamList%n_eam_types = 0
    eamList%currentAddition = 0
  end subroutine destroyEAMtypes

  function getEAMCut(atomID) result(cutValue)
    integer, intent(in) :: atomID
    integer :: eamID
    real(kind=dp) :: cutValue
    
    eamID = EAMList%atidToEAMType(atomID)
    cutValue = EAMList%EAMParams(eamID)%eam_rcut
  end function getEAMCut


  subroutine setCutoffEAM(rcut)
    real(kind=dp) :: rcut
    EAM_rcut = rcut
  end subroutine setCutoffEAM


  subroutine deallocate_EAMType(thisEAMType)
    type (EAMtype), pointer :: thisEAMType

    call deleteSpline(thisEAMType%F)
    call deleteSpline(thisEAMType%rho)
    call deleteSpline(thisEAMType%phi)
    call deleteSpline(thisEAMType%Z)

  end subroutine deallocate_EAMType

  !! Calculates rho_r
  subroutine calc_eam_prepair_rho(atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i)
    integer :: Atid1, Atid2
    real(kind = dp), dimension(3) :: d
    real(kind = dp), intent(inout)               :: r
    real(kind = dp), intent(inout)               :: rijsq
    ! value of electron density rho do to atom i at atom j
    real(kind = dp), intent(inout) :: rho_i_at_j
    ! value of electron density rho do to atom j at atom i
    real(kind = dp), intent(inout) :: rho_j_at_i
    real(kind = dp) :: rci, rcj
    integer :: eam_err
   
    integer :: myid_atom1 ! EAM atid
    integer :: myid_atom2 

    ! check to see if we need to be cleaned at the start of a force loop

    Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
    Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
    
    if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
       
       ! get cutoff for atom 1
       rci = EAMList%EAMParams(myid_atom1)%eam_rcut
       ! get type specific cutoff for atom 2
       rcj = EAMList%EAMParams(myid_atom2)%eam_rcut
       
       if (r.lt.rci) then
          call lookupEAMSpline(EAMList%EAMParams(myid_atom1)%rho, r, &
               rho_i_at_j)
       endif

       if (r.lt.rcj) then
          call lookupEAMSpline(EAMList%EAMParams(myid_atom2)%rho, r, &
               rho_j_at_i)
       endif
    endif
  end subroutine calc_eam_prepair_rho


  !! Calculate the functional F(rho) for all local atoms
  subroutine calc_eam_preforce_Frho(nlocal, pot, particle_pot)
    integer :: nlocal
    real(kind=dp) :: pot
    integer :: i, j
    integer :: atom
    real(kind=dp) :: U,U1
    real( kind = dp ), dimension(nlocal) :: particle_pot
    integer :: atype1
    integer :: me, atid1

!    cleanme = .true.
    !! Calculate F(rho) and derivative 
    do atom = 1, nlocal
       atid1 = atid(atom)
       me = eamList%atidToEAMtype(atid1)

       ! Myid is set to -1 for non EAM atoms.
       ! Punt if we are a non-EAM atom type.
       if (me == -1) then
          frho(atom) = 0.0_dp
          dfrhodrho(atom) = 0.0_dp
       else

       call lookupEAMSpline1d(EAMList%EAMParams(me)%F, rho(atom), &
            u, u1)
       
          frho(atom) = u
          dfrhodrho(atom) = u1
       endif
       pot = pot + u
       particle_pot(atom) = particle_pot(atom) + u

    enddo

  end subroutine calc_eam_preforce_Frho
  
  !! Does EAM pairwise Force calculation.  
  subroutine do_eam_pair(atid1, atid2, d, rij, r2, sw, vpair, &
       fpair, pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, &
       fshift_i, fshift_j)
    !Arguments    
    integer, intent(in) ::  atid1, atid2
    real( kind = dp ), intent(in) :: rij, r2
    real( kind = dp ) :: pot, sw, vpair
    real( kind = dp ), dimension(3) :: f1
    real( kind = dp ), intent(in), dimension(3) :: d
    real( kind = dp ), intent(inout), dimension(3) :: fpair
    real( kind = dp ), intent(inout) :: dfrhodrho_i, dfrhodrho_j
    real( kind = dp ), intent(inout) :: rho_i, rho_j
    real( kind = dp ), intent(inout):: fshift_i, fshift_j

    real( kind = dp ) :: drdx, drdy, drdz
    real( kind = dp ) :: phab, pha, dvpdr
    real( kind = dp ) :: rha, drha, dpha
    real( kind = dp ) :: rhb, drhb, dphb
    real( kind = dp ) :: dudr
    real( kind = dp ) :: rci, rcj
    real( kind = dp ) :: drhoidr, drhojdr
    real( kind = dp ) :: Fx, Fy, Fz
    real( kind = dp ) :: r, phb
    real( kind = dp ) :: u1, u2

    integer :: id1, id2
    integer  :: mytype_atom1
    integer  :: mytype_atom2

    phab = 0.0E0_DP
    dvpdr = 0.0E0_DP
    drha = 0.0
    drhb = 0.0

    if (rij .lt. EAM_rcut) then

       mytype_atom1 = EAMList%atidToEAMType(atid1)
       mytype_atom2 = EAMList%atidTOEAMType(atid2)


       ! get cutoff for atom 1
       rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
       ! get type specific cutoff for atom 2
       rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut

       drdx = d(1)/rij
       drdy = d(2)/rij
       drdz = d(3)/rij

       if (rij.lt.rci) then

          ! Calculate rho and drho for atom1

          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%rho, &
               rij, rha, drha)
          
          ! Calculate Phi(r) for atom1.
          
          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%phi, &
               rij, pha, dpha)

       endif

       if (rij.lt.rcj) then      

          ! Calculate rho and drho for atom2

          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%rho, &
               rij, rhb, drhb)

          ! Calculate Phi(r) for atom2.

          call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%phi, &
               rij, phb, dphb)

       endif

       if (rij.lt.rci) then 
          phab = phab + 0.5E0_DP*(rhb/rha)*pha
          dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + &
               pha*((drhb/rha) - (rhb*drha/rha/rha)))
       endif

       if (rij.lt.rcj) then
          phab = phab + 0.5E0_DP*(rha/rhb)*phb
          dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + &
               phb*((drha/rhb) - (rha*drhb/rhb/rhb)))
       endif

       drhoidr = drha
       drhojdr = drhb

       dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr        

       fx = dudr * drdx
       fy = dudr * drdy
       fz = dudr * drdz

       ! 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.
       
       call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%F, &
            rho_i-rhb, fshift_i, u1)
       call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%F, &
            rho_j-rha, fshift_j, u2)            

       pot = pot + phab

       f1(1) = f1(1) + fx
       f1(2) = f1(2) + fy
       f1(3) = f1(3) + fz

       vpair = vpair + phab

    endif
  end subroutine do_eam_pair

  subroutine lookupEAMSpline(cs, xval, yval)
    
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval
    real( kind = DP ) ::  dx
    integer :: i, j
    !
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.
    
    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
    
    dx = xval - cs%x(j)
    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))
    
    return
  end subroutine lookupEAMSpline

  subroutine lookupEAMSpline1d(cs, xval, yval, dydx)
    
    implicit none

    type (cubicSpline), intent(in) :: cs
    real( kind = DP ), intent(in)  :: xval
    real( kind = DP ), intent(out) :: yval, dydx
    real( kind = DP ) :: dx
    integer :: i, j
    
    !  Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains 
    !  or is nearest to xval.


    j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
    
    dx = xval - cs%x(j)
    yval = cs%y(j) + dx*(cs%b(j) + dx*(cs%c(j) + dx*cs%d(j)))

    dydx = cs%b(j) + dx*(2.0d0 * cs%c(j) + 3.0d0 * dx * cs%d(j))
       
    return
  end subroutine lookupEAMSpline1d

end module eam
Revision:	1604
Committed:	Mon Aug 8 18:53:40 2011 UTC (14 years ago) by jmichalk
File size:	15095 byte(s)
Log Message:	This commit should allow EADM simulations to be run on the cluster. Main additions include EADM_FF.hpp/cpp as well as a system dipole correlation option in DynamicProps.
#	User	Rev	Content
1	gezelter	246	!!
2	chuckv	1388	!! Copyright (c) 2005, 2009 The University of Notre Dame. All Rights Reserved.
3	gezelter	246	!!
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	gezelter	1390	!! 1. Redistributions of source code must retain the above copyright
10	gezelter	246	!! notice, this list of conditions and the following disclaimer.
11			!!
12	gezelter	1390	!! 2. Redistributions in binary form must reproduce the above copyright
13	gezelter	246	!! 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	gezelter	1390	!! 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	gezelter	115	module eam
42	gezelter	960	use definitions
43	gezelter	115	use simulation
44			use force_globals
45			use status
46			use atype_module
47	chrisfen	578	use vector_class
48	gezelter	938	use interpolation
49	gezelter	115	implicit none
50			PRIVATE
51	chuckv	656	#define __FORTRAN90
52			#include "UseTheForce/DarkSide/fInteractionMap.h"
53	gezelter	115
54			logical, save :: EAM_FF_initialized = .false.
55			integer, save :: EAM_Mixing_Policy
56			real(kind = dp), save :: EAM_rcut
57			logical, save :: haveRcut = .false.
58
59			character(len = statusMsgSize) :: errMesg
60			integer :: eam_err
61
62			character(len = 200) :: errMsg
63			character(len=*), parameter :: RoutineName = "EAM MODULE"
64	gezelter	507	!! Logical that determines if eam arrays should be zeroed
65	gezelter	115	logical :: cleanme = .true.
66
67			type, private :: EAMtype
68			integer :: eam_atype
69			real( kind = DP ) :: eam_lattice
70			real( kind = DP ) :: eam_rcut
71			integer :: eam_atype_map
72	gezelter	938	type(cubicSpline) :: rho
73			type(cubicSpline) :: Z
74			type(cubicSpline) :: F
75			type(cubicSpline) :: phi
76	gezelter	115	end type EAMtype
77
78
79			type, private :: EAMTypeList
80			integer :: n_eam_types = 0
81			integer :: currentAddition = 0
82	gezelter	507
83	gezelter	115	type (EAMtype), pointer :: EAMParams(:) => null()
84	chuckv	577	integer, pointer :: atidToEAMType(:) => null()
85	gezelter	115	end type EAMTypeList
86
87			type (eamTypeList), save :: EAMList
88
89			!! standard eam stuff
90
91
92			public :: setCutoffEAM
93			public :: do_eam_pair
94			public :: newEAMtype
95			public :: calc_eam_prepair_rho
96			public :: calc_eam_preforce_Frho
97	chuckv	472	public :: destroyEAMTypes
98	chrisfen	578	public :: getEAMCut
99	chrisfen	943	public :: lookupEAMSpline
100			public :: lookupEAMSpline1d
101	gezelter	115
102			contains
103
104			subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
105	gezelter	938	eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho, c_ident, status)
106	gezelter	115	real (kind = dp ) :: lattice_constant
107			integer :: eam_nrho
108			real (kind = dp ) :: eam_drho
109			integer :: eam_nr
110			real (kind = dp ) :: eam_dr
111			real (kind = dp ) :: rcut
112	gezelter	938	real (kind = dp ), dimension(eam_nr) :: eam_Z_r, rvals
113			real (kind = dp ), dimension(eam_nr) :: eam_rho_r, eam_phi_r
114			real (kind = dp ), dimension(eam_nrho) :: eam_F_rho, rhovals
115	chrisfen	579	integer :: c_ident
116	gezelter	115	integer :: status
117
118	chuckv	577	integer :: nAtypes,nEAMTypes,myATID
119	gezelter	115	integer :: maxVals
120			integer :: alloc_stat
121	gezelter	938	integer :: current, j
122	gezelter	115	integer,pointer :: Matchlist(:) => null()
123
124			status = 0
125
126			!! Assume that atypes has already been set and get the total number of types in atypes
127			!! Also assume that every member of atypes is a EAM model.
128
129			! check to see if this is the first time into
130			if (.not.associated(EAMList%EAMParams)) then
131	chuckv	577	call getMatchingElementList(atypes, "is_EAM", .true., nEAMtypes, MatchList)
132			EAMList%n_eam_types = nEAMtypes
133			allocate(EAMList%EAMParams(nEAMTypes))
134			nAtypes = getSize(atypes)
135			allocate(EAMList%atidToEAMType(nAtypes))
136	gezelter	1578	EAMList%atidToEAMType = -1
137	gezelter	115	end if
138
139			EAMList%currentAddition = EAMList%currentAddition + 1
140			current = EAMList%currentAddition
141
142	chrisfen	579	myATID = getFirstMatchingElement(atypes, "c_ident", c_ident)
143	chuckv	577	EAMList%atidToEAMType(myATID) = current
144	gezelter	507
145	chrisfen	579	EAMList%EAMParams(current)%eam_atype = c_ident
146	gezelter	115	EAMList%EAMParams(current)%eam_lattice = lattice_constant
147			EAMList%EAMParams(current)%eam_rcut = rcut
148
149	gezelter	938	! Build array of r values
150			do j = 1, eam_nr
151			rvals(j) = real(j-1,kind=dp) * eam_dr
152			end do
153
154			! Build array of rho values
155			do j = 1, eam_nrho
156			rhovals(j) = real(j-1,kind=dp) * eam_drho
157			end do
158
159			! convert from eV to kcal / mol:
160			do j = 1, eam_nrho
161			eam_F_rho(j) = eam_F_rho(j) * 23.06054E0_DP
162			end do
163
164			! precompute the pair potential and get it into kcal / mol:
165			eam_phi_r(1) = 0.0E0_DP
166			do j = 2, eam_nr
167			eam_phi_r(j) = 331.999296E0_DP * (eam_Z_r(j)**2) / rvals(j)
168			enddo
169
170	gezelter	939	call newSpline(EAMList%EAMParams(current)%rho, rvals, eam_rho_r, .true.)
171			call newSpline(EAMList%EAMParams(current)%Z, rvals, eam_Z_r, .true.)
172			call newSpline(EAMList%EAMParams(current)%F, rhovals, eam_F_rho, .true.)
173			call newSpline(EAMList%EAMParams(current)%phi, rvals, eam_phi_r, .true.)
174	gezelter	115	end subroutine newEAMtype
175
176
177	chuckv	472	! kills all eam types entered and sets simulation to uninitalized
178			subroutine destroyEAMtypes()
179			integer :: i
180			type(EAMType), pointer :: tempEAMType=>null()
181	gezelter	115
182	chuckv	472	do i = 1, EAMList%n_eam_types
183			tempEAMType => eamList%EAMParams(i)
184			call deallocate_EAMType(tempEAMType)
185			end do
186			if(associated( eamList%EAMParams)) deallocate( eamList%EAMParams)
187			eamList%EAMParams => null()
188	gezelter	507
189	chuckv	472	eamList%n_eam_types = 0
190			eamList%currentAddition = 0
191			end subroutine destroyEAMtypes
192
193	chrisfen	578	function getEAMCut(atomID) result(cutValue)
194			integer, intent(in) :: atomID
195	chrisfen	579	integer :: eamID
196	chrisfen	578	real(kind=dp) :: cutValue
197
198	chrisfen	579	eamID = EAMList%atidToEAMType(atomID)
199			cutValue = EAMList%EAMParams(eamID)%eam_rcut
200	chrisfen	578	end function getEAMCut
201
202	gezelter	115
203	gezelter	938	subroutine setCutoffEAM(rcut)
204	gezelter	115	real(kind=dp) :: rcut
205			EAM_rcut = rcut
206			end subroutine setCutoffEAM
207
208	gezelter	507
209	gezelter	115	subroutine deallocate_EAMType(thisEAMType)
210			type (EAMtype), pointer :: thisEAMType
211
212	gezelter	938	call deleteSpline(thisEAMType%F)
213			call deleteSpline(thisEAMType%rho)
214			call deleteSpline(thisEAMType%phi)
215			call deleteSpline(thisEAMType%Z)
216	gezelter	507
217	gezelter	115	end subroutine deallocate_EAMType
218
219	gezelter	507	!! Calculates rho_r
220	gezelter	1464	subroutine calc_eam_prepair_rho(atid1, atid2, d, r, rijsq, rho_i_at_j, rho_j_at_i)
221			integer :: Atid1, Atid2
222	gezelter	115	real(kind = dp), dimension(3) :: d
223			real(kind = dp), intent(inout) :: r
224			real(kind = dp), intent(inout) :: rijsq
225			! value of electron density rho do to atom i at atom j
226	chuckv	1388	real(kind = dp), intent(inout) :: rho_i_at_j
227	gezelter	115	! value of electron density rho do to atom j at atom i
228	chuckv	1388	real(kind = dp), intent(inout) :: rho_j_at_i
229	gezelter	1573	real(kind = dp) :: rci, rcj
230	gezelter	115	integer :: eam_err
231	gezelter	1386
232	chuckv	592	integer :: myid_atom1 ! EAM atid
233			integer :: myid_atom2
234	gezelter	507
235			! check to see if we need to be cleaned at the start of a force loop
236	gezelter	115
237	chuckv	592	Myid_atom1 = Eamlist%atidtoeamtype(Atid1)
238			Myid_atom2 = Eamlist%atidtoeamtype(Atid2)
239	gezelter	1386
240	gezelter	115	if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
241	gezelter	1386
242	gezelter	1573	! get cutoff for atom 1
243			rci = EAMList%EAMParams(myid_atom1)%eam_rcut
244			! get type specific cutoff for atom 2
245			rcj = EAMList%EAMParams(myid_atom2)%eam_rcut
246
247			if (r.lt.rci) then
248			call lookupEAMSpline(EAMList%EAMParams(myid_atom1)%rho, r, &
249			rho_i_at_j)
250			endif
251	gezelter	115
252	gezelter	1573	if (r.lt.rcj) then
253			call lookupEAMSpline(EAMList%EAMParams(myid_atom2)%rho, r, &
254			rho_j_at_i)
255			endif
256	gezelter	507	endif
257	gezelter	115	end subroutine calc_eam_prepair_rho
258
259
260			!! Calculate the functional F(rho) for all local atoms
261	gezelter	1285	subroutine calc_eam_preforce_Frho(nlocal, pot, particle_pot)
262	gezelter	115	integer :: nlocal
263			real(kind=dp) :: pot
264	gezelter	938	integer :: i, j
265	gezelter	115	integer :: atom
266	gezelter	938	real(kind=dp) :: U,U1
267	gezelter	1285	real( kind = dp ), dimension(nlocal) :: particle_pot
268	gezelter	115	integer :: atype1
269	gezelter	938	integer :: me, atid1
270	gezelter	115
271	chuckv	1388	! cleanme = .true.
272	gezelter	507	!! Calculate F(rho) and derivative
273	gezelter	115	do atom = 1, nlocal
274	chuckv	592	atid1 = atid(atom)
275	jmichalk	1604	me = eamList%atidToEAMtype(atid1)
276
277			! Myid is set to -1 for non EAM atoms.
278	gezelter	1578	! Punt if we are a non-EAM atom type.
279			if (me == -1) then
280			frho(atom) = 0.0_dp
281			dfrhodrho(atom) = 0.0_dp
282	jmichalk	1604	else
283
284			call lookupEAMSpline1d(EAMList%EAMParams(me)%F, rho(atom), &
285			u, u1)
286	gezelter	938
287	gezelter	1578	frho(atom) = u
288			dfrhodrho(atom) = u1
289			endif
290	gezelter	115	pot = pot + u
291	gezelter	1285	particle_pot(atom) = particle_pot(atom) + u
292	gezelter	939
293	gezelter	115	enddo
294
295			end subroutine calc_eam_preforce_Frho
296	gezelter	938
297	gezelter	115	!! Does EAM pairwise Force calculation.
298	gezelter	1464	subroutine do_eam_pair(atid1, atid2, d, rij, r2, sw, vpair, &
299	gezelter	1419	fpair, pot, f1, rho_i, rho_j, dfrhodrho_i, dfrhodrho_j, &
300	gezelter	1464	fshift_i, fshift_j)
301	gezelter	115	!Arguments
302	gezelter	1464	integer, intent(in) :: atid1, atid2
303	gezelter	115	real( kind = dp ), intent(in) :: rij, r2
304			real( kind = dp ) :: pot, sw, vpair
305	gezelter	1386	real( kind = dp ), dimension(3) :: f1
306	gezelter	115	real( kind = dp ), intent(in), dimension(3) :: d
307			real( kind = dp ), intent(inout), dimension(3) :: fpair
308	chuckv	1388	real( kind = dp ), intent(inout) :: dfrhodrho_i, dfrhodrho_j
309			real( kind = dp ), intent(inout) :: rho_i, rho_j
310			real( kind = dp ), intent(inout):: fshift_i, fshift_j
311	gezelter	115
312	gezelter	938	real( kind = dp ) :: drdx, drdy, drdz
313			real( kind = dp ) :: phab, pha, dvpdr
314			real( kind = dp ) :: rha, drha, dpha
315			real( kind = dp ) :: rhb, drhb, dphb
316	gezelter	115	real( kind = dp ) :: dudr
317	gezelter	938	real( kind = dp ) :: rci, rcj
318			real( kind = dp ) :: drhoidr, drhojdr
319			real( kind = dp ) :: Fx, Fy, Fz
320			real( kind = dp ) :: r, phb
321	chuckv	1388	real( kind = dp ) :: u1, u2
322	gezelter	115
323	gezelter	938	integer :: id1, id2
324	gezelter	115	integer :: mytype_atom1
325			integer :: mytype_atom2
326
327			phab = 0.0E0_DP
328			dvpdr = 0.0E0_DP
329	gezelter	1574	drha = 0.0
330			drhb = 0.0
331	gezelter	115
332			if (rij .lt. EAM_rcut) then
333
334	chuckv	592	mytype_atom1 = EAMList%atidToEAMType(atid1)
335			mytype_atom2 = EAMList%atidTOEAMType(atid2)
336
337
338	gezelter	115	! get cutoff for atom 1
339			rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
340			! get type specific cutoff for atom 2
341			rcj = EAMList%EAMParams(mytype_atom2)%eam_rcut
342	gezelter	507
343	gezelter	115	drdx = d(1)/rij
344			drdy = d(2)/rij
345			drdz = d(3)/rij
346	gezelter	507
347	gezelter	115	if (rij.lt.rci) then
348	gezelter	938
349			! Calculate rho and drho for atom1
350
351	chrisfen	943	call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%rho, &
352	gezelter	938	rij, rha, drha)
353
354			! Calculate Phi(r) for atom1.
355
356	chrisfen	943	call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%phi, &
357	gezelter	938	rij, pha, dpha)
358
359	gezelter	115	endif
360
361			if (rij.lt.rcj) then
362	gezelter	507
363	gezelter	938	! Calculate rho and drho for atom2
364
365	chrisfen	943	call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%rho, &
366	gezelter	938	rij, rhb, drhb)
367
368			! Calculate Phi(r) for atom2.
369
370	chrisfen	943	call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%phi, &
371	gezelter	938	rij, phb, dphb)
372
373	gezelter	115	endif
374
375			if (rij.lt.rci) then
376			phab = phab + 0.5E0_DP(rhb/rha)pha
377			dvpdr = dvpdr + 0.5E0_DP((rhb/rha)dpha + &
378			pha((drhb/rha) - (rhbdrha/rha/rha)))
379			endif
380	gezelter	507
381	gezelter	115	if (rij.lt.rcj) then
382			phab = phab + 0.5E0_DP(rha/rhb)phb
383			dvpdr = dvpdr + 0.5E0_DP((rha/rhb)dphb + &
384			phb((drha/rhb) - (rhadrhb/rhb/rhb)))
385			endif
386	gezelter	507
387	gezelter	115	drhoidr = drha
388			drhojdr = drhb
389
390	gezelter	1574	dudr = drhojdrdfrhodrho_i + drhoidrdfrhodrho_j + dvpdr
391	gezelter	115
392			fx = dudr * drdx
393			fy = dudr * drdy
394			fz = dudr * drdz
395
396	gezelter	1464	! particle_pot is the difference between the full potential
397			! and the full potential without the presence of a particular
398			! particle (atom1).
399			!
400			! This reduces the density at other particle locations, so
401			! we need to recompute the density at atom2 assuming atom1
402			! didn't contribute. This then requires recomputing the
403			! density functional for atom2 as well.
404			!
405			! Most of the particle_pot heavy lifting comes from the
406			! pair interaction, and will be handled by vpair.
407
408			call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom1)%F, &
409			rho_i-rhb, fshift_i, u1)
410			call lookupEAMSpline1d(EAMList%EAMParams(mytype_atom2)%F, &
411			rho_j-rha, fshift_j, u2)
412	gezelter	115
413	gezelter	1386	pot = pot + phab
414	gezelter	938
415	gezelter	1386	f1(1) = f1(1) + fx
416			f1(2) = f1(2) + fy
417			f1(3) = f1(3) + fz
418	gezelter	507
419	gezelter	1386	vpair = vpair + phab
420	gezelter	507
421			endif
422	gezelter	115	end subroutine do_eam_pair
423
424	chrisfen	943	subroutine lookupEAMSpline(cs, xval, yval)
425
426			implicit none
427
428			type (cubicSpline), intent(in) :: cs
429			real( kind = DP ), intent(in) :: xval
430			real( kind = DP ), intent(out) :: yval
431			real( kind = DP ) :: dx
432			integer :: i, j
433			!
434			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
435			! or is nearest to xval.
436
437	gezelter	960	j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
438	chrisfen	943
439			dx = xval - cs%x(j)
440			yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
441
442			return
443			end subroutine lookupEAMSpline
444
445			subroutine lookupEAMSpline1d(cs, xval, yval, dydx)
446
447			implicit none
448
449			type (cubicSpline), intent(in) :: cs
450			real( kind = DP ), intent(in) :: xval
451			real( kind = DP ), intent(out) :: yval, dydx
452			real( kind = DP ) :: dx
453			integer :: i, j
454
455			! Find the interval J = [ cs%x(J), cs%x(J+1) ] that contains
456			! or is nearest to xval.
457
458
459	gezelter	960	j = MAX(1, MIN(cs%n-1, int((xval-cs%x(1)) * cs%dx_i) + 1))
460	chrisfen	943
461			dx = xval - cs%x(j)
462			yval = cs%y(j) + dx(cs%b(j) + dx(cs%c(j) + dx*cs%d(j)))
463
464			dydx = cs%b(j) + dx(2.0d0 cs%c(j) + 3.0d0 * dx * cs%d(j))
465
466			return
467			end subroutine lookupEAMSpline1d
468
469	gezelter	115	end module eam