--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/04/14 18:19:15 492 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2004/05/07 21:35:05 1150 @@ -1,292 +1,747 @@ -module calc_eam - use definitions, ONLY : DP +module eam + use definitions, ONLY : DP,default_error + use simulation use force_globals -#ifdef MPI + 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 + logical, save :: haveRcut = .false. + character(len = statusMsgSize) :: errMesg + integer :: eam_err - !! standard eam stuff - integer :: n_eam_atypes - integer, allocatable, dimension(:) :: eam_atype - real( kind = DP ), allocatable, dimension(:) :: eam_dr - integer, allocatable, dimension(:) :: eam_nr - integer, allocatable, dimension(:) :: eam_nrho - real( kind = DP ), allocatable, dimension(:) :: eam_lattice - real( kind = DP ), allocatable, dimension(:) :: eam_drho - integer , allocatable, dimension(:) :: eam_atype_map - real( kind = DP ), allocatable, dimension(:,:) :: eam_rvals - real( kind = DP ), allocatable, dimension(:,:) :: eam_rhovals - real( kind = DP ), allocatable, dimension(:) :: eam_rcut - real( kind = DP ), allocatable, dimension(:,:) :: eam_F_rho - real( kind = DP ), allocatable, dimension(:,:) :: eam_Z_r - real( kind = DP ), allocatable, dimension(:,:) :: eam_rho_r - real( kind = DP ), allocatable, dimension(:,:) :: eam_phi_r - real( kind = DP ), allocatable, dimension(:,:) :: eam_F_rho_pp - real( kind = DP ), allocatable, dimension(:,:) :: eam_Z_r_pp - real( kind = DP ), allocatable, dimension(:,:) :: eam_rho_r_pp - real( kind = DP ), allocatable, dimension(:,:) :: eam_phi_r_pp + 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 - real( kind = DP ), private :: time0,time1,time2,time3 - integer, private :: eam_err + !! Arrays for derivatives used in force calculation + real( kind = dp), dimension(:), allocatable :: frho + real( kind = dp), dimension(:), allocatable :: rho - private :: mass_weight - private :: allocate_eam_atype,allocate_eam_module,deallocate_eam_module - private :: read_eam_pot, get_eam_sizes + real( kind = dp), dimension(:), allocatable :: dfrhodrho + real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho -contains - subroutine initialize_eam() +!! Arrays for MPI storage +#ifdef IS_MPI + real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_col + real( kind = dp),save, dimension(:), allocatable :: dfrhodrho_row + real( kind = dp),save, dimension(:), allocatable :: frho_row + real( kind = dp),save, dimension(:), allocatable :: frho_col + real( kind = dp),save, dimension(:), allocatable :: rho_row + real( kind = dp),save, dimension(:), allocatable :: rho_col + real( kind = dp),save, dimension(:), allocatable :: rho_tmp + real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_col + real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row +#endif + type, private :: EAMTypeList + integer :: n_eam_types = 0 + integer :: currentAddition = 0 + + type (EAMtype), pointer :: EAMParams(:) => null() + end type EAMTypeList - character(len=80) :: eam_pot_file - integer :: i, j, max_size, prev_max_size - integer :: number_rho, number_r - integer :: eam_unit - integer :: this_error - character(len=300) :: msg - integer, external :: nfiles - !for mpi + type (eamTypeList), save :: EAMList + !! standard eam stuff -#ifdef MPI - if (node == 0) & - n_eam_atypes = nfiles(trim(eam_pot_dir)//char(0)) - call mpi_bcast(n_eam_atypes,1,mpi_integer,0,mpi_comm_world,mpi_err) - if (n_eam_atypes == -1) then - call error("INITIALIZE_EAM","NO EAM potentials found!") - endif - write(msg,'(a5,i4,a12,i5,a14)') 'Node: ',node,' reading ...', & - n_eam_atypes, ' eam atom types' - call info('INITIALIZE_EAM', trim(msg)) -#else - n_eam_atypes = nfiles(trim(eam_pot_dir)//char(0)) - if (n_eam_atypes == -1) then - call error("INITIALIZE_EAM","NO EAM potentials found!") - endif + public :: init_EAM_FF + public :: setCutoffEAM + public :: do_eam_pair + public :: newEAMtype + public :: calc_eam_prepair_rho + public :: calc_eam_preforce_Frho + public :: clean_EAM - write(msg,'(a12,i5,a14)') ' Reading ...', & - n_eam_atypes, ' eam atom types' - call info('INITIALIZE_EAM', trim(msg)) -#endif +contains - call allocate_eam_atype(n_eam_atypes) + 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 - !! get largest number of data points for any potential -#ifdef MPI - if (node == 0) then -#endif - prev_max_size = 0 - do i = 1, n_eam_atypes - call getfilename(i, eam_pot_file) - max_size = max(get_eam_sizes( & - trim(eam_pot_dir) // '/' // eam_pot_file), & - prev_max_size) - prev_max_size = max_size - end do -#ifdef MPI + + !! 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 mpi_bcast(max_size,1,mpi_integer,0,mpi_comm_world,mpi_err) -#endif + call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat) + if (alloc_stat /= 0) then + status = -1 + return + end if - call allocate_eam_module(n_eam_atypes,max_size) - allocate(eam_atype_map(get_max_atype())) + ! 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 -#ifdef MPI - if (node == 0) then -#endif - do i = 1, n_eam_atypes - call getfilename(i, eam_pot_file) - call read_eam_pot(i,trim(eam_pot_dir) // '/' // eam_pot_file, & - this_error) + end subroutine newEAMtype - do j = 1, eam_nr(i) - eam_rvals(j,i) = dble(j-1)*eam_dr(i) - enddo - do j = 1, eam_nrho(i) - eam_rhovals(j,i) = dble(j-1)*eam_drho(i) - enddo + subroutine init_EAM_FF(status) + integer :: status + integer :: i,j + real(kind=dp) :: current_rcut_max + integer :: alloc_stat + integer :: number_r, number_rho + + + status = 0 + 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: - do j = 1, eam_nrho(i) - eam_F_rho(j,i) = eam_F_rho(j,i)*23.06054E0_DP - enddo + 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: - eam_phi_r(1,i) = 0.0E0_DP - do j = 2, eam_nr(i) - eam_phi_r(j,i) = (eam_Z_r(j,i)**2)/eam_rvals(j,i) - eam_phi_r(j,i) = eam_phi_r(j,i)*331.999296E0_DP + 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 -#ifdef MPI - call info('INITIALIZE_EAM','NODE 0: Distributing spline arrays') - endif + - call mpi_bcast(this_error,n_eam_atypes,mpi_integer,0, & - mpi_comm_world,mpi_err) - if (this_error /= 0) then - call error('INITIALIZE_EAM',"Cannot read eam files") - endif + 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 - call mpi_bcast(eam_atype,n_eam_atypes,mpi_integer,0, & - mpi_comm_world,mpi_err) +! current_rcut_max = EAMList%EAMParams(1)%eam_rcut + !! find the smallest rcut for any eam atype +! do i = 2, EAMList%currentAddition +! current_rcut_max =max(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) +! end do - !! distribute values to cluster...... - call mpi_bcast(eam_nr,n_eam_atypes,mpi_integer,& - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_nrho,n_eam_atypes,mpi_integer,& - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_rvals,n_eam_atypes*max_size,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_rcut,n_eam_atypes,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_rhovals,n_eam_atypes*max_size,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) +! 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 + write(*,*) "allocateEAM failed" + status = -1 + return + endif + + end subroutine init_EAM_FF - !! distribute arrays - call mpi_bcast(eam_rho_r,n_eam_atypes*max_size,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_Z_r,n_eam_atypes*max_size,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_F_rho,n_eam_atypes*max_size,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) - call mpi_bcast(eam_phi_r,n_eam_atypes*max_size,mpi_double_precision, & - 0,mpi_comm_world,mpi_err) +!! 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 +#ifdef IS_MPI + integer :: nrow + integer :: ncol #endif - call info('INITIALIZE_EAM', 'creating splines') + integer :: alloc_stat - do i = 1, n_eam_atypes - number_r = eam_nr(i) - number_rho = eam_nrho(i) - call eam_spline(i, number_r, eam_rvals, eam_rho_r, eam_rho_r_pp, & - 0.0E0_DP, 0.0E0_DP, 'N') - call eam_spline(i, number_r, eam_rvals, eam_Z_r, eam_Z_r_pp, & - 0.0E0_DP, 0.0E0_DP, 'N') - call eam_spline(i, number_rho, eam_rhovals, eam_F_rho, eam_F_rho_pp, & - 0.0E0_DP, 0.0E0_DP, 'N') - call eam_spline(i, number_r, eam_rvals, eam_phi_r, eam_phi_r_pp, & - 0.0E0_DP, 0.0E0_DP, 'N') - enddo + status = 0 +#ifdef IS_MPI + nrow = getNrow(plan_row) + ncol = getNcol(plan_col) +#endif - do i = 1, n_eam_atypes - eam_atype_map(eam_atype(i)) = i - end do + 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 - call info('INITIALIZE_EAM','Done creating splines') + if (allocated(rho_tmp)) deallocate(rho_tmp) + allocate(rho_tmp(nlocal),stat=alloc_stat) + if (alloc_stat /= 0) then + status = -1 + return + end if - return - end subroutine initialize_eam + 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 - subroutine allocate_eam_atype(n_size_atype) - integer, intent(in) :: n_size_atype - allocate(eam_atype(n_size_atype)) - allocate(eam_drho(n_size_atype)) - allocate(eam_dr(n_size_atype)) - allocate(eam_nr(n_size_atype)) - allocate(eam_nrho(n_size_atype)) - allocate(eam_lattice(n_size_atype)) - allocate(eam_rcut(n_size_atype)) +! Now do column arrays - end subroutine allocate_eam_atype + 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 - subroutine allocate_eam_module(n_size_atype,n_eam_points) - integer, intent(in) :: n_eam_points - integer, intent(in) :: n_size_atype + end subroutine allocateEAM - allocate(eam_rvals(n_eam_points,n_size_atype)) - allocate(eam_rhovals(n_eam_points,n_size_atype)) - allocate(eam_F_rho(n_eam_points,n_size_atype)) - allocate(eam_Z_r(n_eam_points,n_size_atype)) - allocate(eam_rho_r(n_eam_points,n_size_atype)) - allocate(eam_phi_r(n_eam_points,n_size_atype)) - allocate(eam_F_rho_pp(n_eam_points,n_size_atype)) - allocate(eam_Z_r_pp(n_eam_points,n_size_atype)) - allocate(eam_rho_r_pp(n_eam_points,n_size_atype)) - allocate(eam_phi_r_pp(n_eam_points,n_size_atype)) +!! C sets rcut to be the largest cutoff of any atype +!! present in this simulation. Doesn't include all atypes +!! sim knows about, just those in the simulation. + subroutine setCutoffEAM(rcut, status) + real(kind=dp) :: rcut + integer :: status + status = 0 - end subroutine allocate_eam_module + EAM_rcut = rcut - subroutine deallocate_eam_module() + end subroutine setCutoffEAM - deallocate(eam_atype) - deallocate(eam_drho) - deallocate(eam_dr) - deallocate(eam_nr) - deallocate(eam_nrho) - deallocate(eam_lattice) - deallocate(eam_atype_map) - deallocate(eam_rvals) - deallocate(eam_rhovals) - deallocate(eam_rcut) - deallocate(eam_Z_r) - deallocate(eam_rho_r) - deallocate(eam_phi_r) - deallocate(eam_F_rho_pp) - deallocate(eam_Z_r_pp) - deallocate(eam_rho_r_pp) - deallocate(eam_phi_r_pp) - end subroutine deallocate_eam_module + + 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 + rho_tmp = 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 + + + + +#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 +! write(*,*) atom1,atom2,r,rho_i_at_j + 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 + + 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_tmp,plan_col,eam_err) + if (eam_err /= 0 ) then + write(errMsg,*) " Error scattering rho_col into rho" + call handleError(RoutineName,errMesg) + endif - subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress) -!Arguments + rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal) +#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(atom) = u + dfrhodrho(atom) = u1 + d2frhodrhodrho(atom) = 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, sw, vpair, 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 ) :: pot, sw, vpair + real( kind = dp ), dimension(3,nLocal) :: 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 + ! write(*,*) "Frho: ", Frho(atom1) + 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 + - r = dsqrt(rijsq) - efr(1,j) = -rxij - efr(2,j) = -ryij - efr(3,j) = -rzij + 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) - call calc_eam_rho(r, rha, drha, d2rha, atype1) - call calc_eam_phi(r, pha, dpha, d2pha, atype1) - rci = eam_rcut(eam_atype_map(atype1)) -#ifdef MPI - atype2 = ident_col(j) + +! get cutoff for atom 1 + rci = EAMList%EAMParams(mytype_atom1)%eam_rcut +#ifdef IS_MPI + mytype_atom2 = atid_col(atom2) #else - atype2 = ident(j) + 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)) + ! 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) - phab = 0.0E0_DP - dvpdr = 0.0E0_DP - d2vpdrdr = 0.0E0_DP + !! 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) - if (r.lt.rci) then + +! 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))) @@ -295,9 +750,9 @@ contains pha*((d2rhb/rha) - 2.0E0_DP*(drhb*drha/rha/rha) + & (2.0E0_DP*rhb*drha*drha/rha/rha/rha) - (rhb*d2rha/rha/rha))) endif + - - if (r.lt.rcj) then + 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))) @@ -306,294 +761,233 @@ contains phb*((d2rha/rhb) - 2.0E0_DP*(drha*drhb/rhb/rhb) + & (2.0E0_DP*rha*drhb*drhb/rhb/rhb/rhb) - (rha*d2rhb/rhb/rhb))) endif - - -#ifdef MPI - - e_row(i) = e_row(i) + phab*0.5 - e_col(i) = e_col(i) + phab*0.5 -#else - if (do_pot) pot = pot + phab -#endif - + drhoidr = drha drhojdr = drhb d2rhoidrdr = d2rha d2rhojdrdr = d2rhb -#ifdef MPI - dudr = drhojdr*dfrhodrho_row(i)+drhoidr*dfrhodrho_col(j) & + + +#ifdef IS_MPI + dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) & + dvpdr - if (nmflag) then - d2 = d2vpdrdr + & - d2rhoidrdr*dfrhodrho_col(j) + & - d2rhojdrdr*dfrhodrho_row(i) + & - drhoidr*drhoidr*d2frhodrhodrho_col(j) + & - drhojdr*drhojdr*d2frhodrhodrho_row(i) - endif #else - dudr = drhojdr*dfrhodrho(i)+drhoidr*dfrhodrho(j) & + dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) & + dvpdr - - d2 = d2vpdrdr + & - d2rhoidrdr*dfrhodrho(j) + & - d2rhojdrdr*dfrhodrho(i) + & - drhoidr*drhoidr*d2frhodrhodrho(j) + & - drhojdr*drhojdr*d2frhodrhodrho(i) + ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2) #endif + fx = dudr * drdx + fy = dudr * drdy + fz = dudr * drdz - do dim = 1, 3 - drdx1 = efr(dim,j) / r - ftmp = dudr * drdx1 +#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 + vpair = vpair + phab -#ifdef MPI - f_col(dim,j) = f_col(dim,j) - ftmp - f_row(dim,i) = f_row(dim,i) + ftmp + 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 - f(dim,j) = f(dim,j) - ftmp - f(dim,i) = f(dim,i) + ftmp -#endif - if (nmflag) then - idim = 3 * (i-1) + dim - jdim = 3 * (j-1) + dim + if(do_pot) then + pot = pot + phab + end if + vpair = vpair + phab - do dim2 = 1, 3 + 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 - kt1 = d2 * efr(dim,j) * efr(dim2,j)/r/r - kt2 = - dudr * efr(dim,j) * efr(dim2,j)/r/r/r + drhoidr = drha + drhojdr = drhb + d2rhoidrdr = d2rha + d2rhojdrdr = d2rhb - if (dim.eq.dim2) then - kt3 = dudr / r - else - kt3 = 0.0E0_DP - endif +#ifdef IS_MPI + d2 = d2vpdrdr + & + d2rhoidrdr*dfrhodrho_col(atom2) + & + d2rhojdrdr*dfrhodrho_row(atom1) + & + drhoidr*drhoidr*d2frhodrhodrho_col(atom2) + & + drhojdr*drhojdr*d2frhodrhodrho_row(atom1) + +#else - ! The factor of 2 below is to compensate for - ! overcounting. - ! Mass weighting is done separately... + d2 = d2vpdrdr + & + d2rhoidrdr*dfrhodrho(atom2) + & + d2rhojdrdr*dfrhodrho(atom1) + & + drhoidr*drhoidr*d2frhodrhodrho(atom2) + & + drhojdr*drhojdr*d2frhodrhodrho(atom1) +#endif + end if - ktmp = (kt1+kt2+kt3)/2.0E0_DP - idim2 = 3 * (i-1) + dim2 - jdim2 = 3 * (j-1) + dim2 - d(idim, idim2) = d(idim,idim2) + ktmp - d(idim2, idim) = d(idim2,idim) + ktmp - - d(idim, jdim2) = d(idim,jdim2) - ktmp - d(idim2, jdim) = d(idim2,jdim) - ktmp - - d(jdim, idim2) = d(jdim,idim2) - ktmp - d(jdim2, idim) = d(jdim2,idim) - ktmp - - d(jdim, jdim2) = d(jdim,jdim2) + ktmp - d(jdim2, jdim) = d(jdim2,jdim) + ktmp - - enddo - endif - enddo - - endif - enddo -endif - -enddo + + + 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 calc_eam_pair + + end subroutine do_eam_pair -subroutine calc_eam_rho(r, rho, drho, d2rho, atype) - ! include 'headers/sizes.h' + 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 -integer atype, etype, number_r -real( kind = DP ) :: r, rho, drho, d2rho -integer :: i + + ! 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 -etype = eam_atype_map(atype) + 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) + -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 + end subroutine eam_splint -return -end subroutine calc_eam_rho -subroutine calc_eam_frho(dens, u, u1, u2, atype) + subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary) - ! include 'headers/sizes.h' -integer atype, etype, number_rho -real( kind = DP ) :: dens, u, u1, u2 -real( kind = DP ) :: rho_vals + ! 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 -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 + 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 -return -end subroutine calc_eam_frho + 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 -subroutine calc_eam_phi(r, phi, dphi, d2phi, atype) + end subroutine eam_spline -integer atype, etype, number_r -real( kind = DP ) :: r, phi, dphi, d2phi - -etype = eam_atype_map(atype) - -if (r.lt.eam_rcut(etype)) then -number_r = eam_nr(etype) -call eam_splint(etype, number_r, eam_rvals, eam_phi_r, & - eam_phi_r_pp, r, phi, dphi, d2phi) -else -phi = 0.0E0_DP -dphi = 0.0E0_DP -d2phi = 0.0E0_DP -endif - -return -end subroutine calc_eam_phi - - -subroutine eam_splint(atype, nx, xa, ya, yppa, x, y, dy, d2y) - - ! include 'headers/sizes.h' - -real( kind = DP ), dimension(:,:) :: xa -real( kind = DP ), dimension(:,:) :: ya -real( kind = DP ), dimension(:,:) :: yppa -real( kind = DP ) :: x, y, dy, d2y -real( kind = DP ) :: del, h, a, b, c, d - - -integer atype, nx, j - - -! 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,atype)) / (xa(nx,atype) - xa(1,atype))) + 1 - -! check to make sure we're inside the spline range: -if ((j.gt.nx).or.(j.lt.1)) call error('eam_splint', & -'x is outside bounds of spline') - -! check to make sure we haven't screwed up the calculation of j: -if ((x.lt.xa(j,atype)).or.(x.gt.xa(j+1,atype))) then -if (j.ne.nx) then - call error('eam_splint', & - 'x is outside bounding range') -endif -endif - -del = xa(j+1,atype) - x -h = xa(j+1,atype) - xa(j,atype) - -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,atype) + b*ya(j+1,atype) + c*yppa(j,atype) + d*yppa(j+1,atype) - -dy = (ya(j+1,atype)-ya(j,atype))/h & -- (3.0E0_DP*a*a - 1.0E0_DP)*h*yppa(j,atype)/6.0E0_DP & -+ (3.0E0_DP*b*b - 1.0E0_DP)*h*yppa(j+1,atype)/6.0E0_DP - -d2y = a*yppa(j,atype) + b*yppa(j+1,atype) - -return -end subroutine eam_splint - -subroutine eam_spline(atype, nx, xa, ya, yppa, yp1, ypn, boundary) - - ! include 'headers/sizes.h' - - - ! 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, atype, max_array_size - -real( kind = DP ), dimension(:,:) :: xa -real( kind = DP ), dimension(:,:) :: ya -real( kind = DP ), dimension(:,:) :: yppa -real( kind = DP ), allocatable, dimension(:) :: u -real( kind = DP ) :: yp1,ypn,un,qn,sig,p -character boundary - -max_array_size = size(xa,1) -allocate(u(max_array_size)) - - -if ((boundary.eq.'l').or.(boundary.eq.'L').or. & -(boundary.eq.'b').or.(boundary.eq.'B')) then -yppa(1, atype) = -0.5E0_DP -u(1) = (3.0E0_DP/(xa(2,atype)-xa(1,atype)))*((ya(2,atype)-& - ya(1,atype))/(xa(2,atype)-xa(1,atype))-yp1) -else -yppa(1,atype) = 0.0E0_DP -u(1) = 0.0E0_DP -endif - -do i = 2, nx - 1 -sig = (xa(i,atype) - xa(i-1,atype)) / (xa(i+1,atype) - xa(i-1,atype)) -p = sig * yppa(i-1,atype) + 2.0E0_DP -yppa(i,atype) = (sig - 1.0E0_DP) / p -u(i) = (6.0E0_DP*((ya(i+1,atype)-ya(i,atype))/(xa(i+1,atype)-xa(i,atype)) - & - (ya(i,atype)-ya(i-1,atype))/(xa(i,atype)-xa(i-1,atype)))/ & - (xa(i+1,atype)-xa(i-1,atype)) - 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,atype)-xa(nx-1,atype)))* & - (ypn-(ya(nx,atype)-ya(nx-1,atype))/(xa(nx,atype)-xa(nx-1,atype))) -else -qn = 0.0E0_DP -un = 0.0E0_DP -endif - -yppa(nx,atype)=(un-qn*u(nx-1))/(qn*yppa(nx-1,atype)+1.0E0_DP) - -do k = nx-1, 1, -1 -yppa(k,atype)=yppa(k,atype)*yppa(k+1,atype)+u(k) -enddo - -deallocate(u) -return -end subroutine eam_spline - - - - -end module calc_eam +end module eam