module potential
  ! EAM 7
  use algor
  use constants,  only: dp
  use cell,       only: cell_type
  implicit none

  public

  character(len=10),  parameter ::  potential_type  = 'eam'      ! eam or harmonic
  logical,            parameter ::  non_local_correction = .false.
  logical,            save      ::  pot_initialized = .false.

  ! EAM potential:
  ! Variable names are for first use - when no longer required, pointers (correct name) point to these arrays for storage
  ! For both Ni and H:
  real(kind=dp),  dimension(:,:,:), allocatable           ::  rij_vec
  real(kind=dp),  dimension(:,:),   allocatable,  target  ::  rij_mag       ! seperate from rij_vec array for efficiency...
  real(kind=dp),  dimension(:),     allocatable,  target  ::  rho_bar
  real(kind=dp),  dimension(:,:),   allocatable           ::  d_rho_atomic
  real(kind=dp),  dimension(:,:),                 pointer ::  d_phi
  real(kind=dp),  dimension(:),                   pointer ::  d_f_rb

  real(kind=dp),  dimension(:,:),   allocatable,  target  ::  rho_hat
  real(kind=dp),  dimension(:,:),                 pointer ::  d_f_rh
  logical                                                 ::  h_present

  integer,        parameter                     ::  h_num = 1, ni_num = 28  ! faster than lookup in consts

  ! Harmonic potential:
  real(kind=dp)                 ::  harm_period = 100.0_dp                            ! isotropic...
  real(kind=dp),  dimension(3)  ::  harm_spring_const = (/ 1.0_dp, 1.0_dp, 1.0_dp /)  ! can be anisotropic..

  public  ::  pot_get_potential
  public  ::  pot_get_forces

contains

subroutine pot_get_potential(cell, energy)
  use io, only: io_err
  implicit none
  type(cell_type),  intent(in)  ::  cell
  real(kind=dp),    intent(out) ::  energy

  select case (potential_type)
    case ('eam')
      if (pot_initialized .eqv. .false.) call pot_eam_init(cell)
      call pot_eam_potential(cell, energy)
    case default
      call io_err("pot_get_potential: unsupported potential type")
  end select
end subroutine pot_get_potential

subroutine pot_get_forces(cell, force, potential)
  use constants,  only: two_pi, units_natural_to_atomic
  use io,         only: io_err
  implicit none
  type(cell_type),                          intent(in)    ::  cell
  real(kind=dp),  dimension(3,cell%natoms), intent(out)   ::  force
  real(kind=dp),                            intent(out)   ::  potential
  ! local for SHO
  real(kind=dp) ::  omega

  select case (potential_type)
    case ('eam')
      if (pot_initialized .eqv. .false.) call pot_eam_init(cell)
      call pot_eam_forces(cell, force, potential)
    case default
      call io_err("Error in pot_get_force: unsupported potential type")
  end select
end subroutine pot_get_forces

subroutine pot_eam_init(cell)
  use io,   only: io_err
  implicit none
  type(cell_type),  intent(in)  ::  cell
  integer ::  iatom, istat, n_hydrogen

  h_present = .false.
  n_hydrogen = 0

  ! first check that all atoms are either hydrogen or nickel:
  do iatom = 1, cell%natoms
    if ((cell%atom_species(iatom) .ne. h_num) .and. (cell%atom_species(iatom) .ne. ni_num)) then
      call io_err("Error in potential - pot_eam_init: Unexpected atom type.")
    else
      if (cell%atom_species(iatom) .eq. h_num) n_hydrogen = n_hydrogen + 1
    end if
  end do

  if (n_hydrogen .gt. 0) h_present = .true.

  allocate(rij_vec(3, cell%natoms, cell%natoms), stat=istat)!arraysize), stat=istat)
  if (istat .ne. 0) call io_err("Error in pot_eam_init - failed to allocate rij_vec array.")

  allocate(rij_mag(cell%natoms, cell%natoms), stat=istat)!arraysize), stat=istat)
  if (istat .ne. 0) call io_err("Error in pot_eam_init - failed to allocate rij_mag array.")

  allocate(rho_bar(cell%natoms), stat=istat)
  if (istat .ne. 0) call io_err("Error in pot_eam_init - failed to allocate rho_bar array.")

  allocate(d_rho_atomic(cell%natoms, cell%natoms), stat=istat)
  if (istat .ne. 0) call io_err("Error in pot_eam_init - failed to allocate d_rho_atomic array.")

  ! initialize pointers here:
  ! once we are done with data on RHS, we will store LHS data in these arrays instead
  ! pointers save memory and make code easier to follow
  d_phi => rij_mag
  d_f_rb => rho_bar

  allocate(rho_hat(2, max(1,n_hydrogen)), stat=istat)
  if (istat .ne. 0) call io_err("pot_eam_init - could not allocate rho_hat array.")
  d_f_rh => rho_hat

  pot_initialized = .true.
end subroutine pot_eam_init

subroutine pot_eam_potential(cell, potential)
  use constants,  only: units_atomic_to_natural, units_natural_to_atomic
  use cell,       only: cell_min_img
  !use algor,      only: algor_symmetric_matrix_index
  implicit none
  type(cell_type),  intent(in)  ::  cell
  real(kind=dp),    intent(out) ::  potential
  ! local variables
  real(kind=dp) ::  rho_a
  integer ::  iatom, jatom, ihydrogen!, ielement, i

  potential = 0.0_dp  ! we have sums

  ! rij vectors and pairwise terms:
  do iatom = 1, cell%natoms!-1
    do jatom = 1, cell%natoms!iatom+1, cell%natoms
      if (iatom .ne. jatom) then
        !ielement = algor_symmetric_matrix_index(iatom, jatom, cell%natoms)
        ! could instead set ielement = 1 outside loop then do ielement = ielement+1 inside inner loop
        ! followed by the same thing in the outer loop (after the inner loop has completed) to match this.
        ! speed difference is expected to be small, so just call this function to keep things simple.

        ! rij(4, ielement) is magintude of vector
        call cell_min_img(cell, iatom, jatom, rij_vec(:,iatom,jatom), distance=rij_mag(iatom,jatom))!ielement), distance=rij(4,ielement))
        ! make unit vectors but keep original magnitude
        rij_vec(:,iatom,jatom) = rij_vec(:,iatom,jatom)/rij_mag(iatom,jatom)!ielement) = rij(1:3,ielement)/rij(4,ielement)

        ! scale magnitude: length from Bohr to Angstroms
        rij_mag(iatom,jatom) = units_atomic_to_natural(length=rij_mag(iatom,jatom))

        ! Calculate the pairwise potential energy contributions:
        ! Since we are only working on half of the matrix (excluding the diagonal), we don't need to x0.5
        ! Note: no need to store each of these terms as they only contribute to the potential energy
        potential = potential + 0.5_dp*eam_phi(rij_mag(iatom,jatom), cell%atom_species(iatom), cell%atom_species(jatom))
        !eam_phi(rij(4,ielement), cell%atom_species(iatom), cell%atom_species(jatom))
      end if
    end do  !jatom
  end do  !iatom

  ! Get density dependent terms
  rho_bar(:) = 0.0_dp

  ! extra terms for H
  rho_hat(:,:) = 0.0_dp
  ihydrogen = 1

  do iatom = 1, cell%natoms

    do jatom = 1, cell%natoms
      ! i .ne. j
      if (iatom .ne. jatom) then
        ! get index of reduced N^2 array
!        ielement = algor_symmetric_matrix_index(iatom, jatom, cell%natoms)

        rho_a = eam_rho(rij_mag(iatom,jatom), cell%atom_species(jatom))!ielement), cell%atom_species(jatom))
        rho_bar(iatom) = rho_bar(iatom) + rho_a

        if (h_present .and. (cell%atom_species(iatom) .eq. h_num)) then
          rho_hat(1,ihydrogen) = rho_hat(1,ihydrogen) + rho_a**2        ! to power Q1=2
          rho_hat(2,ihydrogen) = rho_hat(2,ihydrogen) + sqrt(rho_a)     ! to power Q2=0.5
        end if

      end if
    end do  !jatom

    ! rho_hat(:,1) = 0.0_dp if hydrogen doesn't exist.. won't affect eam_f_ni anyway
    potential = potential + eam_f(rho_bar(iatom), rho_hat(1,ihydrogen), rho_hat(2,ihydrogen), cell%atom_species(iatom))
    if (h_present .and. (cell%atom_species(iatom) .eq. h_num)) ihydrogen = ihydrogen + 1
  end do !iatom

  ! rescale units into atomic
  potential = units_natural_to_atomic(energy=potential)
end subroutine pot_eam_potential

subroutine pot_eam_forces(cell, force, potential)
  use constants,  only: units_natural_to_atomic
