--- trunk/OOPSE/libmdtools/calc_eam.F90 2003/08/27 16:25:11 730 +++ trunk/OOPSE/libmdtools/calc_eam.F90 2004/05/27 00:48:12 1198 @@ -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 @@ -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 @@ -234,12 +235,14 @@ contains ! 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 @@ -247,19 +250,17 @@ contains subroutine allocateEAM(status) integer, intent(out) :: status - integer :: nlocal #ifdef IS_MPI - integer :: nrow - integer :: ncol + integer :: nAtomsInRow + integer :: nAtomsInCol #endif integer :: alloc_stat - nlocal = getNlocal() - + status = 0 #ifdef IS_MPI - nrow = getNrow(plan_row) - ncol = getNcol(plan_col) + nAtomsInRow = getNatomsInRow(plan_atom_row) + nAtomsInCol = getNatomsInCol(plan_atom_col) #endif if (allocated(frho)) deallocate(frho) @@ -291,26 +292,34 @@ 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) + allocate(frho_row(nAtomsInRow),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) + allocate(rho_row(nAtomsInRow),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) + allocate(dfrhodrho_row(nAtomsInRow),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) + allocate(d2frhodrhodrho_row(nAtomsInRow),stat=alloc_stat) if (alloc_stat /= 0) then status = -1 return @@ -320,25 +329,25 @@ contains ! Now do column arrays if (allocated(frho_col)) deallocate(frho_col) - allocate(frho_col(ncol),stat=alloc_stat) + allocate(frho_col(nAtomsInCol),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) + allocate(rho_col(nAtomsInCol),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) + allocate(dfrhodrho_col(nAtomsInCol),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) + allocate(d2frhodrhodrho_col(nAtomsInCol),stat=alloc_stat) if (alloc_stat /= 0) then status = -1 return @@ -374,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 @@ -522,7 +532,7 @@ contains EAMList%EAMParams(myid_atom2)%eam_rho_r_pp, & r, rho_j_at_i,drho,d2rho) - + #ifdef IS_MPI @@ -531,6 +541,10 @@ contains rho(atom1) = rho(atom1) + rho_j_at_i #endif endif + + + + end subroutine calc_eam_prepair_rho @@ -553,19 +567,22 @@ contains 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) + call scatter(rho_row,rho,plan_atom_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_atom_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) @@ -600,19 +617,19 @@ contains #ifdef IS_MPI !! communicate f(rho) and derivatives back into row and column arrays - call gather(frho,frho_row,plan_row, eam_err) + call gather(frho,frho_row,plan_atom_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) + call gather(dfrhodrho,dfrhodrho_row,plan_atom_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) + call gather(frho,frho_col,plan_atom_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) + call gather(dfrhodrho,dfrhodrho_col,plan_atom_col, eam_err) if (eam_err /= 0) then call handleError("cal_eam_forces()","MPI gather dfrhodrho_col failure") endif @@ -622,8 +639,8 @@ contains if (nmflag) then - call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_row) - call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_col) + call gather(d2frhodrhodrho,d2frhodrhodrho_row,plan_atom_row) + call gather(d2frhodrhodrho,d2frhodrhodrho_col,plan_atom_col) endif #endif @@ -634,14 +651,17 @@ contains !! Does EAM pairwise Force calculation. - subroutine do_eam_pair(atom1,atom2,d,rij,r2,pot,f,do_pot,do_stress) + subroutine do_eam_pair(atom1, atom2, d, rij, r2, sw, vpair, fpair, & + pot, f, do_pot) !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 ), intent(inout), dimension(3) :: fpair + + logical, intent(in) :: do_pot real( kind = dp ) :: drdx,drdy,drdz real( kind = dp ) :: d2 @@ -794,6 +814,23 @@ contains f(3,atom2) = f(3,atom2) - fz #endif + vpair = vpair + phab +#ifdef IS_MPI + id1 = AtomRowToGlobal(atom1) + id2 = AtomColToGlobal(atom2) +#else + id1 = atom1 + id2 = atom2 +#endif + + if (molMembershipList(id1) .ne. molMembershipList(id2)) then + + fpair(1) = fpair(1) + fx + fpair(2) = fpair(2) + fy + fpair(3) = fpair(3) + fz + + endif + if (nmflag) then drhoidr = drha @@ -817,40 +854,8 @@ contains drhojdr*drhojdr*d2frhodrhodrho(atom1) #endif end if - - - - 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 - - + endif end subroutine do_eam_pair