OOPSE/libmdtools/calc_eam.F90

module eam
  use definitions, ONLY : DP,default_error
  use simulation
  use force_globals
  use status
  use atype_module
#ifdef MPI
  use mpiSimulation
#endif
  implicit none
  PRIVATE


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

  character(len = 200) :: errMsg
  character(len=*), parameter :: RoutineName =  "EAM MODULE"
  logical :: cleanme = .true.


  type, private :: EAMtype
     integer           :: eam_atype       
     real( kind = DP ) :: eam_dr          
     integer           :: eam_nr           
     integer           :: eam_nrho          
     real( kind = DP ) :: eam_lattice        
     real( kind = DP ) :: eam_drho      
     real( kind = DP ) :: eam_rcut     
     integer           :: eam_atype_map
     
     real( kind = DP ), pointer, dimension(:) :: eam_rvals        => null()
     real( kind = DP ), pointer, dimension(:) :: eam_rhovals      => null()
     real( kind = DP ), pointer, dimension(:) :: eam_F_rho        => null()
     real( kind = DP ), pointer, dimension(:) :: eam_Z_r          => null()
     real( kind = DP ), pointer, dimension(:) :: eam_rho_r        => null()
     real( kind = DP ), pointer, dimension(:) :: eam_phi_r        => null()
     real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp     => null()
     real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp       => null()
     real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp     => null()
     real( kind = DP ), pointer, dimension(:) :: eam_phi_r_pp     => null()
  end type EAMtype


  !! Arrays for derivatives used in force calculation
  real( kind = dp), dimension(:), allocatable :: frho
  real( kind = dp), dimension(:), allocatable :: rho

  real( kind = dp), dimension(:), allocatable :: dfrhodrho
!  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho


!! Arrays for MPI storage
#ifdef MPI
  real( kind = dp), dimension(:), allocatable :: dfrhodrho_col
  real( kind = dp), dimension(:), allocatable :: dfrhodrho_row
  real( kind = dp), dimension(:), allocatable :: frho_row
  real( kind = dp), dimension(:), allocatable :: frho_col
  real( kind = dp), dimension(:), allocatable :: rho_row
  real( kind = dp), dimension(:), allocatable :: rho_col

#endif

  type, private :: EAMTypeList
     integer           :: n_eam_types = 0
     integer           :: currentAddition = 0
     
     type (EAMtype), pointer  :: EAMParams(:) => null()
  end type EAMTypeList


  type (eamTypeList) :: EAMList

  !! standard eam stuff  


  public :: init_EAM_FF
!  public :: EAM_new_rcut
!  public :: do_EAM_pair
  public :: newEAMtype

  
