--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/07/24 19:57:35 650 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2003/09/05 21:28:52 747 @@ -15,7 +15,7 @@ module eam logical, save :: EAM_FF_initialized = .false. integer, save :: EAM_Mixing_Policy real(kind = dp), save :: EAM_rcut - real(kind = dp), save :: EAM_rcut_orig + logical, save :: haveRcut = .false. character(len = statusMsgSize) :: errMesg integer :: eam_err @@ -60,14 +60,14 @@ module eam !! Arrays for MPI storage #ifdef IS_MPI - real( kind = dp), dimension(:), allocatable :: dfrhodrho_col - real( kind = dp), dimension(:), allocatable :: dfrhodrho_row - real( kind = dp), dimension(:), allocatable :: frho_row - real( kind = dp), dimension(:), allocatable :: frho_col - real( kind = dp), dimension(:), allocatable :: rho_row - real( kind = dp), dimension(:), allocatable :: rho_col - real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_col - real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho_row + 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 :: d2frhodrhodrho_col + real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho_row #endif type, private :: EAMTypeList @@ -78,18 +78,18 @@ module eam end type EAMTypeList - type (eamTypeList) :: EAMList + type (eamTypeList), save :: EAMList !! standard eam stuff public :: init_EAM_FF -! public :: EAM_new_rcut + public :: setCutoffEAM public :: do_eam_pair public :: newEAMtype public :: calc_eam_prepair_rho public :: calc_eam_preforce_Frho - + public :: clean_EAM contains @@ -114,9 +114,10 @@ contains integer :: alloc_stat integer :: current integer,pointer :: Matchlist(:) => null() - type (EAMtype), pointer :: makeEamtype => null() + status = 0 + !! Assume that atypes has already been set and get the total number of types in atypes !! Also assume that every member of atypes is a EAM model. @@ -132,12 +133,11 @@ contains current = EAMList%currentAddition - call allocate_EAMType(eam_nrho,eam_nr,makeEamtype,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 - makeEamtype => EAMList%EAMParams(current) ! this is a possible bug, we assume a correspondence between the vector atypes and ! EAMAtypes @@ -147,6 +147,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 @@ -163,28 +164,40 @@ contains integer :: alloc_stat integer :: number_r, number_rho + 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 @@ -207,20 +220,19 @@ contains EAMList%EAMParams(i)%eam_phi_r_pp, & 0.0E0_DP, 0.0E0_DP, 'N') enddo - - current_rcut_max = EAMList%EAMParams(1)%eam_rcut + + +! current_rcut_max = EAMList%EAMParams(1)%eam_rcut !! find the smallest rcut for any eam atype - do i = 2, EAMList%currentAddition - current_rcut_max = min(current_rcut_max,EAMList%EAMParams(i)%eam_rcut) - end do +! 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 +! 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 @@ -336,10 +348,23 @@ contains 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 - subroutine clean_EAM() + EAM_rcut = rcut -! clean non-IS_MPI first + end subroutine setCutoffEAM + + + + subroutine clean_EAM() + + ! clean non-IS_MPI first frho = 0.0_dp rho = 0.0_dp dfrhodrho = 0.0_dp @@ -359,7 +384,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 @@ -458,10 +483,10 @@ contains integer :: myid_atom2 ! check to see if we need to be cleaned at the start of a force loop - if (cleanme) call clean_EAM - cleanme = .false. - + + + #ifdef IS_MPI myid_atom1 = atid_Row(atom1) myid_atom2 = atid_Col(atom2) @@ -471,9 +496,9 @@ contains #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, & @@ -487,6 +512,7 @@ contains #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 @@ -506,6 +532,7 @@ contains #endif endif + end subroutine calc_eam_prepair_rho @@ -521,10 +548,10 @@ contains integer :: atype1 integer :: me integer :: n_rho_points - ! reset clean forces to be true at top of calc rho. - cleanme = .true. -!! Scatter the electron density from pre-pair calculation back to local atoms + + 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 @@ -562,14 +589,15 @@ contains end if - frho(i) = u - dfrhodrho(i) = u1 - d2frhodrhodrho(i) = u2 + 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) @@ -599,6 +627,7 @@ contains endif #endif + end subroutine calc_eam_preforce_Frho @@ -631,16 +660,14 @@ contains 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 @@ -727,11 +754,11 @@ contains #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 @@ -753,8 +780,11 @@ contains 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 @@ -802,10 +832,7 @@ contains #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