!  use algor,      only: algor_symmetric_matrix_index
  implicit none
  type(cell_type),                          intent(in)  ::  cell
  real(kind=dp),  dimension(3,cell%natoms), intent(out) ::  force
  real(kind=dp),                            intent(out) ::  potential
  ! local variables
  real(kind=dp) :: contrib
  integer ::  iatom, jatom, ihydrogen, jhydrogen, ihat

  ! need to fill out some of the matrices anyway, so potential is basically free...
  call pot_eam_potential(cell, potential)

  ! get everything we need for the forces in one go...
  do iatom = 1, cell%natoms
    ! have rho_bar, need dF/rho_bar         0 is rho_bar derivative
    d_f_rb(iatom) = eam_d_f(rho_bar(iatom), 0, cell%atom_species(iatom))

    do jatom = 1, cell%natoms
      if (iatom .ne. jatom) then
        ! MUST do d_rho_atomic first since d_phi points to the rij_mag array!
        d_rho_atomic(iatom,jatom) = eam_d_rho(rij_mag(iatom,jatom), cell%atom_species(jatom))
        d_phi(iatom, jatom) = eam_d_phi(rij_mag(iatom,jatom), cell%atom_species(iatom), cell%atom_species(jatom))
      end if
    end do
  end do

  if (h_present) then
    do ihydrogen = 1, size(rho_hat,2)
      do ihat = 1, 2
        d_f_rh(ihat,ihydrogen) = eam_d_f_h(rho_hat(ihat,ihydrogen), ihat)   ! always use eam_d_f_h
      end do
    end do
  end if

  ! Can actually calculate the forces now
  force(:,:) = 0.0_dp

  ihydrogen = 0
  do iatom = 1, cell%natoms

    if (h_present .and. (cell%atom_species(iatom) .eq. h_num)) ihydrogen = ihydrogen + 1
    jhydrogen = 0

    do jatom = 1, cell%natoms
      if (iatom .ne. jatom) then

        contrib = d_f_rb(jatom)*d_rho_atomic(jatom,iatom) + d_f_rb(iatom)*d_rho_atomic(iatom,jatom) + d_phi(jatom,iatom)

        if (h_present) then
          if (cell%atom_species(iatom) .eq. h_num) then
            if(d_rho_atomic(iatom,jatom) > 0.0_dp) &
            contrib = contrib + d_f_rh(1,ihydrogen)*2.0_dp*d_rho_atomic(iatom,jatom) &
                            & + d_f_rh(2,ihydrogen)*0.5_dp/sqrt(d_rho_atomic(iatom,jatom))
          end if

          if (cell%atom_species(jatom) .eq. h_num) then
            jhydrogen = jhydrogen + 1
            if(d_rho_atomic(jatom,iatom) > 0.0_dp) &
            contrib = contrib + d_f_rh(1,jhydrogen)*2.0_dp*d_rho_atomic(jatom,iatom) &
                            & + d_f_rh(2,jhydrogen)*0.5_dp/sqrt(d_rho_atomic(jatom,iatom))
          end if
        end if

        force(:,iatom) = force(:,iatom) - contrib*rij_vec(:, iatom, jatom)
      end if
    end do
  end do

  ! rescale units:
  force(:,:) = units_natural_to_atomic(force=force(:,:))
end subroutine pot_eam_forces

function eam_phi(rij, atomi, atomj)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  integer,        intent(in)  ::  atomi, atomj
  real(kind=dp)               ::  eam_phi
  ! associated constants:
  real(kind=dp),  parameter   ::  c = 14.3888_dp  ! eV Angstrom

  eam_phi = c*eam_z(rij, atomi)*eam_z(rij, atomj)/rij
end function eam_phi

function eam_d_phi(rij, atomi, atomj)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  integer,        intent(in)  ::  atomi, atomj
  real(kind=dp)               ::  eam_d_phi
  ! associated constants:
  real(kind=dp),  parameter   ::  c = 14.3888_dp  ! eV Angstrom

  eam_d_phi = c*((eam_d_z(rij,atomi)*eam_z(rij,atomj) + eam_d_z(rij,atomj)*eam_z(rij,atomi))/rij &
  &           - (eam_z(rij,atomi)*eam_z(rij,atomj))/(rij**2))
end function eam_d_phi

function eam_z(rij, atomtype)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  integer,        intent(in)  ::  atomtype
  real(kind=dp)               ::  eam_z

  ! do no error checking for atom type - we do this on initialization
  select case (atomtype)
    case (h_num)
      ! hydrogen
      eam_z = eam_z_h(rij)*eam_s(rij)    ! multiply by smoothing function
    case (ni_num)
      ! nickel
      eam_z = eam_z_ni(rij)*eam_s(rij)   ! multiply by smoothing function
  end select
end function eam_z

function eam_s(rij)
  use constants,  only: pi, two_pi
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_s
  ! parameters:
  real(kind=dp),  parameter   ::  r_c = 5.0_dp    ! Angstrom
  real(kind=dp),  parameter   ::  delta = 5.0_dp  ! Angstrom
  real(kind=dp),  parameter   ::  r_c_p_delta = r_c+delta

  if (rij .le. r_c) then
    eam_s = 1.0_dp
  else if ((rij .gt. r_c) .and. (rij .lt. r_c_p_delta)) then
    eam_s = (rij-r_c_p_delta)/(-delta) - sin(pi*((2.0_dp*rij - 2.0_dp*r_c - delta)/delta))/two_pi
    if (eam_s .lt. 0.0_dp) eam_s = 0.0_dp   ! fix for when the above goes negative (specific values of 9.99...)
  else
    ! rij .ge. r_c+delta
    eam_s = 0.0_dp
  end if
end function eam_s

function eam_d_s(rij)
  use constants,  only: pi
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_d_s
  ! parameters:
  real(kind=dp),  parameter   ::  r_c = 5.0_dp    ! Angstrom
  real(kind=dp),  parameter   ::  delta = 5.0_dp  ! Angstrom
  real(kind=dp),  parameter   ::  r_c_p_delta = r_c+delta

  if (rij .le. r_c) then
    eam_d_s = 0.0_dp
  else if ((rij .gt. r_c) .and. (rij .lt. r_c_p_delta)) then
    eam_d_s = -(1.0_dp + cos(pi*(2.0_dp*rij - 2.0_dp*r_c - delta)/delta))/delta
  else
    ! rij .ge. r_c+delta
    eam_d_s = 0.0_dp
  end if
end function eam_d_s

function eam_d_z(rij, atomtype)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  integer,        intent(in)  ::  atomtype
  real(kind=dp)               ::  eam_d_z

  ! no error checking, assume correct atom types
  select case (atomtype)
    case (h_num)
      eam_d_z = eam_d_z_h(rij)*eam_s(rij) + eam_z_h(rij)*eam_d_s(rij)
    case (ni_num)
      eam_d_z = eam_d_z_ni(rij)*eam_s(rij) + eam_z_ni(rij)*eam_d_s(rij)
  end select
end function eam_d_z

function eam_z_ni(rij)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_z_ni
  ! associated constants:
  real(kind=dp),  parameter   ::  z0 = 10.0_dp  ! Dimensionless
  real(kind=dp),  parameter   ::  a = 0.537_dp  ! Angstrom
  real(kind=dp),  parameter   ::  b = 1.116_dp  ! Angstrom
! real(kind=dp),  parameter   ::  c = 1.0_dp    ! Dimensionless

  ! same as EAM5, but written as 1/a, 1/b (same constants)
! eam_z_ni = z0*(1.0_dp + (rij/b)**c)*algor_exp(-rij/a)
  eam_z_ni = z0*(1.0_dp + (rij/b))*algor_exp(-rij/a)
end function eam_z_ni

function eam_d_z_ni(rij)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_d_z_ni
  ! associated constants:
  real(kind=dp),  parameter   ::  z0 = 10.0_dp  ! Dimensionless
  real(kind=dp),  parameter   ::  a = 0.537_dp  ! Angstrom
  real(kind=dp),  parameter   ::  b = 1.116_dp  ! Angstrom
! real(kind=dp),  parameter   ::  c = 1.0_dp    ! Dimensionless

  ! same as EAM5, but written as 1/a, 1/b (same constants)

  !eam_d_z_ni = z0*algor_exp(-rij/a)*( ((rij/b)**c)*( (c/rij) - (1.0_dp/a) ) - (1.0_dp/a) )
  eam_d_z_ni = z0*algor_exp(-rij/a)*( (-1.0_dp/a) + (1.0_dp/b) - (rij/(a*b)) )
end function eam_d_z_ni

function eam_z_h(rij)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_z_h
  integer ::  j
  ! associated constants:
  real(kind=dp),                parameter   ::  r_sw    =  1.09_dp  ! Angstrom
  real(kind=dp),                parameter   ::  alpha_z = -5.957_dp ! Angstrom ^-5
  real(kind=dp),                parameter   ::  beta_z  =  9.804_dp ! Angstrom ^-4
  real(kind=dp),                parameter   ::  gamma_z = -3.715_dp ! Angstrom ^-3
  real(kind=dp),  dimension(3), parameter   ::  a = (/ 0.374_dp, 0.010_dp, 0.011_dp /)    ! Angstrom
  real(kind=dp),  dimension(3), parameter   ::  b = (/ 0.932_dp, 0.549_dp, 0.377_dp /)    ! Angstrom
  real(kind=dp),  dimension(3), parameter   ::  c = (/ 2.775_dp, 142.678_dp, 95.455_dp /) ! Dimensionless
  real(kind=dp),  dimension(3), parameter   ::  d = (/ 0.00752_dp, 8.0036_dp, 2.401_dp /) ! Dimensionless, Angstrom^-1, Angstrom

  if (rij .le. r_sw) then
    eam_z_h = 1.0_dp + alpha_z*rij**5 + beta_z*rij**4 + gamma_z*rij**3
  else ! rij > r_sw
    eam_z_h = 0.0_dp
    do j = 1, 3
      eam_z_h = eam_z_h + ((rij/b(j))**c(j))*algor_exp(-rij/a(j))
    end do
    eam_z_h = eam_z_h + d(1)*sech(d(2)*(rij-d(3)))
  end if