contains


  subroutine newEAMtype(lattice_constant,eam_nrho,eam_drho,eam_nr,&
       eam_dr,rcut,eam_Z_r,eam_rho_r,eam_F_rho,&
       eam_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
    real (kind = dp ), dimension(eam_nr)   :: eam_rho_r
    real (kind = dp ), dimension(eam_nrho) :: eam_F_rho
    integer                                :: eam_ident
    integer                                :: status

    integer                                :: nAtypes
    integer                                :: maxVals
    integer                                :: alloc_stat
    integer                                :: current
    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., nAtypes, MatchList)
       EAMList%n_eam_types = nAtypes
       allocate(EAMList%EAMParams(nAtypes))
    end if

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

    !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat)
    if (alloc_stat /= 0) then
       status = -1
       return
    end if


    ! this is a possible bug, we assume a correspondence between the vector atypes and
    ! EAMAtypes
      
    EAMList%EAMParams(current)%eam_atype    = eam_ident
    EAMList%EAMParams(current)%eam_lattice  = lattice_constant
    EAMList%EAMParams(current)%eam_nrho     = eam_nrho
    EAMList%EAMParams(current)%eam_drho     = eam_drho
    EAMList%EAMParams(current)%eam_nr       = eam_nr
    EAMList%EAMParams(current)%eam_rcut     = rcut
    EAMList%EAMParams(current)%eam_Z_r      = eam_Z_r
    EAMList%EAMParams(current)%eam_rho_r    = eam_rho_r
    EAMList%EAMParams(current)%eam_F_rho    = eam_F_rho

  end subroutine newEAMtype


  subroutine init_EAM_FF(status)
    integer :: status
    integer :: i,j
    real(kind=dp) :: current_rcut_max
    integer :: alloc_stat
    integer :: number_r, number_rho

       do i = 1, EAMList%currentAddition
  
          EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = &
               real(EAMList%EAMParams(i)%eam_nr,kind=dp)* &
               EAMList%EAMParams(i)%eam_dr

          EAMList%EAMParams(i)%eam_rhovals(1:EAMList%EAMParams(i)%eam_nrho) = &
               real(EAMList%EAMParams(i)%eam_nrho,kind=dp)* &
               EAMList%EAMParams(i)%eam_drho

          ! convert from eV to kcal / mol:
          EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP

          ! precompute the pair potential and get it into kcal / mol:
          EAMList%EAMParams(i)%eam_phi_r = 0.0E0_DP
          do j = 2, EAMList%EAMParams(i)%eam_nr
             EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j)
             EAMList%EAMParams(i)%eam_phi_r(j) =  EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP
          enddo
       end do

       do i = 1,  EAMList%currentAddition
          number_r   = EAMList%EAMParams(i)%eam_nr
          number_rho = EAMList%EAMParams(i)%eam_nrho
          
          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
               EAMList%EAMParams(i)%eam_rho_r, &
               EAMList%EAMParams(i)%eam_rho_r_pp, &
               0.0E0_DP, 0.0E0_DP, 'N')
          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
               EAMList%EAMParams(i)%eam_Z_r, &
               EAMList%EAMParams(i)%eam_Z_r_pp, &
               0.0E0_DP, 0.0E0_DP, 'N')
          call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, &
               EAMList%EAMParams(i)%eam_F_rho, &
               EAMList%EAMParams(i)%eam_F_rho_pp, &
               0.0E0_DP, 0.0E0_DP, 'N')
          call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, &
               EAMList%EAMParams(i)%eam_phi_r, &
               EAMList%EAMParams(i)%eam_phi_r_pp, &
               0.0E0_DP, 0.0E0_DP, 'N')
       enddo
       
       current_rcut_max = EAMList%EAMParams(1)%eam_rcut
       !! find the smallest rcut
       do i = 2, EAMList%currentAddition 
          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
       end do

       EAM_rcut = current_rcut_max
!       do i = 1, EAMList%currentAddition
!          EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i
!       end do


       !! Allocate arrays for force calculation
          call allocateEAM(alloc_stat)
          if (alloc_stat /= 0 ) then
             status = -1
             return
          endif

  end subroutine init_EAM_FF

!! routine checks to see if array is allocated, deallocates array if allocated
!! and then creates the array to the required size
  subroutine allocateEAM(status)
    integer, intent(out) :: status

    integer :: nlocal
#ifdef MPI
    integer :: nrow
    integer :: ncol
#endif
    integer :: alloc_stat


    nlocal = getNlocal()

#ifdef MPI
    nrow = getNrow(plan_row)
    ncol = getNcol(plan_col)
