From 1c595144cdc5cb09d2f2c3c51932b016d3363262 Mon Sep 17 00:00:00 2001 From: Alberto Otero de la Roza Date: Tue, 13 Feb 2024 19:22:04 +0100 Subject: [PATCH] implemented ndiv in list_near_lattice_points --- src/crystalmod.f90 | 4 +++- src/crystalmod@env.f90 | 50 ++++++++++++++++++++++++++++-------------- src/global@proc.F90 | 47 +++++++++++++++++++++++---------------- 3 files changed, 65 insertions(+), 36 deletions(-) diff --git a/src/crystalmod.f90 b/src/crystalmod.f90 index f615adfa..47438a96 100644 --- a/src/crystalmod.f90 +++ b/src/crystalmod.f90 @@ -460,7 +460,8 @@ end subroutine promolecular_env module subroutine find_asterisms_covalent(c) class(crystal), intent(inout) :: c end subroutine find_asterisms_covalent - module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,up2n,nozero) + module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,ndiv,& + up2d,up2n,nozero) class(crystal), intent(inout) :: c real*8, intent(in) :: xp(3) integer, intent(in) :: icrd @@ -468,6 +469,7 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u integer, intent(out) :: nat real*8, allocatable, intent(inout), optional :: dist(:) integer, allocatable, intent(inout), optional :: lvec(:,:) + integer, intent(in), optional :: ndiv(3) real*8, intent(in), optional :: up2d integer, intent(in), optional :: up2n logical, intent(in), optional :: nozero diff --git a/src/crystalmod@env.f90 b/src/crystalmod@env.f90 index 35b7d4ce..ac01ae61 100644 --- a/src/crystalmod@env.f90 +++ b/src/crystalmod@env.f90 @@ -810,6 +810,7 @@ end subroutine find_asterisms_covalent !> regardless of the value of sorted. !> !> Optional input: + !> - ndiv = divide the parent lattice vectors by ndiv; useful for grids !> - nozero = disregard zero-distance lattice points. !> !> Output: @@ -820,9 +821,10 @@ end subroutine find_asterisms_covalent !> - lvec(3,nid) = lattice vector coordinates (crystallographic). !> !> This routine is thread-safe. - module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,up2n,nozero) + module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,ndiv,& + up2d,up2n,nozero) use types, only: realloc - use tools, only: mergesort + use tools, only: mergesort, wscell use tools_math, only: det3, cross use tools_io, only: ferror, faterr use param, only: icrd_cart, icrd_crys @@ -833,6 +835,7 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u integer, intent(out) :: nat real*8, allocatable, intent(inout), optional :: dist(:) integer, allocatable, intent(inout), optional :: lvec(:,:) + integer, intent(in), optional :: ndiv(3) real*8, intent(in), optional :: up2d integer, intent(in), optional :: up2n logical, intent(in), optional :: nozero @@ -842,10 +845,11 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u integer :: nb, nshellb, nsafe, idx integer :: nx(3), nn, i0(3), i1(3), i, j, k, lvecx(3) real*8 :: up2d_2, x(3), lvecini(3), blockcv(3), blockomega, xorigc(3) - real*8 :: xdif(3), dd, dmax + real*8 :: xdif(3), dd, dmax, x0(3) integer, allocatable :: idb(:,:), iaux(:) integer, allocatable :: at_lvec(:,:), iord(:) real*8, allocatable :: at_dist(:) + real*8 :: x2c(3,3), x2xr(3,3), xr2x(3,3), xr2c(3,3), c2xr(3,3) real*8, parameter :: eps = 1d-20 @@ -861,17 +865,31 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u call ferror("list_near_lattice_points","must give one of up2d or up2n",faterr) end if - ! translate the point to the main cell + if (present(ndiv)) then + do i = 1, 3 + x2c(:,i) = c%m_x2c(:,i) / real(ndiv(i),8) + end do + x0 = xp * ndiv + call wscell(x2c,.true.,m_xr2x=xr2x,m_xr2c=xr2c,m_x2xr=x2xr,m_c2xr=c2xr) + else + x0 = xp + xr2c = c%m_xr2c + c2xr = c%m_c2xr + xr2x = c%m_xr2x + x2xr = c%m_x2xr + end if + + ! translate the point to the main reduced cell if (icrd == icrd_crys) then - x = c%x2xr(xp) + x = matmul(x2xr,x0) elseif (icrd == icrd_cart) then - x = c%c2xr(xp) + x = matmul(c2xr,x0) else - x = xp + x = x0 end if lvecini = floor(x) x = x - lvecini - xorigc = c%xr2c(x) + xorigc = matmul(xr2c,x) ! allocate space for lattice points nat = 0 @@ -892,7 +910,7 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u call make_block_shell(c,nshellb,nb,idb,dmax) do i = 1, nb - xdif = matmul(c%m_xr2c,real(idb(:,i),8)) - xorigc + xdif = matmul(xr2c,real(idb(:,i),8)) - xorigc dd = dot_product(xdif,xdif) ! check if we should add the atom to the list @@ -905,7 +923,7 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u call realloc(at_lvec,3,2*nat) end if at_dist(nat) = sqrt(dd) - at_lvec(:,nat) = nint(matmul(c%m_xr2x,lvecini + idb(:,i))) + at_lvec(:,nat) = nint(matmul(xr2x,lvecini + idb(:,i))) end if end do @@ -948,10 +966,10 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u ! calculate the number of blocks in each direction required for satifying ! that the largest sphere in the super-block has radius > dmax. ! r = Vblock / 2 / max(cv(3)/n3,cv(2)/n2,cv(1)/n1) - blockcv(1) = norm2(cross(c%m_xr2c(:,2),c%m_xr2c(:,3))) - blockcv(2) = norm2(cross(c%m_xr2c(:,1),c%m_xr2c(:,3))) - blockcv(3) = norm2(cross(c%m_xr2c(:,1),c%m_xr2c(:,2))) - blockomega = det3(c%m_xr2c) + blockcv(1) = norm2(cross(xr2c(:,2),xr2c(:,3))) + blockcv(2) = norm2(cross(xr2c(:,1),xr2c(:,3))) + blockcv(3) = norm2(cross(xr2c(:,1),xr2c(:,2))) + blockomega = det3(xr2c) nx = ceiling(blockcv / (0.5d0 * blockomega / max(up2d,1d-40))) ! define the search space @@ -970,7 +988,7 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u do j = i0(2), i1(2) do k = i0(3), i1(3) lvecx = (/i,j,k/) - xdif = matmul(c%m_xr2c,real(lvecx,8)) - xorigc + xdif = matmul(xr2c,real(lvecx,8)) - xorigc dd = dot_product(xdif,xdif) ! check if we should add the lattice point to the list @@ -987,7 +1005,7 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u call realloc(at_lvec,3,2*nat) end if at_dist(nat) = sqrt(dd) - at_lvec(:,nat) = nint(matmul(c%m_xr2x,lvecini + lvecx)) + at_lvec(:,nat) = nint(matmul(xr2x,lvecini + lvecx)) end if end do ! k = i0(3), i1(3) end do ! j = i0(2), i1(2) diff --git a/src/global@proc.F90 b/src/global@proc.F90 index 0ee083f9..71bfa60c 100644 --- a/src/global@proc.F90 +++ b/src/global@proc.F90 @@ -71,10 +71,10 @@ module subroutine critic_main() logical :: doref, doname, doflags #endif ! xxxx - real*8 :: xp(3), dmax + real*8 :: xp(3), dmax, xred(3,3) real*8, allocatable :: dist(:) integer, allocatable :: lvec(:,:), eid(:) - integer :: nat, ierr, lvec2(3) + integer :: nat, ierr, lvec2(3), n(3) type(environ) :: e ! Start reading @@ -634,34 +634,43 @@ module subroutine critic_main() elseif (equal(word,'temp')) then ! xxxx - xp = (/-15.3d0,0.1d0,21.4d0/) - ! xp = 0d0 + !! build lattice test + xp = (/10.3d0,0.1d0,-3.4d0/) dmax = 30d0 - ! call e%build_lattice(sy%c%m_x2c,150d0) - ! call e%list_near_atoms(xp,icrd_crys,.true.,nat,ierr,eid=eid,dist=dist,lvec=lvec2,up2n=2,nozero=.true.) + call e%build_lattice(sy%c%m_x2c,150d0) + call e%list_near_atoms(xp,icrd_crys,.true.,nat,ierr,eid=eid,dist=dist,lvec=lvec2,up2n=10,nozero=.true.) + do i = 1, nat + write (*,*) i, e%xr2x(e%at(eid(i))%x) + lvec2, dist(i) + end do + write (*,*) + + call sy%c%list_near_lattice_points(xp,icrd_crys,.true.,nat,dist,lvec,up2n=10,nozero=.true.) + do i = 1, nat + write (*,*) i, lvec(:,i), dist(i) + end do + + ! !! ndiv + ! n = (/100,80,120/) + ! xp = (/0.3d0,0.1d0,0.4d0/) * n + ! dmax = 30d0 + ! xred(:,1) = sy%c%m_x2c(:,1) / n(1) + ! xred(:,2) = sy%c%m_x2c(:,2) / n(2) + ! xred(:,3) = sy%c%m_x2c(:,3) / n(3) + ! call e%build_lattice(xred,5d0) + ! call e%list_near_atoms(xp,icrd_crys,.true.,nat,ierr,eid=eid,dist=dist,lvec=lvec2,up2n=10,nozero=.true.) ! do i = 1, nat ! write (*,*) i, e%xr2x(e%at(eid(i))%x) + lvec2, dist(i) ! end do ! write (*,*) - ! call sy%c%list_near_lattice_points(xp,icrd_crys,.true.,nat,dist,lvec,up2n=2,nozero=.true.) + ! xp = (/0.3d0,0.1d0,0.4d0/) + ! call sy%c%list_near_lattice_points(xp,icrd_crys,.true.,nat,dist=dist,lvec=lvec,& + ! ndiv=n,up2n=10,nozero=.true.) ! do i = 1, nat ! write (*,*) i, lvec(:,i), dist(i) ! end do - call e%build_lattice(sy%c%m_x2c,150d0) - call e%list_near_atoms(xp,icrd_crys,.true.,nat,ierr,eid=eid,dist=dist,lvec=lvec2,up2n=1,nozero=.true.) - do i = 1, nat - write (*,*) i, e%xr2x(e%at(eid(i))%x) + lvec2, dist(i) - end do - write (*,*) - - call sy%c%nearest_lattice_point(xp,icrd_crys,dist(1),lvec2,nozero=.true.) - do i = 1, 1 - write (*,*) i, lvec2, dist(i) - end do - write (*,*) "fin!" stop 1