From ecf6e5e530359957e20e72d3bcd7852ae5d62685 Mon Sep 17 00:00:00 2001
From: Alberto Otero de la Roza <aoterodelaroza@gmail.com>
Date: Tue, 13 Feb 2024 16:03:31 +0100
Subject: [PATCH] implemented nearest_lattice_point

---
 src/crystalmod.f90     |  9 +++++++++
 src/crystalmod@env.f90 | 31 +++++++++++++++++++++++++++++++
 src/global@proc.F90    | 24 ++++++++++++++++++------
 3 files changed, 58 insertions(+), 6 deletions(-)

diff --git a/src/crystalmod.f90 b/src/crystalmod.f90
index aebb703d..f615adfa 100644
--- a/src/crystalmod.f90
+++ b/src/crystalmod.f90
@@ -189,6 +189,7 @@ module crystalmod
      procedure :: promolecular_env
      procedure :: find_asterisms_covalent
      procedure :: list_near_lattice_points
+     procedure :: nearest_lattice_point
 
      ! molecular environments and neighbors (mols)
      procedure :: identify_fragment !< Build an atomic fragment of the crystal
@@ -471,6 +472,14 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u
        integer, intent(in), optional :: up2n
        logical, intent(in), optional :: nozero
      end subroutine list_near_lattice_points
+     module subroutine nearest_lattice_point(c,xp,icrd,dist,lvec,nozero)
+       class(crystal), intent(inout) :: c
+       real*8, intent(in) :: xp(3)
+       integer, intent(in) :: icrd
+       real*8, intent(out) :: dist
+       integer, intent(out), optional :: lvec(3)
+       logical, intent(in), optional :: nozero
+     end subroutine nearest_lattice_point
      module function identify_fragment(c,nat,x0) result(fr)
        class(crystal), intent(in) :: c
        integer, intent(in) :: nat
diff --git a/src/crystalmod@env.f90 b/src/crystalmod@env.f90
index f557a799..35b7d4ce 100644
--- a/src/crystalmod@env.f90
+++ b/src/crystalmod@env.f90
@@ -1020,6 +1020,37 @@ module subroutine list_near_lattice_points(c,xp,icrd,sorted,nat,dist,lvec,up2d,u
 
   end subroutine list_near_lattice_points
 
+  !> Given the point xp (in icrd coordinates), calculate the nearest
+  !> lattice point. The lattice point (in cryst. coords.) is returne
+  !> in lvec and the distance in dist. If nozero, disregard
+  !> zero-distance lattice points. This routine is thread-safe.
+  module subroutine nearest_lattice_point(c,xp,icrd,dist,lvec,nozero)
+    class(crystal), intent(inout) :: c
+    real*8, intent(in) :: xp(3)
+    integer, intent(in) :: icrd
+    real*8, intent(out) :: dist
+    integer, intent(out), optional :: lvec(3)
+    logical, intent(in), optional :: nozero
+
+    integer :: nat
+    integer, allocatable :: eid(:), lvec_(:,:)
+    real*8, allocatable :: dist_(:)
+
+    ! get just one atom
+    call c%list_near_lattice_points(xp,icrd,.false.,nat,dist=dist_,lvec=lvec_,up2n=1,&
+       nozero=nozero)
+
+    ! if no atoms in output, return
+    if (present(lvec)) lvec = 0
+    dist = huge(1d0)
+    if (nat == 0) return
+
+    ! write and finish
+    dist = dist_(1)
+    if (present(lvec)) lvec = lvec_(:,1)
+
+  end subroutine nearest_lattice_point
+
   !xx! private procedures
 
   ! Find the indices for the nth shell of blocks. Sets the number of indices (nb),
diff --git a/src/global@proc.F90 b/src/global@proc.F90
index 8c7116d1..0ee083f9 100644
--- a/src/global@proc.F90
+++ b/src/global@proc.F90
@@ -634,20 +634,32 @@ module subroutine critic_main()
        elseif (equal(word,'temp')) then
           ! xxxx
 
-          ! xp = (/-15.3d0,0.1d0,21.4d0/)
-          xp = 0d0
+          xp = (/-15.3d0,0.1d0,21.4d0/)
+          ! xp = 0d0
           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.)
+          ! 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.)
+          ! 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=2,nozero=.true.)
+          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%list_near_lattice_points(xp,icrd_crys,.true.,nat,dist,lvec,up2n=2,nozero=.true.)
-          do i = 1, nat
-             write (*,*) i, lvec(:,i), dist(i)
+          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!"