#endif

    if (allocated(frho)) deallocate(frho)
    allocate(frho(nlocal),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    if (allocated(rho)) deallocate(rho)
    allocate(rho(nlocal),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    if (allocated(dfrhodrho)) deallocate(dfrhodrho)
    allocate(dfrhodrho(nlocal),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    
#ifdef MPI

    if (allocated(frho_row)) deallocate(frho_row)
    allocate(frho_row(nrow),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    if (allocated(rho_row)) deallocate(rho_row)
    allocate(rho_row(nrow),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    if (allocated(dfrhodrho_row)) deallocate(dfrhodrho_row)
    allocate(dfrhodrho_row(nrow),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if

! Now do column arrays

    if (allocated(frho_col)) deallocate(frho_col)
    allocate(frho_col(ncol),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    if (allocated(rho_col)) deallocate(rho_col)
    allocate(rho_col(ncol),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    if (allocated(dfrhodrho_col)) deallocate(dfrhodrho_col)
    allocate(dfrhodrho_col(ncol),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if

#endif

  end subroutine allocateEAM


  subroutine clean_EAM()

! clean non-MPI first
    frho = 0.0_dp
    rho  = 0.0_dp
    dfrhodrho = 0.0_dp
! clean MPI if needed
#ifdef MPI
    frho_row = 0.0_dp
    frho_col = 0.0_dp
    rho_row  = 0.0_dp
    rho_col  = 0.0_dp
    dfrhodrho_row = 0.0_dp
    dfrhodrho_col = 0.0_dp
#endif
  end subroutine clean_EAM


  subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat)
    integer, intent(in)          :: eam_n_rho
    integer, intent(in)          :: eam_n_r
    type (EAMType), pointer      :: thisEAMType
    integer, optional   :: stat
    integer             :: alloc_stat


    if (present(stat)) stat = 0
    
    allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat)   
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat)   
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat)   
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat)        
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat)      
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat)      
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat)   
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat)   
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat)   
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
    allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat)
    if (alloc_stat /= 0 ) then
       if (present(stat)) stat = -1
       return
    end if
       

  end subroutine allocate_EAMType


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

    ! free Arrays in reverse order of allocation...
    deallocate(thisEAMType%eam_phi_r_pp)      
    deallocate(thisEAMType%eam_rho_r_pp)   
    deallocate(thisEAMType%eam_Z_r_pp)   
    deallocate(thisEAMType%eam_F_rho_pp)   
    deallocate(thisEAMType%eam_phi_r)      
    deallocate(thisEAMType%eam_rho_r)      
    deallocate(thisEAMType%eam_Z_r)   
    deallocate(thisEAMType%eam_F_rho)
    deallocate(thisEAMType%eam_rhovals)
    deallocate(thisEAMType%eam_rvals)
   
  end subroutine deallocate_EAMType

!! Calculates rho_r
  subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq)
    integer :: atom1,atom2
    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) :: rho_i_at_j
    ! value of electron density rho do to atom j at atom i
    real(kind = dp) :: rho_j_at_i

    ! we don't use the derivatives, dummy variables
    real( kind = dp) :: drho,d2rho
    integer :: eam_err
   
    if (cleanme) call clean_EAM
    cleanme = .false.

    call  calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1)

#ifdef  MPI
    rho_col(atom2) = rho_col(atom2) + rho_i_at_j
#else
    rho(atom2) = rho(atom2) + rho_i_at_j
#endif

    call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2)

#ifdef  MPI
    rho_row(atom1) = rho_row(atom1) + rho_j_at_i
#else
    rho(atom1) = rho(atom1) + rho_j_at_i
#endif
  end subroutine calc_eam_prepair_rho

  !! Calculate the functional F(rho) for all atoms
  subroutine calc_eam_prepair_Frho(nlocal,pot)
    integer :: nlocal
    real(kind=dp) :: pot
    integer :: i,j
    real(kind=dp) :: U,U1,U2
    integer :: atype1
    ! reset clean forces to be true at top of calc rho.
    cleanme = .true.

!! Scatter the electron density in pre-pair
#ifdef MPI
    call scatter(rho_row,rho,plan_row,eam_err)
    if (eam_err /= 0 ) then
      write(errMsg,*) " Error scattering rho_row into rho"
      call handleError(RoutineName,errMesg)
   endif       
    call scatter(rho_col,rho,plan_col,eam_err)
    if (eam_err /= 0 ) then
      write(errMsg,*) " Error scattering rho_col into rho"
      call handleError(RoutineName,errMesg)
   endif
#endif

   do i = 1, nlocal
      call calc_eam_frho(rho(i),u,u1,u2,atype1)
      frho(i) = u
      dfrhodrho(i) = u1
!      d2frhodrhodrho(i) = u2
      pot = pot + u
   enddo