end function eam_z_h

function eam_d_z_h(rij)
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_d_z_h
  real(kind=dp)               ::  argument
  integer ::  j
  ! associated constants:
  real(kind=dp),                parameter   ::  r_sw    =  1.09_dp  ! Angstrom
  real(kind=dp),                parameter   ::  alpha_z = -5.957_dp ! Angstrom ^-5
  real(kind=dp),                parameter   ::  beta_z  =  9.804_dp ! Angstrom ^-4
  real(kind=dp),                parameter   ::  gamma_z = -3.715_dp ! Angstrom ^-3
  real(kind=dp),  dimension(3), parameter   ::  a = (/ 0.374_dp, 0.010_dp, 0.011_dp /)    ! Angstrom
  real(kind=dp),  dimension(3), parameter   ::  b = (/ 0.932_dp, 0.549_dp, 0.377_dp /)    ! Angstrom
  real(kind=dp),  dimension(3), parameter   ::  c = (/ 2.775_dp, 142.678_dp, 95.455_dp /) ! Dimensionless
  real(kind=dp),  dimension(3), parameter   ::  d = (/ 0.00752_dp, 8.0036_dp, 2.401_dp /) ! Dimensionless, Angstrom^-1, Angstrom

  if (rij .le. r_sw) then
    eam_d_z_h = 5.0_dp*alpha_z*rij**4 + 4.0_dp*beta_z*rij**3 + 3.0_dp*gamma_z*rij**2
  else ! rij > r_sw
    eam_d_z_h = 0.0_dp
    do j = 1, 3
      eam_d_z_h = eam_d_z_h + algor_exp(-rij/a(j))*((rij/b(j))**c(j))*(c(j)/rij - (1.0_dp/a(j)))
    end do
    argument = d(2)*(rij-d(3))
    eam_d_z_h = eam_d_z_h - d(1)*d(2)*tanh(argument)*sech(argument)
  end if
end function eam_d_z_h

function eam_f(rho_bar, rho_hat_1, rho_hat_2, atomtype)
  implicit none
  real(kind=dp) ::  eam_f
  real(kind=dp),  intent(in)  ::  rho_bar
  real(kind=dp),  intent(in)  ::  rho_hat_1
  real(kind=dp),  intent(in)  ::  rho_hat_2
  integer,        intent(in)  ::  atomtype

  select case (atomtype)
    case (h_num)
      ! hydrogen
      eam_f = eam_f_h(rho_bar, rho_hat_1, rho_hat_2)
    case (ni_num)
      ! nickel
      eam_f = eam_f_ni(rho_bar)
  end select
end function eam_f

function eam_d_f(rho_bar_or_hat, derivative, atomtype)
  implicit none
  real(kind=dp),    intent(in)  ::  rho_bar_or_hat
  integer,          intent(in)  ::  derivative
  integer,          intent(in)  ::  atomtype
  real(kind=dp)                 ::  eam_d_f

  select case (atomtype)
    case (h_num)
      eam_d_f = eam_d_f_h(rho_bar_or_hat, derivative)
    case (ni_num)
      eam_d_f = eam_d_f_ni(rho_bar_or_hat)
  end select
end function eam_d_f

function eam_f_ni(rho_bar)
  implicit none
  real(kind=dp),  intent(in)  ::  rho_bar
  real(kind=dp)               ::  eam_f_ni
  ! associated constants:
  real(kind=dp), parameter    ::  a_ni = -303.289_dp    ! eV Angstrom^3
  real(kind=dp), parameter    ::  alpha_ni = 7.647_dp   ! Angstrom^3
  real(kind=dp), parameter    ::  a_s = -2.015E10_dp    ! eV Angstrom^15
  real(kind=dp), parameter    ::  rho_c = 0.11_dp       ! Angstrom^-3
  real(kind=dp), parameter    ::  b_ni = 87.987_dp      ! eV^1/3 Angstrom^3
  real(kind=dp), parameter    ::  beta_ni = 75.075_dp   ! Angstrom^3
  real(kind=dp), parameter    ::  b_s = -5.384E8_dp     ! eV Angstrom^12
  real(kind=dp), parameter    ::  delta = 0.01_dp       ! Angstrom^-3
  real(kind=dp), parameter    ::  c_ni = -522.774_dp    ! eV Angstrom^3
  real(kind=dp), parameter    ::  gamma_ni = 373.379_dp ! Angstrom^3
  real(kind=dp), parameter    ::  c_s = -4.060E6_dp     ! eV Angstrom^9
  real(kind=dp), parameter    ::  d_ni = 39.421_dp      ! eV^1/5 Angstrom^3
  real(kind=dp), parameter    ::  delta_ni = 56.342_dp  ! Angstrom^3
  real(kind=dp), parameter    ::  d_s = -11.031_dp      ! eV
  real(kind=dp), parameter    ::  rho_c_m_delta = rho_c-delta
  real(kind=dp) ::  rho_bar_m_rho_c

  if ((rho_bar .ge. 0.0_dp) .and. (rho_bar .lt. rho_c_m_delta)) then
    ! EAM6 paper wrong here:
    eam_f_ni = a_ni*rho_bar*algor_exp(-alpha_ni*rho_bar) + ((b_ni*rho_bar)**3)*algor_exp(-beta_ni*rho_bar) &
    &         + c_ni*rho_bar*algor_exp(-gamma_ni*rho_bar) + ((d_ni*rho_bar)**5)*algor_exp(-delta_ni*rho_bar)
  else if ((rho_bar .ge. rho_c_m_delta) .and. (rho_bar .le. rho_c)) then
    rho_bar_m_rho_c = rho_bar - rho_c
    eam_f_ni = a_s*(rho_bar_m_rho_c)**5 + b_s*(rho_bar_m_rho_c)**4 + c_s*(rho_bar_m_rho_c)**3 + d_s
  else
    ! ie: rho_bar > rho_c
    eam_f_ni = d_s
  end if
end function eam_f_ni

function eam_d_f_ni(rho_bar)
  implicit none
  real(kind=dp),  intent(in)  ::  rho_bar
  real(kind=dp)               ::  eam_d_f_ni
  ! associated constants:
  real(kind=dp), parameter    ::  a_ni = -303.289_dp    ! eV Angstrom^3
  real(kind=dp), parameter    ::  alpha_ni = 7.647_dp   ! Angstrom^3
  real(kind=dp), parameter    ::  a_s = -2.015E10_dp    ! eV Angstrom^15
  real(kind=dp), parameter    ::  rho_c = 0.11_dp       ! Angstrom^-3
  real(kind=dp), parameter    ::  b_ni = 87.987_dp      ! eV^1/3 Angstrom^3
  real(kind=dp), parameter    ::  beta_ni = 75.075_dp   ! Angstrom^3
  real(kind=dp), parameter    ::  b_s = -5.384E8_dp     ! eV Angstrom^12
  real(kind=dp), parameter    ::  delta = 0.01_dp       ! Angstrom^-3
  real(kind=dp), parameter    ::  c_ni = -522.774_dp    ! eV Angstrom^3
  real(kind=dp), parameter    ::  gamma_ni = 373.379_dp ! Angstrom^3
  real(kind=dp), parameter    ::  c_s = -4.060E6_dp     ! eV Angstrom^9
  real(kind=dp), parameter    ::  d_ni = 39.421_dp      ! eV^1/5 Angstrom^3
  real(kind=dp), parameter    ::  delta_ni = 56.342_dp  ! Angstrom^3
  real(kind=dp), parameter    ::  rho_c_m_delta = rho_c-delta
  real(kind=dp) ::  rho_bar_m_rho_c

  if ((rho_bar .ge. 0.0_dp) .and. (rho_bar .lt. rho_c_m_delta)) then
    ! EAM6 paper wrong here:
    eam_d_f_ni = a_ni*algor_exp(-alpha_ni*rho_bar)*(1.0_dp - alpha_ni*rho_bar) &
    & + (b_ni**3)*algor_exp(-beta_ni*rho_bar)*(3.0_dp*rho_bar**2 - beta_ni*rho_bar**3) &
    & + c_ni*algor_exp(-gamma_ni*rho_bar)*(1.0_dp - gamma_ni*rho_bar)  &
    & + (d_ni**5)*algor_exp(-delta_ni*rho_bar)*(5.0_dp*rho_bar**4 - delta_ni*rho_bar**5)
  else if ((rho_bar .ge. rho_c_m_delta) .and. (rho_bar .le. rho_c)) then
    rho_bar_m_rho_c = rho_bar - rho_c
    eam_d_f_ni = 5.0_dp*a_s*(rho_bar_m_rho_c)**4 + 4.0_dp*b_s*(rho_bar_m_rho_c)**3 + 3.0_dp*c_s*(rho_bar_m_rho_c)**2
  else
    ! ie: rho_bar > rho_c
    eam_d_f_ni = 0.0_dp
  end if
