OOPSE/libmdtools/calc_eam.F90

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


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

  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.
  logical :: nmflag  = .false.

 
  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 IS_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
  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_col
  real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_row
#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 :: setCutoffEAM
  public :: do_eam_pair
  public :: newEAMtype
  public :: calc_eam_prepair_rho
  public :: calc_eam_preforce_Frho
  

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

    write(*,*) "Adding new eamtype: ",eam_ident
    !! 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_dr       = eam_dr
    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

    if (EAMList%currentAddition == 0) then
       call handleError("init_EAM_FF","No members in EAMList")
       status = -1
       return
    end if


       do i = 1, EAMList%currentAddition

! Build array of r values

          do j = 1,EAMList%EAMParams(i)%eam_nr
             EAMList%EAMParams(i)%eam_rvals(j) = &
                  real(j-1,kind=dp)* &
                  EAMList%EAMParams(i)%eam_dr
              end do
! Build array of rho values
          do j = 1,EAMList%EAMParams(i)%eam_nrho
             EAMList%EAMParams(i)%eam_rhovals(j) = &
                  real(j-1,kind=dp)* &
                  EAMList%EAMParams(i)%eam_drho
          end do
          ! 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(1) = 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 for any eam atype
       do i = 2, EAMList%currentAddition 
          current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut)
       end do

       EAM_rcut = current_rcut_max
       EAM_rcut_orig = 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 IS_MPI
    integer :: nrow
    integer :: ncol
#endif
    integer :: alloc_stat


    nlocal = getNlocal()

#ifdef IS_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

    if (allocated(d2frhodrhodrho)) deallocate(d2frhodrhodrho)
    allocate(d2frhodrhodrho(nlocal),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
    
#ifdef IS_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
    if (allocated(d2frhodrhodrho_row)) deallocate(d2frhodrhodrho_row)
    allocate(d2frhodrhodrho_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
    if (allocated(d2frhodrhodrho_col)) deallocate(d2frhodrhodrho_col)
    allocate(d2frhodrhodrho_col(ncol),stat=alloc_stat)
    if (alloc_stat /= 0) then 
       status = -1
       return
    end if
  
#endif

  end subroutine allocateEAM

  subroutine setCutoffEAM(rcut, status)
    real(kind=dp) :: rcut
    integer :: status

    if (rcut < EAM_rcut) then
       EAM_rcut = rcut
    endif


  end subroutine setCutoffEAM


  subroutine clean_EAM()

! clean non-IS_MPI first
    frho = 0.0_dp
    rho  = 0.0_dp
    dfrhodrho = 0.0_dp
! clean MPI if needed
#ifdef IS_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)               :: 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
   
    integer :: myid_atom1
    integer :: myid_atom2

! check to see if we need to be cleaned at the start of a force loop
    if (cleanme) call clean_EAM
    cleanme = .false.
    

#ifdef IS_MPI
    myid_atom1 = atid_Row(atom1)
    myid_atom2 = atid_Col(atom2)
#else
    myid_atom1 = atid(atom1)
    myid_atom2 = atid(atom2)
#endif

    if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then
 

       call eam_splint(EAMList%EAMParams(myid_atom1)%eam_nr, &
            EAMList%EAMParams(myid_atom1)%eam_rvals, &
            EAMList%EAMParams(myid_atom1)%eam_rho_r, &
            EAMList%EAMParams(myid_atom1)%eam_rho_r_pp, &
            r, rho_i_at_j,drho,d2rho)


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

       if (r.lt.EAMList%EAMParams(myid_atom2)%eam_rcut) then
          call eam_splint(EAMList%EAMParams(myid_atom2)%eam_nr, &
               EAMList%EAMParams(myid_atom2)%eam_rvals, &
               EAMList%EAMParams(myid_atom2)%eam_rho_r, &
               EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, &
               r, rho_j_at_i,drho,d2rho)


#ifdef  IS_MPI
          rho_row(atom1) = rho_row(atom1) + rho_j_at_i
#else
          rho(atom1) = rho(atom1) + 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)
    integer :: nlocal
    real(kind=dp) :: pot
    integer :: i,j
    integer :: atom
    real(kind=dp) :: U,U1,U2
    integer :: atype1
    integer :: me
    integer :: n_rho_points
    ! reset clean forces to be true at top of calc rho.
    cleanme = .true.

!! Scatter the electron density from  pre-pair calculation back to local atoms
#ifdef IS_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


!! Calculate F(rho) and derivative 
    do atom = 1, nlocal
       me = atid(atom)
       n_rho_points = EAMList%EAMParams(me)%eam_nrho
       !  Check to see that the density is not greater than the larges rho we have calculated
       if (rho(atom) < EAMList%EAMParams(me)%eam_rhovals(n_rho_points)) then
          call eam_splint(n_rho_points, &
               EAMList%EAMParams(me)%eam_rhovals, &
               EAMList%EAMParams(me)%eam_f_rho, &
               EAMList%EAMParams(me)%eam_f_rho_pp, &
               rho(atom), & ! Actual Rho
               u, u1, u2)
       else 
          ! Calculate F(rho with the largest available rho value
          call eam_splint(n_rho_points, &
               EAMList%EAMParams(me)%eam_rhovals, &
               EAMList%EAMParams(me)%eam_f_rho, &
               EAMList%EAMParams(me)%eam_f_rho_pp, &
               EAMList%EAMParams(me)%eam_rhovals(n_rho_points), & ! Largest rho
               u,u1,u2)
       end if


       frho(i) = u
       dfrhodrho(i) = u1
       d2frhodrhodrho(i) = u2
       pot = pot + u
    enddo


#ifdef IS_MPI
    !! communicate f(rho) and derivatives back into row and column arrays
    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_preforce_Frho


  !! Does EAM pairwise Force calculation.  
  subroutine do_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 ) :: d2
    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  :: mytype_atom1
    integer  :: mytype_atom2


!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
       mytype_atom1 = atid_row(atom1)
#else
       mytype_atom1 = atid(atom1)
#endif
       
       drdx = d(1)/rij
       drdy = d(2)/rij
       drdz = d(3)/rij
       

       call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
               EAMList%EAMParams(mytype_atom1)%eam_rho_r, &
               EAMList%EAMParams(mytype_atom1)%eam_rho_r_pp, &
               rij, rha,drha,d2rha)

       !! Calculate Phi(r) for atom1.
       call eam_splint(EAMList%EAMParams(mytype_atom1)%eam_nr, &
               EAMList%EAMParams(mytype_atom1)%eam_rvals, &
               EAMList%EAMParams(mytype_atom1)%eam_phi_r, &
               EAMList%EAMParams(mytype_atom1)%eam_phi_r_pp, &
               rij, pha,dpha,d2pha)


