Skip to content

Commit

Permalink
implemented ndiv in list_near_lattice_points
Browse files Browse the repository at this point in the history
  • Loading branch information
@aoterodelaroza
aoterodelaroza committed Feb 13, 2024
1 parent ecf6e5e commit 1c59514
Show file tree
Hide file tree
Showing 3 changed files with 65 additions and 36 deletions.
4 changes: 3 additions & 1 deletion src/crystalmod.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -460,14 +460,16 @@ 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
logical, intent(in) :: sorted
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
Expand Down
50 changes: 34 additions & 16 deletions src/[email protected]
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand Down
47 changes: 28 additions & 19 deletions src/[email protected]
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down

0 comments on commit 1c59514

Please sign in to comment.