end function eam_d_f_ni

function eam_f_h(rho_bar, rho_hat_1, rho_hat_2)
  implicit none
  real(kind=dp),  intent(in)  ::  rho_bar
  real(kind=dp),  intent(in)  ::  rho_hat_1
  real(kind=dp),  intent(in)  ::  rho_hat_2
  real(kind=dp)               ::  eam_f_h
  ! associated constants:
  real(kind=dp),  parameter   ::  alpha_h_1 = -4.803_dp     ! eV Ang^3*gamma_h_1
  real(kind=dp),  parameter   ::  alpha_h_2 = -74.556_dp    ! eV Ang^3*gamm_h_2
  real(kind=dp),  parameter   ::  a_h_1     =  1.841E11_dp  ! eV Ang^15
  real(kind=dp),  parameter   ::  d_h_1     = -1.304E11_dp  ! eV Ang^15
  real(kind=dp),  parameter   ::  rho_cut_3 =  0.009_dp     ! Ang^-3
  real(kind=dp),  parameter   ::  beta_h_1  =  34.46_dp     ! Ang^-3
  real(kind=dp),  parameter   ::  beta_h_2  =  7.2784_dp    ! Ang^-3
  real(kind=dp),  parameter   ::  b_h_1     =  3.276E9_dp   ! eV Ang^12
  real(kind=dp),  parameter   ::  e_h_1     =  2.376E9_dp   ! eV Ang^12
  real(kind=dp),  parameter   ::  delta_3   =  0.007_dp     ! Ang^-3
  real(kind=dp),  parameter   ::  gamma_h_1 =  0.3855_dp    ! Dimensionless
  real(kind=dp),  parameter   ::  gamma_h_2 =  1.0154_dp    ! Dimensionless
  real(kind=dp),  parameter   ::  c_h_1     =  1.510E7_dp   ! eV Ang^9
  real(kind=dp),  parameter   ::  f_h_1     = -1.195E7_dp   ! eV Ang^9
  real(kind=dp),  parameter   ::  rho_cut_3_m_delta_3 = rho_cut_3 - delta_3 ! Ang^-3
  ! Fortran is column major (row, column)
  real(kind=dp),  dimension(3,2), parameter ::  &
  ! xi:       11 (eV),   12 (Ang^6), 13 (Ang-6),   21 (eV),   22 (Ang^3/2),  23 (Ang^-3/2)
  & xi = reshape((/ -0.04084_dp, 3.686E8_dp, 13563.0_dp, -0.00325_dp,   54.772_dp,     0.812_dp /), shape(xi))
  real(kind=dp),  dimension(2) :: rho_hat
  real(kind=dp) ::  delta_rho_5, delta_rho_6
  integer :: k

  ! primary term:
  if ((rho_bar .ge. 0.0_dp) .and. (rho_bar .le. rho_cut_3_m_delta_3)) then
    eam_f_h = (alpha_h_1*rho_bar**gamma_h_1)*algor_exp(-beta_h_1*rho_bar)
  else if ((rho_bar .gt. rho_cut_3_m_delta_3) .and. (rho_bar .le. rho_cut_3)) then
    delta_rho_5 = rho_bar - rho_cut_3
    delta_rho_6 = rho_bar - rho_cut_3_m_delta_3

    eam_f_h = a_h_1*delta_rho_5**5 + b_h_1*delta_rho_5**4 + c_h_1*delta_rho_5**3 &
      &     + d_h_1*delta_rho_6**5 + e_h_1*delta_rho_6**4 + f_h_1*delta_rho_6**3
  else ! rho_bar > rho_cut_3
    eam_f_h = (alpha_h_2*rho_bar**gamma_h_2)*algor_exp(-beta_h_2*rho_bar)
  end if

  if (non_local_correction) then
    ! non-local term:
    rho_hat(:) = (/ rho_hat_1, rho_hat_2 /)
    do k = 1,2
      eam_f_h = eam_f_h + xi(1,k)*sech(xi(2,k)*(rho_hat(k)-xi(3,k)))
    end do
  end if
end function eam_f_h

function eam_d_f_h(rho_bar_or_hat, derivative)
  use io, only: io_err
  implicit none
  real(kind=dp),    intent(in)  ::  rho_bar_or_hat
  integer,          intent(in)  ::  derivative
  real(kind=dp)                 ::  eam_d_f_h
  ! associated constants:
  real(kind=dp),  parameter   ::  alpha_h_1 = -4.803_dp     ! eV Ang^3*gamma_h_1
  real(kind=dp),  parameter   ::  alpha_h_2 = -74.556_dp    ! eV Ang^3*gamm_h_2
  real(kind=dp),  parameter   ::  a_h_1     =  1.841E11_dp  ! eV Ang^15
  real(kind=dp),  parameter   ::  d_h_1     = -1.304E11_dp  ! eV Ang^15
  real(kind=dp),  parameter   ::  rho_cut_3 =  0.009_dp     ! Ang^-3
  real(kind=dp),  parameter   ::  beta_h_1  =  34.46_dp     ! Ang^-3
  real(kind=dp),  parameter   ::  beta_h_2  =  7.2784_dp    ! Ang^-3
  real(kind=dp),  parameter   ::  b_h_1     =  3.276E9_dp   ! eV Ang^12
  real(kind=dp),  parameter   ::  e_h_1     =  2.376E9_dp   ! eV Ang^12
  real(kind=dp),  parameter   ::  delta_3   =  0.007_dp     ! Ang^-3
  real(kind=dp),  parameter   ::  gamma_h_1 =  0.3855_dp    ! Dimensionless
  real(kind=dp),  parameter   ::  gamma_h_2 =  1.0154_dp    ! Dimensionless
  real(kind=dp),  parameter   ::  c_h_1     =  1.510E7_dp   ! eV Ang^9
  real(kind=dp),  parameter   ::  f_h_1     = -1.195E7_dp   ! eV Ang^9
  real(kind=dp),  parameter   ::  rho_cut_3_m_delta_3 = rho_cut_3 - delta_3 ! Ang^-3
  ! Fortran is column major (row, column)
  real(kind=dp),  dimension(3,2), parameter ::  &
  ! xi:       11 (eV),   12 (Ang^6), 13 (Ang-6),   21 (eV),   22 (Ang^3/2),  23 (Ang^-3/2)
  & xi = reshape((/ -0.04084_dp, 3.686E8_dp, 13563.0_dp, -0.00325_dp,   54.772_dp,     0.812_dp /), shape(xi))
  integer :: k
  real(kind=dp) :: rho_bar
  real(kind=dp) :: delta_rho_5, delta_rho_6
  real(kind=dp) :: argument

  k = 0
  select case (derivative)
    case (0)    ! rho_bar
      ! primary term:
      rho_bar = rho_bar_or_hat
      if ((rho_bar .ge. 0.0_dp) .and. (rho_bar .le. rho_cut_3_m_delta_3)) then
        eam_d_f_h = alpha_h_1*algor_exp(-beta_h_1*rho_bar)*(gamma_h_1*rho_bar**(gamma_h_1 - 1.0_dp) - beta_h_1*rho_bar**gamma_h_1)
      else if ((rho_bar .gt. rho_cut_3_m_delta_3) .and. (rho_bar .le. rho_cut_3)) then
        delta_rho_5 = rho_bar - rho_cut_3
        delta_rho_6 = rho_bar - rho_cut_3_m_delta_3

        eam_d_f_h = (5.0_dp*a_h_1*delta_rho_5**4 + 4.0_dp*b_h_1*delta_rho_5**3 + 3.0_dp*c_h_1*delta_rho_5**2)*(-rho_cut_3) &
         &   + (5.0_dp*d_h_1*delta_rho_6**4 + 4.0_dp*e_h_1*delta_rho_6**3 + 3.0_dp*f_h_1*delta_rho_6**2)*(-rho_cut_3_m_delta_3)

      else ! rho_bar > rho_cut_3
        eam_d_f_h = alpha_h_2*algor_exp(-beta_h_2*rho_bar)*(gamma_h_2*rho_bar**(gamma_h_2 - 1.0_dp) - beta_h_2*rho_bar**gamma_h_2)
      end if
    case (1,2)  ! rho_hat_1/rho_hat_2
      k = derivative
    case default
      call io_err("Error in potential - eam_d_f_h: Unexpected derivative")
  end select

  ! non-local term:
  if (k .ne. 0) then
    if (non_local_correction) then
      ! k is either 1, 2 or 0
      argument = xi(2,k)*(rho_bar_or_hat-xi(3,k))
      eam_d_f_h = -xi(1,k)*xi(2,k)*tanh(argument)*sech(argument)
    else
      eam_d_f_h = 0.0_dp
    end if
  end if