#ifdef MPI
    !! communicate f(rho) and derivatives
    call gather(frho,frho_row,plan_row, eam_err)
    if (eam_err /=  0) then
       call handleError("cal_eam_forces()","MPI gather frho_row failure")
    endif
    call gather(dfrhodrho,dfrhodrho_row,plan_row, eam_err)
    if (eam_err /=  0) then
       call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure")
    endif
    call gather(frho,frho_col,plan_col, eam_err)
    if (eam_err /=  0) then
       call handleError("cal_eam_forces()","MPI gather frho_col failure")
    endif
    call gather(dfrhodrho,dfrhodrho_col,plan_col, eam_err)
    if (eam_err /=  0) then
       call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure")
    endif

    if (nmflag) then
       call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row)
       call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col)
    endif
#endif


  end subroutine calc_eam_prepair_Frho


  subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress)
!Arguments    
    integer, intent(in) ::  atom1, atom2
    real( kind = dp ), intent(in) :: rij, r2
    real( kind = dp ) :: pot
    real( kind = dp ), dimension(3,getNlocal()) :: f
    real( kind = dp ), intent(in), dimension(3) :: d
    logical, intent(in) :: do_pot, do_stress

    real( kind = dp ) :: drdx,drdy,drdz
    real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr
    real( kind = dp ) :: rha,drha,d2rha, dpha
    real( kind = dp ) :: rhb,drhb,d2rhb, dphb
    real( kind = dp ) :: dudr
    real( kind = dp ) :: rci,rcj
    real( kind = dp ) :: drhoidr,drhojdr
    real( kind = dp ) :: d2rhoidrdr
    real( kind = dp ) :: d2rhojdrdr
    real( kind = dp ) :: Fx,Fy,Fz
    real( kind = dp ) :: r,d2pha,phb,d2phb
    integer :: id1,id2
    integer  :: atype1,atype2


!Local Variables
    

    phab = 0.0E0_DP
    dvpdr = 0.0E0_DP
    d2vpdrdr = 0.0E0_DP
       

    if (rij .lt. EAM_rcut) then
#ifdef IS_MPI
!!!!! FIX ME
       atype1 = atid_row(atom1)
#else
       atype1 = atid(atom1)
#endif

       drdx = d(1)/rij
       drdy = d(2)/rij
       drdz = d(3)/rij
       
       
       call calc_eam_rho(r, rha, drha, d2rha, atype1)
       call calc_eam_phi(r, pha, dpha, d2pha, atype1)
  !     rci = eam_rcut(eam_atype_map(atom1))
#ifdef IS_MPI
       atype2 = atid_col(atom2)
#else
       atype2 = atid(atom2)
#endif
       
       call calc_eam_rho(r, rhb, drhb, d2rhb, atype2)
       call calc_eam_phi(r, phb, dphb, d2phb, atype2)
  !     rcj = eam_rcut(eam_atype_map(atype2))

       if (r.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)))
          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rhb/rha)*d2pha + &
               2.0E0_DP*dpha*((drhb/rha) - (rhb*drha/rha/rha)) + &
               pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + &
               (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha)))
       endif
       

       if (r.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)))
          d2vpdrdr = d2vpdrdr + 0.5E0_DP*((rha/rhb)*d2phb + &
               2.0E0_DP*dphb*((drha/rhb) - (rha*drhb/rhb/rhb)) + &
               phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + &
               (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb)))
       endif
       
       drhoidr = drha
       drhojdr = drhb

       d2rhoidrdr = d2rha
       d2rhojdrdr = d2rhb


#ifdef MPI
       dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) &
            + dvpdr
    
#else
       dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) &
            + dvpdr

#endif

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


#ifdef MPI
       if (do_pot) then
          pot_Row(atom1) = pot_Row(atom1) + phab*0.5
          pot_Col(atom2) = pot_Col(atom2) + phab*0.5
       end if

       f_Row(1,atom1) = f_Row(1,atom1) + fx
       f_Row(2,atom1) = f_Row(2,atom1) + fy
       f_Row(3,atom1) = f_Row(3,atom1) + fz

       f_Col(1,atom2) = f_Col(1,atom2) - fx
       f_Col(2,atom2) = f_Col(2,atom2) - fy
       f_Col(3,atom2) = f_Col(3,atom2) - fz
