--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/08/07 00:47:33 669 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2004/01/05 22:49:14 898 @@ -51,11 +51,11 @@ module eam !! Arrays for derivatives used in force calculation - real( kind = dp),save, dimension(:), allocatable :: frho - real( kind = dp),save, dimension(:), allocatable :: rho + real( kind = dp), dimension(:), allocatable :: frho + real( kind = dp), dimension(:), allocatable :: rho - real( kind = dp),save, dimension(:), allocatable :: dfrhodrho - real( kind = dp),save, dimension(:), allocatable :: d2frhodrhodrho + real( kind = dp), dimension(:), allocatable :: dfrhodrho + real( kind = dp), dimension(:), allocatable :: d2frhodrhodrho !! Arrays for MPI storage @@ -66,6 +66,7 @@ module eam 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 @@ -78,7 +79,7 @@ module eam end type EAMTypeList - type (eamTypeList) :: EAMList + type (eamTypeList), save :: EAMList !! standard eam stuff @@ -89,7 +90,7 @@ module eam public :: newEAMtype public :: calc_eam_prepair_rho public :: calc_eam_preforce_Frho - + public :: clean_EAM contains @@ -164,14 +165,15 @@ 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 ! Build array of r values @@ -221,7 +223,6 @@ contains 0.0E0_DP, 0.0E0_DP, 'N') enddo - ! current_rcut_max = EAMList%EAMParams(1)%eam_rcut !! find the smallest rcut for any eam atype ! do i = 2, EAMList%currentAddition @@ -233,15 +234,15 @@ contains ! 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 @@ -249,7 +250,6 @@ contains subroutine allocateEAM(status) integer, intent(out) :: status - integer :: nlocal #ifdef IS_MPI integer :: nrow integer :: ncol @@ -257,8 +257,7 @@ contains integer :: alloc_stat - nlocal = getNlocal() - + status = 0 #ifdef IS_MPI nrow = getNrow(plan_row) ncol = getNcol(plan_col) @@ -293,6 +292,14 @@ contains #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 @@ -365,6 +372,7 @@ contains subroutine clean_EAM() + ! clean non-IS_MPI first frho = 0.0_dp rho = 0.0_dp @@ -375,6 +383,7 @@ contains 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 @@ -484,15 +493,10 @@ contains integer :: myid_atom2 ! check to see if we need to be cleaned at the start of a force loop - - if (cleanme) then - call clean_EAM - cleanme = .false. - end if - + #ifdef IS_MPI myid_atom1 = atid_Row(atom1) myid_atom2 = atid_Col(atom2) @@ -502,9 +506,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, & @@ -518,6 +522,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 @@ -527,7 +532,7 @@ contains EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, & r, rho_j_at_i,drho,d2rho) - + #ifdef IS_MPI @@ -537,8 +542,11 @@ contains #endif endif - + + + + end subroutine calc_eam_prepair_rho @@ -554,24 +562,27 @@ 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 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 + !! Calculate F(rho) and derivative do atom = 1, nlocal me = atid(atom) @@ -595,13 +606,14 @@ 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 @@ -632,6 +644,7 @@ contains endif #endif + end subroutine calc_eam_preforce_Frho @@ -643,7 +656,7 @@ contains 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 ), dimension(3,nLocal) :: f real( kind = dp ), intent(in), dimension(3) :: d logical, intent(in) :: do_pot, do_stress @@ -663,12 +676,11 @@ contains 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 @@ -837,10 +849,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