end function eam_d_f_h

function eam_rho(rij, atomtype)
  implicit none
  real(kind=dp)               ::  eam_rho
  real(kind=dp),  intent(in)  ::  rij
  integer,        intent(in)  ::  atomtype

  ! no error checking for atom type - we do this on initialization
  select case (atomtype)
    case (h_num)
      ! Hydrogen
      eam_rho = eam_rho_h(rij)*eam_s(rij)
    case (ni_num)
      ! Nickel
      eam_rho = eam_rho_ni(rij)*eam_s(rij)
  end select
end function eam_rho

function eam_d_rho(rij, atomtype)
  implicit none
  real(kind=dp)               ::  eam_d_rho
  real(kind=dp),  intent(in)  ::  rij
  integer,        intent(in)  ::  atomtype

  ! no error checking for atom type - we do this on initialization
  select case (atomtype)
    case (h_num)
      ! Hydrogen
      eam_d_rho = eam_d_rho_h(rij)*eam_s(rij) + eam_rho_h(rij)*eam_d_s(rij)
    case (ni_num)
      ! Nickel
      eam_d_rho = eam_d_rho_ni(rij)*eam_s(rij) + eam_rho_ni(rij)*eam_d_s(rij)
  end select
end function eam_d_rho

function eam_rho_h(rij)
  use constants,  only: pi, units_atomic_to_natural
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_rho_h
  real(kind=dp)               ::  bohr_radius

  bohr_radius = units_atomic_to_natural(length=1.0_dp)
  eam_rho_h = algor_exp(-2.0_dp*rij/bohr_radius)/(pi*bohr_radius**3)
end function eam_rho_h

function eam_d_rho_h(rij)
  use constants,  only: pi, units_atomic_to_natural
  implicit none
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp)               ::  eam_d_rho_h
  real(kind=dp)               ::  bohr_radius

  bohr_radius = units_atomic_to_natural(length=1.0_dp)
  eam_d_rho_h = -2.0_dp*algor_exp(-2.0_dp*rij/bohr_radius)/(pi*bohr_radius**4)
end function eam_d_rho_h

function eam_rho_ni(rij)
  implicit none
  real(kind=dp)               ::  eam_rho_ni
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp),  parameter   ::  n_s = 2.0_dp
  real(kind=dp),  parameter   ::  n = 10.0_dp
  real(kind=dp),  parameter   ::  n_m_ns = n-n_s

  eam_rho_ni = n_s*eam_rho_shell_ni(rij,'s') + n_m_ns*eam_rho_shell_ni(rij,'d')
end function eam_rho_ni

function eam_d_rho_ni(rij)
  implicit none
  real(kind=dp)               ::  eam_d_rho_ni
  real(kind=dp),  intent(in)  ::  rij
  real(kind=dp),  parameter   ::  n_s = 2.0_dp
  real(kind=dp),  parameter   ::  n = 10.0_dp
  real(kind=dp),  parameter   ::  n_m_ns = n-n_s

  eam_d_rho_ni = n_s*eam_d_rho_shell_ni(rij,'s') + n_m_ns*eam_d_rho_shell_ni(rij,'d')
end function eam_d_rho_ni

function eam_rho_shell_ni(rij, shelltype)
  use constants,  only: pi
  use io,         only: io_err
  implicit none
  real(kind=dp)                 ::  eam_rho_shell_ni
  real(kind=dp),    intent(in)  ::  rij
  character(len=1), intent(in)  ::  shelltype
  integer ::  i
  ! parameters:
  real(kind=dp),  parameter   ::  over_4pi = 1.0_dp/(4.0_dp*pi)
  ! s shell
  real(kind=dp),  dimension(8), parameter ::  c_s = (/ -0.00389_dp, -0.02991_dp, -0.03189_dp, 0.15289_dp, -0.20048_dp, &
  & -0.05423_dp, 0.49292_dp, 0.61875_dp /)
  ! d shell
  real(kind=dp),  dimension(2), parameter ::  c_d = (/ 0.42120_dp, 0.70658_dp /)

  eam_rho_shell_ni = 0.0_dp

  select case (shelltype)
    case ('s')
      do i = 1,8
        eam_rho_shell_ni = eam_rho_shell_ni + c_s(i)*eam_r_shell_ni(rij,i,'s')
      end do
    case ('d')
      do i = 1,2
        eam_rho_shell_ni = eam_rho_shell_ni + c_d(i)*eam_r_shell_ni(rij,i,'d')
      end do
    case default
      call io_err("Error in potential: eam_rho_shell_ni - unexpected shell type")
  end select

  eam_rho_shell_ni = over_4pi*eam_rho_shell_ni**2
end function eam_rho_shell_ni

function eam_d_rho_shell_ni(rij, shelltype)
  use constants,  only: two_pi
  use io,         only: io_err
  implicit none
  real(kind=dp)                 ::  eam_d_rho_shell_ni
  real(kind=dp),    intent(in)  ::  rij
  character(len=1), intent(in)  ::  shelltype
  integer ::  i
  real(kind=dp) :: ciri
  ! parameters:
  real(kind=dp),  parameter   ::  over_2pi = 1.0_dp/two_pi
  ! s shell
  real(kind=dp),  dimension(8), parameter ::  c_s = (/ -0.00389_dp, -0.02991_dp, -0.03189_dp, 0.15289_dp, -0.20048_dp, &
  & -0.05423_dp, 0.49292_dp, 0.61875_dp /)
  ! d shell
  real(kind=dp),  dimension(2), parameter ::  c_d = (/ 0.42120_dp, 0.70658_dp /)

  eam_d_rho_shell_ni = 0.0_dp
  ciri = 0.0_dp

  select case (shelltype)
    case ('s')
      do i = 1,8
        ciri = ciri + c_s(i)*eam_r_shell_ni(rij,i,'s')
        eam_d_rho_shell_ni = eam_d_rho_shell_ni + c_s(i)*eam_d_r_shell_ni(rij,i,'s')
      end do
    case ('d')
      do i = 1,2
        ciri = ciri + c_d(i)*eam_r_shell_ni(rij,i,'d')
        eam_d_rho_shell_ni = eam_d_rho_shell_ni + c_d(i)*eam_d_r_shell_ni(rij,i,'d')
      end do
    case default
      call io_err("Error in potential: eam_d_rho_shell_ni - unexpected shell type")
  end select

  ! might need abs here
  eam_d_rho_shell_ni = over_2pi*eam_d_rho_shell_ni*ciri
end function eam_d_rho_shell_ni

function eam_r_shell_ni(rij, i, shelltype)
  use io,     only: io_err
  implicit none
  real(kind=dp)                 ::  eam_r_shell_ni
  real(kind=dp),    intent(in)  ::  rij
  integer,          intent(in)  ::  i
  character(len=1), intent(in)  ::  shelltype
  integer ::  n_i
  ! parameters:                               zeta is in Ang^-1
  ! s shell:
  real(kind=dp),  dimension(8), parameter ::  zeta_s = (/ 54.88885_dp, 38.48431_dp, 27.42703_dp, 20.88204_dp, 10.95707_dp, &
  & 7.31958_dp, 3.92650_dp, 2.15289_dp /)
  ! d shell:
  real(kind=dp),  dimension(2), parameter ::  zeta_d = (/ 12.67582_dp, 5.43253_dp /)
  real(kind=dp),  dimension(4), parameter ::  sqrt_factorial_2ni = (/ sqrt(2.0_dp),   sqrt(24.0_dp), &
                                                                  &   sqrt(720.0_dp), sqrt(40320.0_dp) /)

  select case (shelltype)
    case ('s')
      ! assume i is between 1 and 8 (inclusive) - don't do error checking here
      ! make use of integer division
      n_i = (i+1)/2

      eam_r_shell_ni = ((2.0_dp*zeta_s(i))**(real(n_i,dp)+0.5_dp))*(rij**(n_i-1))*(algor_exp(-zeta_s(i)*rij)) &
      & / sqrt_factorial_2ni(n_i)
    case ('d')
      ! again, don't check i - assume everything is okay
      !n_i = 3 ! always
      ! same functional form as above, except we know some values
      eam_r_shell_ni = ((2.0_dp*zeta_d(i))**3.5_dp)*(rij**2)*(algor_exp(-zeta_d(i)*rij))/sqrt_factorial_2ni(3)
    case default
      call io_err("Error in potential: eam_r_shell_ni - unexpected shell type")
  end select
end function eam_r_shell_ni

