--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/07/17 19:25:51 631 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2003/07/23 22:13:59 648 @@ -4,7 +4,7 @@ module eam use force_globals use status use atype_module -#ifdef MPI +#ifdef IS_MPI use mpiSimulation #endif implicit none @@ -14,12 +14,17 @@ module eam 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 @@ -49,18 +54,19 @@ module eam real( kind = dp), dimension(:), allocatable :: rho real( kind = dp), dimension(:), allocatable :: dfrhodrho -! real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho + real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho !! Arrays for MPI storage -#ifdef MPI +#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 @@ -78,9 +84,10 @@ module eam public :: init_EAM_FF ! public :: EAM_new_rcut -! public :: do_EAM_pair + public :: do_eam_pair public :: newEAMtype - + public :: calc_eam_prepair_rho + public :: calc_eam_preforce_Frho contains @@ -106,6 +113,7 @@ contains integer :: alloc_stat integer :: current integer,pointer :: Matchlist(:) => null() + type (EAMtype), pointer :: makeEamtype => null() status = 0 !! Assume that atypes has already been set and get the total number of types in atypes @@ -123,13 +131,13 @@ contains current = EAMList%currentAddition - !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat) + call allocate_EAMType(eam_nrho,eam_nr,makeEamtype,stat=alloc_stat) if (alloc_stat /= 0) then status = -1 return end if + makeEamtype => EAMList%EAMParams(current) - ! this is a possible bug, we assume a correspondence between the vector atypes and ! EAMAtypes @@ -154,6 +162,8 @@ contains integer :: alloc_stat integer :: number_r, number_rho + + do i = 1, EAMList%currentAddition EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = & @@ -198,12 +208,13 @@ contains enddo current_rcut_max = EAMList%EAMParams(1)%eam_rcut - !! find the smallest 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 @@ -224,7 +235,7 @@ contains integer, intent(out) :: status integer :: nlocal -#ifdef MPI +#ifdef IS_MPI integer :: nrow integer :: ncol #endif @@ -233,7 +244,7 @@ contains nlocal = getNlocal() -#ifdef MPI +#ifdef IS_MPI nrow = getNrow(plan_row) ncol = getNcol(plan_col) #endif @@ -250,14 +261,22 @@ contains 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 MPI +#ifdef IS_MPI if (allocated(frho_row)) deallocate(frho_row) allocate(frho_row(nrow),stat=alloc_stat) @@ -277,7 +296,14 @@ contains 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) @@ -298,7 +324,13 @@ contains 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 @@ -306,12 +338,12 @@ contains subroutine clean_EAM() -! clean non-MPI first +! clean non-IS_MPI first frho = 0.0_dp rho = 0.0_dp dfrhodrho = 0.0_dp ! clean MPI if needed -#ifdef MPI +#ifdef IS_MPI frho_row = 0.0_dp frho_col = 0.0_dp rho_row = 0.0_dp @@ -421,38 +453,78 @@ contains 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. + - call calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1) +#ifdef IS_MPI + myid_atom1 = atid_Row(atom1) + myid_atom2 = atid_Col(atom2) +#else + myid_atom1 = atid(atom1) + myid_atom2 = atid(atom2) +#endif -#ifdef MPI - rho_col(atom2) = rho_col(atom2) + rho_i_at_j + 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 + rho(atom2) = rho(atom2) + rho_i_at_j #endif + endif - call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2) + 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 MPI - rho_row(atom1) = rho_row(atom1) + rho_j_at_i + + + +#ifdef IS_MPI + rho_row(atom1) = rho_row(atom1) + rho_j_at_i #else - rho(atom1) = rho(atom1) + rho_j_at_i + rho(atom1) = rho(atom1) + rho_j_at_i #endif + endif + end subroutine calc_eam_prepair_rho - !! Calculate the functional F(rho) for all atoms - subroutine calc_eam_prepair_Frho(nlocal,pot) + + + + !! 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 in pre-pair -#ifdef MPI +!! 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" @@ -465,16 +537,40 @@ contains 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 +!! 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") @@ -492,28 +588,33 @@ contains 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 - end subroutine calc_eam_prepair_Frho - - subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress) -!Arguments + !! 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 @@ -524,8 +625,10 @@ contains real( kind = dp ) :: d2rhojdrdr real( kind = dp ) :: Fx,Fy,Fz real( kind = dp ) :: r,d2pha,phb,d2phb + integer :: id1,id2 - integer :: atype1,atype2 + integer :: mytype_atom1 + integer :: mytype_atom2 !Local Variables @@ -535,35 +638,64 @@ contains 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) + mytype_atom1 = atid_row(atom1) #else - atype1 = atid(atom1) + mytype_atom1 = 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)) + + 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 - atype2 = atid_col(atom2) + mytype_atom2 = atid_col(atom2) #else - atype2 = atid(atom2) + mytype_atom2 = 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 + ! 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))) @@ -574,7 +706,7 @@ contains endif - if (r.lt.rcj) then + 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))) @@ -591,7 +723,7 @@ contains d2rhojdrdr = d2rhb -#ifdef MPI +#ifdef IS_MPI dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) & + dvpdr @@ -606,7 +738,7 @@ contains fz = dudr * drdz -#ifdef MPI +#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 @@ -615,157 +747,114 @@ contains 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 MPI + + + 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) + 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) + + + + 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 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 + + 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, dy, d2y + 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(dble(nx-1) * (x - xa(1)) / (xa(nx) - xa(1))) + 1 + 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(default_error,*) "EAM_splint: x is outside bounds of spline" + 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(default_error,*) "EAM_splint: x is outside bounding range" + write(errMSG,*) "EAM_splint:",x," x is outside bounding range" + call handleError(routineName,errMSG) endif endif @@ -778,20 +867,21 @@ end subroutine calc_eam_pair 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) + + 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 @@ -805,9 +895,13 @@ end subroutine calc_eam_pair real( kind = DP ), dimension(:) :: yppa real( kind = DP ), dimension(size(xa)) :: u real( kind = DP ) :: yp1,ypn,un,qn,sig,p - character boundary + 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