--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/04/14 21:16:37 497 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2003/07/17 19:25:51 631 @@ -1,257 +1,509 @@ -module calc_eam - use definitions, ONLY : DP +module eam + use definitions, ONLY : DP,default_error + use simulation use force_globals + use status + use atype_module #ifdef MPI use mpiSimulation #endif + implicit none PRIVATE logical, save :: EAM_FF_initialized = .false. integer, save :: EAM_Mixing_Policy - integer, save :: EAM_rcut + real(kind = dp), save :: EAM_rcut + character(len = 200) :: errMsg + character(len=*), parameter :: RoutineName = "EAM MODULE" + logical :: cleanme = .true. + + type, private :: EAMtype + integer :: eam_atype + real( kind = DP ) :: eam_dr + integer :: eam_nr + integer :: eam_nrho + real( kind = DP ) :: eam_lattice + real( kind = DP ) :: eam_drho + real( kind = DP ) :: eam_rcut + integer :: eam_atype_map + + real( kind = DP ), pointer, dimension(:) :: eam_rvals => null() + real( kind = DP ), pointer, dimension(:) :: eam_rhovals => null() + real( kind = DP ), pointer, dimension(:) :: eam_F_rho => null() + real( kind = DP ), pointer, dimension(:) :: eam_Z_r => null() + real( kind = DP ), pointer, dimension(:) :: eam_rho_r => null() + real( kind = DP ), pointer, dimension(:) :: eam_phi_r => null() + real( kind = DP ), pointer, dimension(:) :: eam_F_rho_pp => null() + real( kind = DP ), pointer, dimension(:) :: eam_Z_r_pp => null() + real( kind = DP ), pointer, dimension(:) :: eam_rho_r_pp => null() + real( kind = DP ), pointer, dimension(:) :: eam_phi_r_pp => null() + end type EAMtype - !! 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 + !! Arrays for derivatives used in force calculation + real( kind = dp), dimension(:), allocatable :: frho + real( kind = dp), dimension(:), allocatable :: rho - public :: init_EAM_FF - public :: EAM_new_rcut - public :: do_EAM_pair + real( kind = dp), dimension(:), allocatable :: dfrhodrho +! real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho +!! Arrays for MPI storage +#ifdef MPI + real( kind = dp), dimension(:), allocatable :: dfrhodrho_col + real( kind = dp), dimension(:), allocatable :: dfrhodrho_row + real( kind = dp), dimension(:), allocatable :: frho_row + real( kind = dp), dimension(:), allocatable :: frho_col + real( kind = dp), dimension(:), allocatable :: rho_row + real( kind = dp), dimension(:), allocatable :: rho_col -contains - subroutine init_EAM_FF() +#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) :: 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 :: EAM_new_rcut +! public :: do_EAM_pair + public :: newEAMtype - 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())) -#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) + ! this is a possible bug, we assume a correspondence between the vector atypes and + ! EAMAtypes + + EAMList%EAMParams(current)%eam_atype = eam_ident + EAMList%EAMParams(current)%eam_lattice = lattice_constant + EAMList%EAMParams(current)%eam_nrho = eam_nrho + EAMList%EAMParams(current)%eam_drho = eam_drho + EAMList%EAMParams(current)%eam_nr = eam_nr + EAMList%EAMParams(current)%eam_rcut = rcut + EAMList%EAMParams(current)%eam_Z_r = eam_Z_r + EAMList%EAMParams(current)%eam_rho_r = eam_rho_r + EAMList%EAMParams(current)%eam_F_rho = eam_F_rho - do j = 1, eam_nr(i) - eam_rvals(j,i) = dble(j-1)*eam_dr(i) - enddo + end subroutine newEAMtype - 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 + + do i = 1, EAMList%currentAddition + + EAMList%EAMParams(i)%eam_rvals(1:EAMList%EAMParams(i)%eam_nr) = & + real(EAMList%EAMParams(i)%eam_nr,kind=dp)* & + EAMList%EAMParams(i)%eam_dr + + EAMList%EAMParams(i)%eam_rhovals(1:EAMList%EAMParams(i)%eam_nrho) = & + real(EAMList%EAMParams(i)%eam_nrho,kind=dp)* & + EAMList%EAMParams(i)%eam_drho + ! convert from eV to kcal / mol: - 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 = 0.0E0_DP + do j = 2, EAMList%EAMParams(i)%eam_nr + EAMList%EAMParams(i)%eam_phi_r(j) = (EAMList%EAMParams(i)%eam_Z_r(j)**2)/EAMList%EAMParams(i)%eam_rvals(j) + EAMList%EAMParams(i)%eam_phi_r(j) = EAMList%EAMParams(i)%eam_phi_r(j)*331.999296E0_DP enddo + end do + do i = 1, EAMList%currentAddition + number_r = EAMList%EAMParams(i)%eam_nr + number_rho = EAMList%EAMParams(i)%eam_nrho + + call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, & + EAMList%EAMParams(i)%eam_rho_r, & + EAMList%EAMParams(i)%eam_rho_r_pp, & + 0.0E0_DP, 0.0E0_DP, 'N') + call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, & + EAMList%EAMParams(i)%eam_Z_r, & + EAMList%EAMParams(i)%eam_Z_r_pp, & + 0.0E0_DP, 0.0E0_DP, 'N') + call eam_spline(number_rho, EAMList%EAMParams(i)%eam_rhovals, & + EAMList%EAMParams(i)%eam_F_rho, & + EAMList%EAMParams(i)%eam_F_rho_pp, & + 0.0E0_DP, 0.0E0_DP, 'N') + call eam_spline(number_r, EAMList%EAMParams(i)%eam_rvals, & + EAMList%EAMParams(i)%eam_phi_r, & + EAMList%EAMParams(i)%eam_phi_r_pp, & + 0.0E0_DP, 0.0E0_DP, 'N') + enddo + + current_rcut_max = EAMList%EAMParams(1)%eam_rcut + !! find the smallest rcut + do i = 2, EAMList%currentAddition + current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) end do -#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 + EAM_rcut = current_rcut_max +! do i = 1, EAMList%currentAddition +! EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i +! end do - call mpi_bcast(eam_atype,n_eam_atypes,mpi_integer,0, & - mpi_comm_world,mpi_err) - !! 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) + !! Allocate arrays for force calculation + call allocateEAM(alloc_stat) + if (alloc_stat /= 0 ) then + status = -1 + return + endif - !! 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) + end subroutine init_EAM_FF -#endif - call info('INITIALIZE_EAM', 'creating splines') - - 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 - - do i = 1, n_eam_atypes - eam_atype_map(eam_atype(i)) = i - end do +!! routine checks to see if array is allocated, deallocates array if allocated +!! and then creates the array to the required size + subroutine allocateEAM(status) + integer, intent(out) :: status + integer :: nlocal +#ifdef MPI + integer :: nrow + integer :: ncol +#endif + integer :: alloc_stat - call info('INITIALIZE_EAM','Done creating splines') + nlocal = getNlocal() - return - end subroutine initialize_eam +#ifdef MPI + nrow = getNrow(plan_row) + ncol = getNcol(plan_col) +#endif + if (allocated(frho)) deallocate(frho) + allocate(frho(nlocal),stat=alloc_stat) + if (alloc_stat /= 0) then + status = -1 + return + end if + if (allocated(rho)) deallocate(rho) + allocate(rho(nlocal),stat=alloc_stat) + if (alloc_stat /= 0) then + status = -1 + return + end if + if (allocated(dfrhodrho)) deallocate(dfrhodrho) + allocate(dfrhodrho(nlocal),stat=alloc_stat) + if (alloc_stat /= 0) then + status = -1 + return + end if + +#ifdef MPI - subroutine allocate_eam_atype(n_size_atype) - integer, intent(in) :: n_size_atype + 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 - 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 - subroutine allocate_eam_module(n_size_atype,n_eam_points) - integer, intent(in) :: n_eam_points - integer, intent(in) :: n_size_atype +#endif - 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)) + end subroutine allocateEAM - end subroutine allocate_eam_module - subroutine deallocate_eam_module() + subroutine clean_EAM() - 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) +! clean non-MPI first + frho = 0.0_dp + rho = 0.0_dp + dfrhodrho = 0.0_dp +! clean MPI if needed +#ifdef MPI + frho_row = 0.0_dp + frho_col = 0.0_dp + rho_row = 0.0_dp + rho_col = 0.0_dp + dfrhodrho_row = 0.0_dp + dfrhodrho_col = 0.0_dp +#endif + end subroutine clean_EAM - end subroutine deallocate_eam_module + subroutine allocate_EAMType(eam_n_rho,eam_n_r,thisEAMType,stat) + integer, intent(in) :: eam_n_rho + integer, intent(in) :: eam_n_r + type (EAMType), pointer :: thisEAMType + integer, optional :: stat + integer :: alloc_stat + + + + if (present(stat)) stat = 0 + + allocate(thisEAMType%eam_rvals(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_rhovals(eam_n_rho),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_F_rho(eam_n_rho),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_Z_r(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_rho_r(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_phi_r(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_F_rho_pp(eam_n_rho),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_Z_r_pp(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_rho_r_pp(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + allocate(thisEAMType%eam_phi_r_pp(eam_n_r),stat=alloc_stat) + if (alloc_stat /= 0 ) then + if (present(stat)) stat = -1 + return + end if + + + end subroutine allocate_EAMType + + + subroutine deallocate_EAMType(thisEAMType) + type (EAMtype), pointer :: thisEAMType + + ! free Arrays in reverse order of allocation... + deallocate(thisEAMType%eam_phi_r_pp) + deallocate(thisEAMType%eam_rho_r_pp) + deallocate(thisEAMType%eam_Z_r_pp) + deallocate(thisEAMType%eam_F_rho_pp) + deallocate(thisEAMType%eam_phi_r) + deallocate(thisEAMType%eam_rho_r) + deallocate(thisEAMType%eam_Z_r) + deallocate(thisEAMType%eam_F_rho) + deallocate(thisEAMType%eam_rhovals) + deallocate(thisEAMType%eam_rvals) + + end subroutine deallocate_EAMType + +!! Calculates rho_r + subroutine calc_eam_prepair_rho(atom1,atom2,d,r,rijsq) + integer :: atom1,atom2 + real(kind = dp), dimension(3) :: d + real(kind = dp), intent(inout) :: r + real(kind = dp), intent(inout) :: rijsq + ! value of electron density rho do to atom i at atom j + real(kind = dp) :: rho_i_at_j + ! value of electron density rho do to atom j at atom i + real(kind = dp) :: rho_j_at_i + + ! we don't use the derivatives, dummy variables + real( kind = dp) :: drho,d2rho + integer :: eam_err + + if (cleanme) call clean_EAM + cleanme = .false. + + call calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1) + +#ifdef MPI + rho_col(atom2) = rho_col(atom2) + rho_i_at_j +#else + rho(atom2) = rho(atom2) + rho_i_at_j +#endif + + call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2) + +#ifdef MPI + rho_row(atom1) = rho_row(atom1) + rho_j_at_i +#else + rho(atom1) = rho(atom1) + rho_j_at_i +#endif + end subroutine calc_eam_prepair_rho + + !! Calculate the functional F(rho) for all atoms + subroutine calc_eam_prepair_Frho(nlocal,pot) + integer :: nlocal + real(kind=dp) :: pot + integer :: i,j + real(kind=dp) :: U,U1,U2 + integer :: atype1 + ! reset clean forces to be true at top of calc rho. + cleanme = .true. + +!! Scatter the electron density in pre-pair +#ifdef MPI + call scatter(rho_row,rho,plan_row,eam_err) + if (eam_err /= 0 ) then + write(errMsg,*) " Error scattering rho_row into rho" + call handleError(RoutineName,errMesg) + endif + call scatter(rho_col,rho,plan_col,eam_err) + if (eam_err /= 0 ) then + write(errMsg,*) " Error scattering rho_col into rho" + call handleError(RoutineName,errMesg) + endif +#endif + + do i = 1, nlocal + call calc_eam_frho(rho(i),u,u1,u2,atype1) + frho(i) = u + dfrhodrho(i) = u1 +! d2frhodrhodrho(i) = u2 + pot = pot + u + enddo + +#ifdef MPI + !! communicate f(rho) and derivatives + call gather(frho,frho_row,plan_row, eam_err) + if (eam_err /= 0) then + call handleError("cal_eam_forces()","MPI gather frho_row failure") + endif + call gather(dfrhodrho,dfrhodrho_row,plan_row, eam_err) + if (eam_err /= 0) then + call handleError("cal_eam_forces()","MPI gather dfrhodrho_row failure") + endif + call gather(frho,frho_col,plan_col, eam_err) + if (eam_err /= 0) then + call handleError("cal_eam_forces()","MPI gather frho_col failure") + endif + call gather(dfrhodrho,dfrhodrho_col,plan_col, eam_err) + if (eam_err /= 0) then + call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure") + endif + + if (nmflag) then + call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row) + call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col) + endif +#endif + + + end subroutine calc_eam_prepair_Frho + + + + subroutine calc_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress) !Arguments integer, intent(in) :: atom1, atom2 @@ -261,34 +513,56 @@ contains real( kind = dp ), intent(in), dimension(3) :: d logical, intent(in) :: do_pot, do_stress + real( kind = dp ) :: drdx,drdy,drdz + real( kind = dp ) :: phab,pha,dvpdr,d2vpdrdr + real( kind = dp ) :: rha,drha,d2rha, dpha + real( kind = dp ) :: rhb,drhb,d2rhb, dphb + real( kind = dp ) :: dudr + real( kind = dp ) :: rci,rcj + real( kind = dp ) :: drhoidr,drhojdr + real( kind = dp ) :: d2rhoidrdr + real( kind = dp ) :: d2rhojdrdr + real( kind = dp ) :: Fx,Fy,Fz + real( kind = dp ) :: r,d2pha,phb,d2phb + integer :: id1,id2 + integer :: atype1,atype2 + + !Local Variables + + phab = 0.0E0_DP + dvpdr = 0.0E0_DP + d2vpdrdr = 0.0E0_DP + + if (rij .lt. EAM_rcut) then +#ifdef IS_MPI +!!!!! FIX ME + atype1 = atid_row(atom1) +#else + atype1 = atid(atom1) +#endif - r = dsqrt(rijsq) - efr(1,j) = -rxij - efr(2,j) = -ryij - efr(3,j) = -rzij - - + 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(atype1)) -#ifdef MPI - atype2 = ident_col(j) + ! rci = eam_rcut(eam_atype_map(atom1)) +#ifdef IS_MPI + atype2 = atid_col(atom2) #else - atype2 = ident(j) + atype2 = atid(atom2) #endif - + call calc_eam_rho(r, rhb, drhb, d2rhb, atype2) call calc_eam_phi(r, phb, dphb, d2phb, atype2) - rcj = eam_rcut(eam_atype_map(atype2)) + ! rcj = eam_rcut(eam_atype_map(atype2)) - phab = 0.0E0_DP - dvpdr = 0.0E0_DP - d2vpdrdr = 0.0E0_DP - if (r.lt.rci) then phab = phab + 0.5E0_DP*(rhb/rha)*pha dvpdr = dvpdr + 0.5E0_DP*((rhb/rha)*dpha + & @@ -298,8 +572,8 @@ 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 phab = phab + 0.5E0_DP*(rha/rhb)*phb dvpdr = dvpdr + 0.5E0_DP*((rha/rhb)*dphb + & @@ -309,294 +583,269 @@ 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) & + 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) #endif + fx = dudr * drdx + fy = dudr * drdy + fz = dudr * drdz - do dim = 1, 3 - drdx1 = efr(dim,j) / r - ftmp = dudr * drdx1 - #ifdef MPI - f_col(dim,j) = f_col(dim,j) - ftmp - f_row(dim,i) = f_row(dim,i) + ftmp -#else - f(dim,j) = f(dim,j) - ftmp - f(dim,i) = f(dim,i) + ftmp -#endif + if (do_pot) then + pot_Row(atom1) = pot_Row(atom1) + phab*0.5 + pot_Col(atom2) = pot_Col(atom2) + phab*0.5 + end if - if (nmflag) then - idim = 3 * (i-1) + dim - jdim = 3 * (j-1) + dim + 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 - do dim2 = 1, 3 + 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 - kt1 = d2 * efr(dim,j) * efr(dim2,j)/r/r - kt2 = - dudr * efr(dim,j) * efr(dim2,j)/r/r/r + f(1,atom1) = f(1,atom1) + fx + f(2,atom1) = f(2,atom1) + fy + f(3,atom1) = f(3,atom1) + fz - if (dim.eq.dim2) then - kt3 = dudr / r - else - kt3 = 0.0E0_DP - endif + f(1,atom2) = f(1,atom2) - fx + f(2,atom2) = f(2,atom2) - fy + f(3,atom2) = f(3,atom2) - fz +#endif - ! The factor of 2 below is to compensate for - ! overcounting. - ! Mass weighting is done separately... - 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 + if (do_stress) then - d(jdim, idim2) = d(jdim,idim2) - ktmp - d(jdim2, idim) = d(jdim2,idim) - ktmp +#ifdef MPI + id1 = tagRow(atom1) + id2 = tagColumn(atom2) +#else + id1 = atom1 + id2 = atom2 +#endif - d(jdim, jdim2) = d(jdim,jdim2) + ktmp - d(jdim2, jdim) = d(jdim2,jdim) + ktmp + if (molMembershipList(id1) .ne. molMembershipList(id2)) then + + tau_Temp(1) = tau_Temp(1) + fx * d(1) + tau_Temp(2) = tau_Temp(2) + fx * d(2) + tau_Temp(3) = tau_Temp(3) + fx * d(3) + tau_Temp(4) = tau_Temp(4) + fy * d(1) + tau_Temp(5) = tau_Temp(5) + fy * d(2) + tau_Temp(6) = tau_Temp(6) + fy * d(3) + tau_Temp(7) = tau_Temp(7) + fz * d(1) + tau_Temp(8) = tau_Temp(8) + fz * d(2) + tau_Temp(9) = tau_Temp(9) + fz * d(3) + virial_Temp = virial_Temp + & + (tau_Temp(1) + tau_Temp(5) + tau_Temp(9)) - enddo endif - enddo - + endif + endif - enddo -endif -enddo +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 calc_eam_pair + subroutine eam_splint(nx, xa, ya, yppa, x, y, dy, d2y) -subroutine calc_eam_rho(r, rho, drho, d2rho, atype) - ! include 'headers/sizes.h' + integer :: atype, nx, j + real( kind = DP ), dimension(:) :: xa + real( kind = DP ), dimension(:) :: ya + real( kind = DP ), dimension(:) :: yppa + real( kind = DP ) :: x, y, dy, d2y + real( kind = DP ) :: del, h, a, b, c, d -integer atype, etype, number_r -real( kind = DP ) :: r, rho, drho, d2rho -integer :: i -etype = eam_atype_map(atype) + + ! 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 -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 + ! check to make sure we're inside the spline range: + if ((j.gt.nx).or.(j.lt.1)) then + write(default_error,*) "EAM_splint: x is outside bounds of spline" + endif + ! check to make sure we haven't screwed up the calculation of j: + if ((x.lt.xa(j)).or.(x.gt.xa(j+1))) then + if (j.ne.nx) then + write(default_error,*) "EAM_splint: x is outside bounding range" + endif + endif -return -end subroutine calc_eam_rho + 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) -subroutine calc_eam_frho(dens, u, u1, u2, atype) + end subroutine eam_splint - ! include 'headers/sizes.h' + subroutine eam_spline(nx, xa, ya, yppa, yp1, ypn, boundary) -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 + ! yp1 and ypn are the first derivatives of y at the two endpoints + ! if boundary is 'L' the lower derivative is used + ! if boundary is 'U' the upper derivative is used + ! if boundary is 'B' then both derivatives are used + ! if boundary is anything else, then both derivatives are assumed to be 0 + + integer :: nx, i, k, max_array_size + + real( kind = DP ), dimension(:) :: xa + real( kind = DP ), dimension(:) :: ya + real( kind = DP ), dimension(:) :: yppa + real( kind = DP ), dimension(size(xa)) :: u + real( kind = DP ) :: yp1,ypn,un,qn,sig,p + character boundary + + + if ((boundary.eq.'l').or.(boundary.eq.'L').or. & + (boundary.eq.'b').or.(boundary.eq.'B')) then + yppa(1) = -0.5E0_DP + u(1) = (3.0E0_DP/(xa(2)-xa(1)))*((ya(2)-& + ya(1))/(xa(2)-xa(1))-yp1) + else + yppa(1) = 0.0E0_DP + u(1) = 0.0E0_DP + endif + + do i = 2, nx - 1 + sig = (xa(i) - xa(i-1)) / (xa(i+1) - xa(i-1)) + p = sig * yppa(i-1) + 2.0E0_DP + yppa(i) = (sig - 1.0E0_DP) / p + u(i) = (6.0E0_DP*((ya(i+1)-ya(i))/(xa(i+1)-xa(i)) - & + (ya(i)-ya(i-1))/(xa(i)-xa(i-1)))/ & + (xa(i+1)-xa(i-1)) - sig * u(i-1))/p + enddo + + if ((boundary.eq.'u').or.(boundary.eq.'U').or. & + (boundary.eq.'b').or.(boundary.eq.'B')) then + qn = 0.5E0_DP + un = (3.0E0_DP/(xa(nx)-xa(nx-1)))* & + (ypn-(ya(nx)-ya(nx-1))/(xa(nx)-xa(nx-1))) + else + qn = 0.0E0_DP + un = 0.0E0_DP + endif -subroutine calc_eam_phi(r, phi, dphi, d2phi, atype) + yppa(nx)=(un-qn*u(nx-1))/(qn*yppa(nx-1)+1.0E0_DP) + + do k = nx-1, 1, -1 + yppa(k)=yppa(k)*yppa(k+1)+u(k) + enddo + end subroutine eam_spline -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