function eam_d_r_shell_ni(rij, i, shelltype)
  use algor,  only: algor_factorial
  use io,     only: io_err
  implicit none
  real(kind=dp)                 ::  eam_d_r_shell_ni
  real(kind=dp),    intent(in)  ::  rij
  integer,          intent(in)  ::  i
  character(len=1), intent(in)  ::  shelltype
  integer ::  n_i
  ! parameters:                               zeta is in Ang^-1
  ! s shell:
  real(kind=dp),  dimension(8), parameter ::  zeta_s = (/ 54.88885_dp, 38.48431_dp, 27.42703_dp, 20.88204_dp, 10.95707_dp, &
  & 7.31958_dp, 3.92650_dp, 2.15289_dp /)
  ! d shell:
  real(kind=dp),  dimension(2), parameter ::  zeta_d = (/ 12.67582_dp, 5.43253_dp /)
  real(kind=dp),  dimension(4), parameter ::  sqrt_factorial_2ni = (/ sqrt(2.0_dp),   sqrt(24.0_dp), &
                                                                  &   sqrt(720.0_dp), sqrt(40320.0_dp) /)
  select case (shelltype)
    case ('s')
      ! assume i is between 1 and 8 (inclusive) - don't do error checking here
      ! make use of integer division
      n_i = (i+1)/2

      eam_d_r_shell_ni = ((2.0_dp*zeta_s(i))**(real(n_i,dp)+0.5_dp))* &
      & (algor_exp(-zeta_s(i)*rij)*(real((n_i-1),dp)*(rij**(n_i-2)) - zeta_s(i)*(rij**(n_i-1)))) &
      & / sqrt_factorial_2ni(n_i)
    case ('d')
      ! again, don't check i - assume everything is okay
      !n_i = 3 ! always
      ! same functional form as above, except we know some values
      eam_d_r_shell_ni = ((2.0_dp*zeta_d(i))**3.5_dp)* &
      & (algor_exp(-zeta_d(i)*rij)*(2.0_dp*rij - zeta_d(i)*(rij**2)))/sqrt_factorial_2ni(3)
    case default
      call io_err("Error in potential: eam_d_r_shell_ni - unexpected shell type")
  end select
end function eam_d_r_shell_ni

function sech(x)
  implicit none
  real(kind=dp),  intent(in)  ::  x
  real(kind=dp)               ::  sech
  if(x > -500) then
    sech = 1.0_dp/cosh(x)
  else
    sech = 0
  end if
end function sech
end module potential


    Subroutine open_ascii(IBI,CSTAT,Cfile)
    implicit none
!   Open ascii-file (text file)
    INTEGER IBI
    CHARACTER*(*) Cfile,CSTAT
!
!   Open requested file
!
     OPEN (UNIT=IBI,FILE=Cfile,STATUS=CSTAT,ERR=1900)
!
     RETURN
!*********************************************************
 1900 WRITE (*,*) 'FILE..:',trim(cfile)
      WRITE (*,*) 'Error opening requested file.',' Program stopped.'
    STOP
    END


text:
!
! Copyright (c) Stanford University, The Regents of the University of
!               California, and others.
!
! All Rights Reserved.
!
! See Copyright-SimVascular.txt for additional details.
!
! Permission is hereby granted, free of charge, to any person obtaining
! a copy of this software and associated documentation files (the
! "Software"), to deal in the Software without restriction, including
! without limitation the rights to use, copy, modify, merge, publish,
! distribute, sublicense, and/or sell copies of the Software, and to
! permit persons to whom the Software is furnished to do so, subject
! to the following conditions:
!
! The above copyright notice and this permission notice shall be included
! in all copies or substantial portions of the Software.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
! IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
! TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
! PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
! OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
! EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
! PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
! PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
! LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
! NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
!
!--------------------------------------------------------------------
!
!     Parameters used for svFSI.
!
!--------------------------------------------------------------------

!     Maximum possible nsd
      INTEGER(KIND=IKIND), PARAMETER :: maxnsd = 3
!--------------------------------------------------------------------
!     Maximum number of properties
      INTEGER(KIND=IKIND), PARAMETER :: maxNProp = 20
!--------------------------------------------------------------------
!     Gauss point coordinates, upto 5 points
      REAL(KIND=RKIND), PARAMETER :: gW(5,5) =
     2   RESHAPE( (/2._RKIND, 0._RKIND, 0._RKIND, 0._RKIND, 0._RKIND,
     3   1._RKIND, 1._RKIND, 0._RKIND, 0._RKIND, 0._RKIND,
     4   0.5555555555555556_RKIND, 0.8888888888888889_RKIND,
     4   0.5555555555555556_RKIND, 0._RKIND, 0._RKIND,
     5   0.3478548451374538_RKIND, 0.6521451548625462_RKIND,
     5   0.6521451548625462_RKIND, 0.3478548451374538_RKIND, 0._RKIND,
     6   0.2369268850561890_RKIND, 0.4786286704993665_RKIND,
     6   0.5688888888888889_RKIND, 0.4786286704993665_RKIND,
     6   0.2369268850561890_RKIND/) , (/5, 5/) )
!--------------------------------------------------------------------
!     Gauss point weights, upto 5 points
      REAL(KIND=RKIND), PARAMETER :: gXi(5,5) =
     2   RESHAPE( (/0._RKIND, 0._RKIND, 0._RKIND, 0._RKIND, 0._RKIND,
     3   -0.57735026918962584_RKIND, 0.57735026918962584_RKIND,
     3   0._RKIND, 0._RKIND, 0._RKIND,
     4   -0.7745966692414834_RKIND, 0._RKIND, 0.7745966692414834_RKIND,
     4   0._RKIND, 0._RKIND,
     5   -0.86113631159405257_RKIND, -0.33998104358485631_RKIND,
     5   0.33998104358485631_RKIND, 0.86113631159405257_RKIND, 0._RKIND,
     6   -0.90617984593866396_RKIND, -0.53846931010568311_RKIND,
     6   0._RKIND, 0.53846931010568311_RKIND,
     6   0.90617984593866396_RKIND/) , (/5, 5/) )
!--------------------------------------------------------------------
      CHARACTER, PARAMETER :: delimiter = "/"
!--------------------------------------------------------------------
!     Possible physical properties. Current maxNPror is 20.
!     When adding more properties, remember to increase maxNProp
!     Density of fluid, density of solid, viscosity of solid, elasticity
!     modulus, Poisson's ratio, conductivity, internal force (X, Y, Z),
!     stabilization coeff. for backflow divergence, external source,
!     damping, shell thickness, stabilization coeffs. for USTRUCT
!     (momentum, continuity)
      INTEGER(KIND=IKIND), PARAMETER :: prop_NA = 0, fluid_density = 1,
     2   solid_density = 2, solid_viscosity = 3, elasticity_modulus = 4,
     3   poisson_ratio = 5, conductivity = 6, f_x = 7, f_y = 8, f_z = 9,
     4   backflow_stab = 10, source_term = 11, damping = 12,
     5   shell_thickness = 13, ctau_M = 14, ctau_C = 15
!--------------------------------------------------------------------
!     Types of accepted elements
!     Point, Line (linear), Line (quadratic), Triangle (linear),
!     Triangle (quadratic), Quads (bilinear), Quads (serendipity),
!     Quads (biquadratic), Tetrahedron (linear), Tets (quadratic),
!     Hexgonal bricks (trilinear), Hex (quadratic/serendipity),
!     Hex (triquadratic), Wedge, NURBS
      INTEGER(KIND=IKIND), PARAMETER :: eType_NA = 100, eType_PNT = 101,
     2   eType_LIN1 = 102, eType_LIN2 = 103, eType_TRI3 = 104,
     3   eType_TRI6 = 105, eType_QUD4 = 106, eType_QUD8 = 107,
     4   eType_QUD9 = 108, eType_TET4 = 109, eType_TET10 = 110,
     5   eType_HEX8 = 111, eType_HEX20 = 112, eType_HEX27 = 113,
     6   eType_WDG = 114, eType_NRB = 115
!--------------------------------------------------------------------
!     Types of equations that are included in this solver
!     Fluid equation (Navier-Stokes), nonlinear structure (pure d), heat
!     equation, linear elasticity, heat in fluid (advection-diffusion),
!     fluid-structure-interaction, mesh motion, Shell mechanics,
!     Coupled-Momentum-Method, Cardiac Electro-Physiology,
!     Nonlinear structure (v-p), Stokes equations
      INTEGER(KIND=IKIND), PARAMETER :: phys_NA = 200, phys_fluid = 201,
     2   phys_struct = 202, phys_heatS = 203, phys_lElas = 204,
     3   phys_heatF = 205, phys_FSI = 206, phys_mesh = 207,
     4   phys_shell = 208, phys_CMM = 209, phys_CEP = 210,
     5   phys_ustruct = 211, phys_stokes = 212
!--------------------------------------------------------------------
!     Differenty type of coupling for cplBC
!     Not-available, implicit, semi-implicit, and explicit
      INTEGER(KIND=IKIND), PARAMETER :: cplBC_NA = 400, cplBC_I = 401,
     2   cplBC_SI = 402, cplBC_E = 403
!--------------------------------------------------------------------
!     cplBC type of coupling to between 3D and OD-LPN models:
!     Dirichlet type coupling, Neumann type coupling
      INTEGER(KIND=IKIND), PARAMETER :: cplBC_Dir = 66112,
     2   cplBC_Neu = 66113
