--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/07/17 19:25:51 631 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2003/12/17 20:13:33 882 @@ -4,7 +4,8 @@ module eam use force_globals use status use atype_module -#ifdef MPI + use Vector_class +#ifdef IS_MPI use mpiSimulation #endif implicit none @@ -14,12 +15,17 @@ module eam 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 + 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 +55,20 @@ 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 - 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 - +#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 @@ -71,18 +79,19 @@ module eam end type EAMTypeList - type (eamTypeList) :: EAMList + type (eamTypeList), save :: EAMList !! standard eam stuff public :: init_EAM_FF -! public :: EAM_new_rcut -! public :: do_EAM_pair + public :: setCutoffEAM + public :: do_eam_pair public :: newEAMtype + public :: calc_eam_prepair_rho + public :: calc_eam_preforce_Frho + public :: clean_EAM - - contains @@ -106,8 +115,10 @@ contains integer :: alloc_stat integer :: current integer,pointer :: Matchlist(:) => null() + status = 0 + !! 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. @@ -123,13 +134,12 @@ contains current = EAMList%currentAddition - !call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat) + call allocate_EAMType(eam_nrho,eam_nr,EAMList%EAMParams(current),stat=alloc_stat) if (alloc_stat /= 0) then status = -1 return end if - ! this is a possible bug, we assume a correspondence between the vector atypes and ! EAMAtypes @@ -138,6 +148,7 @@ contains 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 @@ -154,26 +165,41 @@ contains 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 - - 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 +! Build array of r values + do j = 1,EAMList%EAMParams(i)%eam_nr + EAMList%EAMParams(i)%eam_rvals(j) = & + real(j-1,kind=dp)* & + EAMList%EAMParams(i)%eam_dr + end do +! Build array of rho values + do j = 1,EAMList%EAMParams(i)%eam_nrho + EAMList%EAMParams(i)%eam_rhovals(j) = & + real(j-1,kind=dp)* & + EAMList%EAMParams(i)%eam_drho + end do ! convert from eV to kcal / mol: EAMList%EAMParams(i)%eam_F_rho = EAMList%EAMParams(i)%eam_F_rho * 23.06054E0_DP ! precompute the pair potential and get it into kcal / mol: - EAMList%EAMParams(i)%eam_phi_r = 0.0E0_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 + do i = 1, EAMList%currentAddition number_r = EAMList%EAMParams(i)%eam_nr @@ -196,26 +222,27 @@ contains 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 - EAM_rcut = current_rcut_max +! 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 + +! EAM_rcut = current_rcut_max +! EAM_rcut_orig = current_rcut_max ! do i = 1, EAMList%currentAddition ! EAMList%EAMParam(i)s%eam_atype_map(eam_atype(i)) = i ! end do - - !! Allocate arrays for force calculation - call allocateEAM(alloc_stat) - if (alloc_stat /= 0 ) then - status = -1 - return - endif - + + call allocateEAM(alloc_stat) + if (alloc_stat /= 0 ) then + write(*,*) "allocateEAM failed" + status = -1 + return + endif + end subroutine init_EAM_FF !! routine checks to see if array is allocated, deallocates array if allocated @@ -224,7 +251,7 @@ contains integer, intent(out) :: status integer :: nlocal -#ifdef MPI +#ifdef IS_MPI integer :: nrow integer :: ncol #endif @@ -232,8 +259,8 @@ contains nlocal = getNlocal() - -#ifdef MPI + status = 0 +#ifdef IS_MPI nrow = getNrow(plan_row) ncol = getNcol(plan_col) #endif @@ -250,15 +277,31 @@ 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(rho_tmp)) deallocate(rho_tmp) + allocate(rho_tmp(nlocal),stat=alloc_stat) + if (alloc_stat /= 0) then + status = -1 + return + end if + + if (allocated(frho_row)) deallocate(frho_row) allocate(frho_row(nrow),stat=alloc_stat) if (alloc_stat /= 0) then @@ -277,7 +320,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,24 +348,44 @@ 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 + +!! 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 + EAM_rcut = rcut - subroutine clean_EAM() + end subroutine setCutoffEAM -! clean non-MPI first + + + 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 MPI +#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 @@ -326,7 +396,7 @@ contains 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 + type (EAMType) :: thisEAMType integer, optional :: stat integer :: alloc_stat @@ -421,60 +491,134 @@ contains real( kind = dp) :: drho,d2rho integer :: eam_err - if (cleanme) call clean_EAM - cleanme = .false. + integer :: myid_atom1 + integer :: myid_atom2 - call calc_eam_rho(r,rho_i_at_j,drho,d2rho,atom1) +! check to see if we need to be cleaned at the start of a force loop + -#ifdef MPI - rho_col(atom2) = rho_col(atom2) + rho_i_at_j + + +#ifdef IS_MPI + myid_atom1 = atid_Row(atom1) + myid_atom2 = atid_Col(atom2) #else - rho(atom2) = rho(atom2) + rho_i_at_j + myid_atom1 = atid(atom1) + myid_atom2 = atid(atom2) #endif - call calc_eam_rho(r,rho_j_at_i,drho,d2rho,atom2) + if (r.lt.EAMList%EAMParams(myid_atom1)%eam_rcut) then -#ifdef MPI - rho_row(atom1) = rho_row(atom1) + rho_j_at_i + + + 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(atom1) = rho(atom1) + rho_j_at_i + 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 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 - ! reset clean forces to be true at top of calc rho. - cleanme = .true. + integer :: me + integer :: n_rho_points -!! Scatter the electron density in pre-pair -#ifdef MPI + + cleanme = .true. + !! Scatter the electron density from pre-pair calculation back to local atoms +#ifdef IS_MPI call scatter(rho_row,rho,plan_row,eam_err) if (eam_err /= 0 ) then write(errMsg,*) " Error scattering rho_row into rho" call handleError(RoutineName,errMesg) endif - call scatter(rho_col,rho,plan_col,eam_err) + 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 + + rho(1:nlocal) = rho(1:nlocal) + rho_tmp(1:nlocal) #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(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") @@ -492,28 +636,34 @@ 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,46 +674,75 @@ 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 + ! 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 - 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 +753,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,14 +770,14 @@ contains d2rhojdrdr = d2rhb -#ifdef MPI +#ifdef IS_MPI dudr = drhojdr*dfrhodrho_row(atom1)+drhoidr*dfrhodrho_col(atom2) & + dvpdr - + #else dudr = drhojdr*dfrhodrho(atom1)+drhoidr*dfrhodrho(atom2) & + dvpdr - + ! write(*,*) "Atom1,Atom2, dfrhodrho(atom1) dfrhodrho(atom2): ", atom1,atom2,dfrhodrho(atom1),dfrhodrho(atom2) #endif fx = dudr * drdx @@ -606,7 +785,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 +794,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 + if(do_pot) then + pot = pot + phab + end if + 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 do_eam_pair -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 - - 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 +914,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 +942,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