! get cutoff for atom 1
       rci = EAMList%EAMParams(mytype_atom1)%eam_rcut
#ifdef IS_MPI
       mytype_atom2 = atid_col(atom2)
#else
       mytype_atom2 = atid(atom2)
#endif

       ! Calculate rho,drho and d2rho for atom1
       call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
               EAMList%EAMParams(mytype_atom2)%eam_rho_r, &
               EAMList%EAMParams(mytype_atom2)%eam_rho_r_pp, &
               rij, rhb,drhb,d2rhb)

       !! Calculate Phi(r) for atom2.
       call eam_splint(EAMList%EAMParams(mytype_atom2)%eam_nr, &
               EAMList%EAMParams(mytype_atom2)%eam_rvals, &
               EAMList%EAMParams(mytype_atom2)%eam_phi_r, &
               EAMList%EAMParams(mytype_atom2)%eam_phi_r_pp, &
               rij, phb,dphb,d2phb)


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


       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)))
          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 (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)))
          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 IS_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 IS_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 (nmflag) then

          drhoidr = drha
          drhojdr = drhb
          d2rhoidrdr = d2rha
          d2rhojdrdr = d2rhb

#ifdef IS_MPI
          d2 = d2vpdrdr + &
               d2rhoidrdr*dfrhodrho_col(atom2) + &
               d2rhojdrdr*dfrhodrho_row(atom1) + &
               drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + &
               drhojdr*drhojdr*d2frhodrhodrho_row(atom1)
               
#else

          d2 = d2vpdrdr + &
               d2rhoidrdr*dfrhodrho(atom2) + &
               d2rhojdrdr*dfrhodrho(atom1) + &
               drhoidr*drhoidr*d2frhodrhodrho(atom2) + &
               drhojdr*drhojdr*d2frhodrhodrho(atom1)
#endif
       end if


       if (do_stress) then
          
#ifdef IS_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) - d(1) * fx
             tau_Temp(2) = tau_Temp(2) - d(1) * fy
             tau_Temp(3) = tau_Temp(3) - d(1) * fz
             tau_Temp(4) = tau_Temp(4) - d(2) * fx
             tau_Temp(5) = tau_Temp(5) - d(2) * fy
             tau_Temp(6) = tau_Temp(6) - d(2) * fz
             tau_Temp(7) = tau_Temp(7) - d(3) * fx
             tau_Temp(8) = tau_Temp(8) - d(3) * fy
             tau_Temp(9) = tau_Temp(9) - d(3) * fz

             virial_Temp = virial_Temp + &
                  (tau_Temp(1) + tau_Temp(5) + tau_Temp(9))

          endif
       endif   
    endif

    
  end subroutine do_eam_pair


  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
    real( kind = DP ) :: dy, d2y
    real( kind = DP ) :: del, h, a, b, c, d
    integer :: pp_arraySize

 
    ! 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(real((nx-1),kind=dp) * (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(errMSG,*) "EAM_splint: x is outside bounds of spline"
       call handleError(routineName,errMSG)
    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(errMSG,*) "EAM_splint:",x," x is outside bounding range"
       call handleError(routineName,errMSG)
       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(len=*) :: boundary
    
    ! make sure the sizes match
    if ((nx /= size(xa)) .or. (nx /= size(ya))) then
       call handleWarning("EAM_SPLINE","Array size mismatch")
    end if

    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:	657
Committed:	Wed Jul 30 21:17:01 2003 UTC (20 years, 11 months ago) by chuckv
File size:	29161 byte(s)
Log Message:	More bug fixes for eam.