!--------------------------------------------------------------------
!     boundary conditions types. Items of this list can be combined
!     BCs from imposing perspective can be Neu/Dir/per
!     BCs time dependence can be std/ustd/cpl/gen/res
!     BCs spatial distribution can be flat/para/ud
!     Beside these nodes at the boundary perimeter can be set to
!     zero and flux through surface can be assigned instead of nodal
!     values.
!     Dirichlet, Neumann, Traction, CMM, Robin, steady, unsteady,
!     coupled, general (combination of ud/ustd), resistance, imposed
!     flux, zero out perimeter, impose BC on the integral of state
!     variable or D (instead of Y), flat profile, parabolic profile,
!     user defined profile, backflow stabilization, BCs for shells
!     (fixed, hinged, free, symmetric), undeforming Neu, RCR-Neu
      INTEGER(KIND=IKIND), PARAMETER :: bType_Dir = 0, bType_Neu = 1,
     2   bType_trac = 2, bType_CMM = 3, bType_Robin = 4, bType_std = 5,
     3   bType_ustd = 6, bType_cpl = 7, bType_gen = 8, bType_res = 9,
     4   bType_flx = 10, bType_zp = 11, bType_impD = 12, bType_flat =13,
     5   bType_para = 14, bType_ud = 15, bType_bfs = 16, bType_fix = 17,
     6   bType_hing = 18, bType_free = 19, bType_symm = 20,
     7   bType_undefNeu = 21, bType_RCR = 22
!--------------------------------------------------------------------
!     Body force types: volumetric (default), traction, Neumann
!     (pressure based), time dependence (steady, unsteady, spatially
!     varying, general)
      INTEGER(KIND=IKIND), PARAMETER :: bfType_vol = 0, bfType_trac = 1,
     2   bfType_Neu = 2, bfType_std = 3, bfType_ustd = 4,
     3   bfType_spl = 5, bfType_gen = 6
!--------------------------------------------------------------------
!     Possible senarios for the output, followed by the possible outputs
!     Undefined output, extract from A, extract from Y, extract from D,
!     simple integral, WSS, traction, vorticity, vortex (lambda_ci),
!     strain invariants (fluid), energy flux, heat flux, absolute
!     velocity (for FSI), fiber directions, fiber alignment, 2nd Piola-
!     Kirchhoff stress, Cauchy stress, von Mises stress, Jacobian,
!     Def. grad. tensor, Green strain, divergence of velocity, viscosity
      INTEGER(KIND=IKIND), PARAMETER :: outGrp_NA = 500, outGrp_A = 501,
     2   outGrp_Y = 502, outGrp_D = 503, outGrp_I = 504, outGrp_WSS =
     3   505, outGrp_trac = 506, outGrp_vort = 507, outGrp_vortex = 508,
     4   outGrp_stInv = 509, outGrp_eFlx = 510, outGrp_hFlx = 511,
     5   outGrp_absV = 512, outGrp_fN = 513, outGrp_fA = 514,
     6   outGrp_stress = 515, outGrp_cauchy = 516, outGrp_mises = 517,
     7   outGrp_J = 518, outGrp_F = 519, outGrp_strain = 520,
     8   outGrp_divV = 521, outGrp_Visc = 522
!--------------------------------------------------------------------
      INTEGER(KIND=IKIND), PARAMETER :: out_velocity = 599,
     2   out_pressure = 598, out_temperature = 597, out_voltage = 596,
     3   out_acceleration = 595, out_displacement = 594, out_integ =593,
     4   out_WSS = 592, out_traction = 591, out_vorticity = 590,
     5   out_vortex = 589, out_strainInv = 588, out_energyFlux = 587,
     6   out_heatFlux = 586, out_absVelocity = 585, out_fibDir = 584,
     7   out_fibAlign = 583, out_stress = 582, out_cauchy = 581,
     8   out_mises = 580, out_jacobian = 579, out_defGrad = 578,
     9   out_strain = 577, out_divergence = 576, out_viscosity = 575
!--------------------------------------------------------------------
!     Mesher choice for remeshing for moving wall problems
      INTEGER(KIND=IKIND), PARAMETER :: RMSH_TETGEN = 1,
     2   RMSH_MESHSIM = 2
!--------------------------------------------------------------------
!     Type of constitutive model (isochoric) for structure equation:
!     St.Venant-Kirchhoff, modified St.Venant-Kirchhoff, NeoHookean,
!     Mooney-Rivlin, modified Holzapfel-Gasser-Ogden with dispersion,
!     Linear model (S = mu*I), Guccione (1995), Holzapfel & Ogden model
!     for myocardium (2009)
      INTEGER(KIND=IKIND), PARAMETER :: stIso_NA = 600,
     2   stIso_StVK = 601, stIso_mStVK = 602, stIso_nHook = 603,
     3   stIso_MR = 604, stIso_HGO = 605, stIso_lin = 606,
     4   stIso_Gucci = 607, stIso_HO = 608
!--------------------------------------------------------------------
!     Type of constitutive model (volumetric) for structure eqn:
!     Quadratic, Simo-Taylor91, Miehe94
      INTEGER(KIND=IKIND), PARAMETER :: stVol_NA = 650,
     2   stVol_Quad = 651, stVol_ST91 = 652, stVol_M94 = 653
!--------------------------------------------------------------------
!     Type of fluid viscosity: constant, Carreau-Yasuda shear-thinning
!     model, Cassons non-Newtonian model
      INTEGER(KIND=IKIND), PARAMETER :: viscType_NA = 699,
     2   viscType_Const = 698, viscType_CY = 697, viscType_Cass = 696
!--------------------------------------------------------------------
!     Preconditioner definitions
      INTEGER(KIND=IKIND), PARAMETER :: PREC_NONE = 700,
     2   PREC_FSILS = 701, PREC_TRILINOS_DIAGONAL = 702,
     3   PREC_TRILINOS_BLOCK_JACOBI = 703, PREC_TRILINOS_ILU = 704,
     4   PREC_TRILINOS_ILUT = 705, PREC_TRILINOS_IC = 706,
     5   PREC_TRILINOS_ICT = 707, PREC_TRILINOS_ML = 708,
     6   PREC_RCS = 709
!--------------------------------------------------------------------
!     Solver definitions
      INTEGER(KIND=IKIND), PARAMETER :: lSolver_NA = 799,
     2   lSolver_CG = 798, lSolver_GMRES=797, lSolver_NS=796,
     3   lSolver_BICGS = 795
!--------------------------------------------------------------------
!     Contact model
      INTEGER(KIND=IKIND), PARAMETER :: cntctM_NA = 800,
     2   cntctM_penalty = 801
!--------------------------------------------------------------------
!     IB method: traditional immersed finite element (IFEM)
      INTEGER(KIND=IKIND), PARAMETER :: ibMthd_NA = 850,
     2  ibMthd_IFEM = 851, ibMthd_FEIBs = 852

!     IB coupling: explicit/implicit
      INTEGER(KIND=IKIND), PARAMETER :: ibCpld_NA = 899, ibCpld_E = 898,
     2   ibCpld_I = 897

!     IB interpolation: direct extraction / L2 projection
      INTEGER(KIND=IKIND), PARAMETER :: ibIntrp_NA = 900,
     2   ibIntrp_DI = 901, ibIntrp_L2 = 902
!--------------------------------------------------------------------
!#######################################################################

ext: f

max_stars_repo_name: ktbolt/svFSI

max_issues_repo_name: ktbolt/svFSI

max_forks_repo_name: ktbolt/svFSI


program one_plus_one

use add_dummy, only : add_one

implicit none (type, external)

if (add_one(1) /= 2) error stop "1+1 = 2"

end program


text:
      SUBROUTINE scalfor(gcx, axm, bxm, axd, bxd, cx, iflag)
      USE vmec_main
      USE vmec_params
      USE vmec_dim, ONLY: ns, nrzt
      USE realspace, ONLY: wint, ru0
      IMPLICIT NONE
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(in) :: iflag
      REAL(rprec), DIMENSION(ns,0:ntor,0:mpol1,ntmax),
     1  INTENT(inout) :: gcx
      REAL(rprec), DIMENSION(ns+1,2), INTENT(in) ::
     1  axm, bxm, axd, bxd
      REAL(rprec), DIMENSION(ns), INTENT(in) :: cx
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      REAL(rprec), PARAMETER :: ftol_edge = 1.e-9_dp, c1p5 = 1.5_dp,
     1      fac = 0.25_dp, edge_pedestal = 0.05_dp
      INTEGER :: m , mp, n, js, jmax, jmin4(0:mnsize-1)
      REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE :: ax, bx, dx
      REAL(rprec) :: mult_fac
!      LOGICAL :: ledge
C-----------------------------------------------
      ALLOCATE (ax(ns,0:ntor,0:mpol1), bx(ns,0:ntor,0:mpol1),
     1          dx(ns,0:ntor,0:mpol1))

      jmax = ns
      IF (ivac .lt. 1) jmax = ns1