#else
       if(do_pot) pot = pot + phab

       f(1,atom1) = f(1,atom1) + fx
       f(2,atom1) = f(2,atom1) + fy
       f(3,atom1) = f(3,atom1) + fz

       f(1,atom2) = f(1,atom2) - fx
       f(2,atom2) = f(2,atom2) - fy
       f(3,atom2) = f(3,atom2) - fz
#endif


       if (do_stress) then

#ifdef MPI
          id1 = tagRow(atom1)
          id2 = tagColumn(atom2)
#else
          id1 = atom1
          id2 = atom2
#endif

          if (molMembershipList(id1) .ne. molMembershipList(id2)) then
             
             tau_Temp(1) = tau_Temp(1) + fx * d(1)
             tau_Temp(2) = tau_Temp(2) + fx * d(2)
             tau_Temp(3) = tau_Temp(3) + fx * d(3)
             tau_Temp(4) = tau_Temp(4) + fy * d(1)
             tau_Temp(5) = tau_Temp(5) + fy * d(2)
             tau_Temp(6) = tau_Temp(6) + fy * d(3)
             tau_Temp(7) = tau_Temp(7) + fz * d(1)
             tau_Temp(8) = tau_Temp(8) + fz * d(2)
             tau_Temp(9) = tau_Temp(9) + fz * d(3)
             virial_Temp = virial_Temp + &
                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))

          endif
       endif
       
    endif


end subroutine calc_eam_pair