!     FOR SOME 3D PLASMAS, THIS SOMETIME HELPS (CHOOSE mult_fac =1 otherwise)
!     TO AVOID JACOBIAN RESETS BY GIVING A SMOOTH TRANSITION FROM FIXED TO FREE ITERATIONS
!      mult_fac = 1._dp/(1._dp + 10*(fsqr+fsqz))
!      gcx(ns,:,:,:) = mult_fac*gcx(ns,:,:,:)

      DO m = 0, mpol1
         mp = MOD(m,2) + 1
         DO n = 0, ntor
            DO js = jmin2(m), jmax
               ax(js,n,m) = -(axm(js+1,mp) + bxm(js+1,mp)*m**2)
               bx(js,n,m) = -(axm(js,mp) + bxm(js,mp)*m**2)
               dx(js,n,m) = -(axd(js,mp) + bxd(js,mp)*m**2
     1                    + cx(js)*(n*nfp)**2)
            END DO

            IF (m .eq. 1) THEN
               dx(2,n,m) = dx(2,n,m) + bx(2,n,m)
               DO js = jmin2(m), jmax
                   ax(js,n,m) = c1p5*ax(js,n,m)
                   bx(js,n,m) = c1p5*bx(js,n,m)
                   dx(js,n,m) = c1p5*dx(js,n,m)
               END DO
            END IF
         END DO
      END DO

      IF (jmax .ge. ns) THEN
!
!     SMALL EDGE PEDESTAL NEEDED TO IMPROVE CONVERGENCE
!     IN PARTICULAR, NEEDED TO ACCOUNT FOR POTENTIAL ZERO
!     EIGENVALUE DUE TO NEUMANN (GRADIENT) CONDITION AT EDGE
!
!!         dx(ns,0,0) = (1+edge_pedestal)*dx(ns,0,0)
!!           dx(ns,0,0) = dx(ns,0,0) + alfR*ohs**2
         dx(ns,:,0:1)     = (1+edge_pedestal)  *dx(ns,:,0:1)
         dx(ns,:,2:mpol1) = (1+2*edge_pedestal)*dx(ns,:,2:mpol1)
!
!     STABILIZATION ALGORITHM FOR ZC_00(NS)
!     FOR UNSTABLE CASE, HAVE TO FLIP SIGN OF -FAC -> +FAC FOR CONVERGENCE
!     COEFFICIENT OF < Ru (R Pvac)> ~ -fac*(z-zeq) WHERE fac (EIGENVALUE, OR
!     FIELD INDEX) DEPENDS ON THE EQUILIBRIUM MAGNETIC FIELD AND CURRENT,
!     AND zeq IS THE EQUILIBRIUM EDGE VALUE OF Z00
          mult_fac = MIN(fac, fac*hs*15)
          IF (iflag .eq. 1) THEN
!
!     METHOD 1: SUBTRACT (INSTABILITY) Pedge ~ fac*z/hs FROM PRECONDITIONER AT EDGE
!
             dx(ns,0,0) = dx(ns,0,0)*(1-mult_fac)/(1+edge_pedestal)
          END IF

      ENDIF

!
!     ACCELERATE (IMPROVE) CONVERGENCE OF FREE BOUNDARY. THIS WAS ADDED
!     TO DEAL WITH CASES WHICH MIGHT OTHERWISE DIVERGE. BY DECREASING THE
!     FSQ TOLERANCE LEVEL WHERE THIS KICKS IN (FTOL_EDGE), THE USER CAN
!     TURN-OFF THIS FEATURE
!
!     DIAGONALIZE (DX DOMINANT) AND REDUCE FORCE (DX ENHANCED) AT EDGE
!     TO IMPROVE CONVERGENCE FOR N != 0 TERMS
!

!      ledge = .false.
!      IF ((fsqr+fsqz) .lt. ftol_edge) ledge = .true.
!      IF ((iter2-iter1).lt.400 .or. ivac.lt.1) ledge = .false.

!      IF (ledge) THEN
!         dx(ns,1:,1:) = 3*dx(ns,1:,1:)
!      END IF

!     FOR DATA MATCHING MODE (0 <= IRESIDUE < 3),
!     MAGNETIC AXIS IS FIXED SO JMIN3(0) => 2 FOR M=0,N=0

      jmin4 = jmin3
      IF (iresidue.ge.0 .and. iresidue.lt.3) jmin4(0) = 2

      CALL tridslv (ax, dx, bx, gcx, jmin4, jmax, mnsize-1, ns, ntmax)

      DEALLOCATE (ax, bx, dx)

      END SUBROUTINE scalfor


text:
*dk,chklun
      subroutine chklun(ifile, iflaga, ignore)
C
C#######################################################################
C
C      PURPOSE -
C
C         THIS ROUTINE CHECK TO SEE IF A FILE HAS BEEN OPENED AND
C              ASSIGN A LOGICAL UNIT NUMBER.
C
C      INPUT ARGUMENTS -
C
C         ifile    - THE FILE NAME.
C         IGNORE   - IGNORE THIS FIELD.
C
C      OUTPUT ARGUMENTS -
C
C         iflaga - INDICATES IF THE FILE (OR UNIT) NUMBER IS OPEN.
C                  = 0 ==> FILE IS NOT OPEN.
C                  = 1 ==> FILE IS OPEN.
C
C      CHANGE HISTORY -
C
C         $Log: chklun.f,v $
C         Revision 2.00  2007/11/03 00:49:10  spchu
C         Import to CVS
C
CPVCS
CPVCS       Rev 1.0   03/17/95 19:33:42   het
CPVCS    Original version converted from the CRAY
C
C#######################################################################
C
      character ifile*(*)
      integer iflaga, ignore
      logical opend

c     The I/O library does not have entry points for handling
c     type INTEGER*8 I/O control list specifiers
      integer*4 iunit

      pointer (ipival, ival(1))
      integer icharlnf
C
C#######################################################################
C
      iflaga=-1
      ierror=-1
      if(len(ifile).eq.0) then
         lenmax=0
C        *** IF THIS HAPPENS, THIS IS "PROABALY" AN INTEGER FIELD WITH A
C               WITH A CHARACTER CONTAINED WITHIN, MAKE AN ASSUMPTION.
      else
         lenmax=icharlnf(ifile)
      endif
      if(lenmax.eq.0) then
         ipival=loc(ifile)
         iunit=ival(1)
         inquire(unit=iunit,opened=opend,err=9999)
      else
         inquire(file=ifile(1:lenmax),opened=opend,err=9999)
      endif
      if(opend.eqv..false.) then
         iflaga=0
      else
         iflaga=1
      endif
      goto 9999
 9999 continue
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
module bc_state_artificial_viscosity_symmetry
#include <messenger.h>
    use mod_kinds,              only: rk,ik
    use mod_constants,          only: ONE, TWO, TEN
    use type_bc_state,          only: bc_state_t
    use type_chidg_worker,      only: chidg_worker_t
    use type_properties,        only: properties_t
    use type_point,             only: point_t
    use mpi_f08,                only: mpi_comm
    use DNAD_D
    implicit none

    !> Extrapolation boundary condition
    !!      - Extrapolate interior variables to be used for calculating the boundary flux.
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/8/2016
    !!
    !-------------------------------------------------------------------------------------------
    type, public, extends(bc_state_t) :: artificial_viscosity_symmetry_t

    contains

        procedure   :: init
        procedure   :: compute_bc_state

    end type artificial_viscosity_symmetry_t
    !*******************************************************************************************

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   8/29/2016
    !!
    !--------------------------------------------------------------------------------
    subroutine init(self)
        class(artificial_viscosity_symmetry_t),   intent(inout) :: self

        !
        ! Set name, family
        !
        call self%set_name('Artificial Viscosity Symmetry')
        call self%set_family('Symmetry')

    end subroutine init
    !********************************************************************************

    !>
    !!
    !!  @author Nathan A. Wukie
    !!  @date   9/8/2016
    !!
    !!
    !-----------------------------------------------------------------------------------------
    subroutine compute_bc_state(self,worker,prop,bc_COMM)
        class(artificial_viscosity_symmetry_t), intent(inout)   :: self
        type(chidg_worker_t),                   intent(inout)   :: worker
        class(properties_t),                    intent(inout)   :: prop
        type(mpi_comm),                         intent(in)      :: bc_COMM

        ! Storage at quadrature nodes
        type(AD_D), allocatable, dimension(:)   ::      &
            eps_m,  deps_dx_m,  deps_dy_m,  deps_dz_m

        !
        ! Interpolate interior solution to quadrature nodes
        !
        eps_m      = worker%get_field('Artificial Viscosity' , 'value', 'face interior')
        deps_dx_m  = worker%get_field('Artificial Viscosity' , 'grad1', 'face interior')
        deps_dy_m  = worker%get_field('Artificial Viscosity' , 'grad2', 'face interior')
        deps_dz_m  = worker%get_field('Artificial Viscosity' , 'grad3', 'face interior')

        !
        ! Store computed boundary state
        !
        call worker%store_bc_state('Artificial Viscosity' , eps_m,     'value')
        call worker%store_bc_state('Artificial Viscosity' , deps_dx_m, 'grad1')
        call worker%store_bc_state('Artificial Viscosity' , deps_dy_m, 'grad2')
        call worker%store_bc_state('Artificial Viscosity' , deps_dz_m, 'grad3')

    end subroutine compute_bc_state
    !******************************************************************************************

end module bc_state_artificial_viscosity_symmetry

max_stars_repo_name: wanglican/ChiDG

max_issues_repo_name: haohb/ChiDG

max_forks_repo_name: haohb/ChiDG