!!$subroutine calc_eam_rho(r, rho, drho, d2rho, atype)
!!$
!!$  !  include 'headers/sizes.h'
!!$
!!$
!!$integer atype, etype, number_r
!!$real( kind = DP )  :: r, rho, drho, d2rho
!!$integer :: i
!!$
!!$
!!$etype = eam_atype_map(atype)
!!$
!!$if (r.lt.eam_rcut(etype)) then
!!$number_r = eam_nr(etype)
!!$call eam_splint(etype, number_r, eam_rvals, eam_rho_r, &
!!$   eam_rho_r_pp, r, rho, drho, d2rho)
!!$else
!!$rho = 0.0E0_DP
!!$drho = 0.0E0_DP
!!$d2rho = 0.0E0_DP
!!$endif
!!$
!!$return
!!$end subroutine calc_eam_rho
!!$
!!$subroutine calc_eam_frho(dens, u, u1, u2, atype)
!!$
!!$  ! include 'headers/sizes.h'
!!$
!!$integer atype, etype, number_rho
!!$real( kind = DP ) :: dens, u, u1, u2
!!$real( kind = DP ) :: rho_vals
!!$
!!$etype = eam_atype_map(atype)
!!$number_rho = eam_nrho(etype)
!!$if (dens.lt.eam_rhovals(number_rho, etype)) then
!!$call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
!!$   eam_f_rho_pp, dens, u, u1, u2)
!!$else
!!$rho_vals = eam_rhovals(number_rho,etype)
!!$call eam_splint(etype, number_rho, eam_rhovals, eam_f_rho, &
!!$   eam_f_rho_pp, rho_vals, u, u1, u2)
!!$endif
!!$
!!$return
!!$end subroutine calc_eam_frho
!!$
!!$subroutine calc_eam_phi(r, phi, dphi, d2phi, atype)
!!$
!!$
!!$
!!$
!!$integer atype, etype, number_r
!!$real( kind = DP ) :: r, phi, dphi, d2phi
!!$
!!$etype = eam_atype_map(atype)
!!$
!!$if (r.lt.eam_rcut(etype)) then
!!$number_r = eam_nr(etype)
!!$call eam_splint(etype, number_r, eam_rvals, eam_phi_r, &
!!$   eam_phi_r_pp, r, phi, dphi, d2phi)
!!$else
!!$phi = 0.0E0_DP
!!$dphi = 0.0E0_DP
!!$d2phi = 0.0E0_DP
!!$endif
!!$
!!$return
!!$end subroutine calc_eam_phi


  subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y)


    integer :: atype, nx, j
    real( kind = DP ), dimension(:) :: xa
    real( kind = DP ), dimension(:) :: ya
    real( kind = DP ), dimension(:) :: yppa
    real( kind = DP ) :: x, y, dy, d2y
    real( kind = DP ) :: del, h, a, b, c, d


    ! this spline code assumes that the x points are equally spaced
    ! do not attempt to use this code if they are not.
    
    
    ! find the closest point with a value below our own:
    j = FLOOR(dble(nx-1) * (x - xa(1)) / (xa(nx) - xa(1))) + 1

    ! check to make sure we're inside the spline range:
    if ((j.gt.nx).or.(j.lt.1)) then
       write(default_error,*) "EAM_splint: x is outside bounds of spline"
    endif
    ! check to make sure we haven't screwed up the calculation of j:
    if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then
       if (j.ne.nx) then
        write(default_error,*) "EAM_splint: x is outside bounding range"
       endif
    endif

    del = xa(j+1) - x
    h = xa(j+1) - xa(j)
    
    a = del / h
    b = 1.0E0_DP - a
    c = a*(a*a - 1.0E0_DP)*h*h/6.0E0_DP
    d = b*(b*b - 1.0E0_DP)*h*h/6.0E0_DP
    
    y = a*ya(j) + b*ya(j+1) + c*yppa(j) + d*yppa(j+1)
    
    dy = (ya(j+1)-ya(j))/h &
         - (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j)/6.0E0_DP &
         + (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1)/6.0E0_DP
    
    d2y = a*yppa(j) + b*yppa(j+1)

  end subroutine eam_splint

  subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary)

  
    ! yp1 and ypn are the first derivatives of y at the two endpoints
    ! if boundary is 'L' the lower derivative is used
    ! if boundary is 'U' the upper derivative is used
    ! if boundary is 'B' then both derivatives are used
    ! if boundary is anything else, then both derivatives are assumed to be 0
    
    integer :: nx, i, k, max_array_size
    
    real( kind = DP ), dimension(:)        :: xa
    real( kind = DP ), dimension(:)        :: ya
    real( kind = DP ), dimension(:)        :: yppa
    real( kind = DP ), dimension(size(xa)) :: u
    real( kind = DP ) :: yp1,ypn,un,qn,sig,p
    character boundary
    
    
    if ((boundary.eq.'l').or.(boundary.eq.'L').or. &
         (boundary.eq.'b').or.(boundary.eq.'B')) then
       yppa(1) = -0.5E0_DP
       u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-&
            ya(1))/(xa(2)-xa(1))-yp1)
    else
       yppa(1) = 0.0E0_DP
       u(1)  = 0.0E0_DP
    endif
    
    do i = 2, nx - 1
       sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1))
       p = sig * yppa(i-1) + 2.0E0_DP
       yppa(i) = (sig - 1.0E0_DP) / p
       u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - &
            (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ &
            (xa(i+1)-xa(i-1)) - sig * u(i-1))/p
    enddo
    
    if ((boundary.eq.'u').or.(boundary.eq.'U').or. &
         (boundary.eq.'b').or.(boundary.eq.'B')) then
       qn = 0.5E0_DP
       un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* &
            (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1)))
    else
       qn = 0.0E0_DP
       un = 0.0E0_DP
    endif

    yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP)
    
    do k = nx-1, 1, -1
       yppa(k)=yppa(k)*yppa(k+1)+u(k)
    enddo

  end subroutine eam_spline


end module eam
Revision:	631
Committed:	Thu Jul 17 19:25:51 2003 UTC (22 years, 3 months ago) by chuckv
Original Path:	*trunk/OOPSE/libmdtools/calc_eam.F90*
File size:	24873 byte(s)
Log Message:	Added massive changes for eam....