remove calculating the rnn2 from the crystal input

src/bisect@proc.f90

@@ -583,7 +583,7 @@ module subroutine sphereintegrals(line)
 
     character(len=:), allocatable :: word
     integer :: meth, nr, ntheta, nphi, np, cpid
-    real*8 :: r0, rend
+    real*8 :: r0, rend, rnn2
     integer :: linmin, linmax
     integer :: lp, i, j, n, nn, k
     real*8 :: sprop(sy%npropi)
@@ -724,10 +724,11 @@ module subroutine sphereintegrals(line)
 
     do i = linmin, linmax
        if ((sy%f(sy%iref)%cp(i)%typ /= sy%f(sy%iref)%typnuc .and. i>sy%c%nneq)) cycle
+       rnn2 = sy%c%get_rnn2(i)
 
        if (nr > 1) then
           if (rend < 0d0) then
-             h = log(sy%c%at(i)%rnn2 * abs(rend) / r0) / (nr - 1)
+             h = log(rnn2 * abs(rend) / r0) / (nr - 1)
           else
              h = log(rend / r0) / (nr - 1)
           end if
@@ -742,7 +743,7 @@ module subroutine sphereintegrals(line)
           string(r0,'e',decimal=6)
        if (rend < 0d0) then
           write (uout,'("  Final radius (rend,",A,"): ",A)') iunitname0(iunit), &
-             string(sy%c%at(i)%rnn2 * abs(rend),'e',decimal=6)
+             string(rnn2 * abs(rend),'e',decimal=6)
        else
           write (uout,'("  Final radius (rend,",A,"): ",A)') iunitname0(iunit), string(rend,'e',decimal=6)
        end if
@@ -755,7 +756,7 @@ module subroutine sphereintegrals(line)
        do n = 1, nr
           if (nr == 1) then
              if (rend < 0d0) then
-                r = sy%c%at(i)%rnn2 * abs(rend)
+                r = rnn2 * abs(rend)
              else
                 r = r0
              end if

src/crystalmod.f90

@@ -170,6 +170,7 @@ module crystalmod
      procedure :: build_env
      procedure :: list_near_atoms
      procedure :: nearest_atom
+     procedure :: get_rnn2
      procedure :: identify_atom
      procedure :: promolecular_atom
      procedure :: find_asterisms_covalent
@@ -398,6 +399,11 @@ module subroutine nearest_atom(c,xp,icrd,nid,dist,distmax,lvec,nid0,id0,iz0,ispc
        integer, intent(in), optional :: ispc0
        logical, intent(in), optional :: nozero
      end subroutine nearest_atom
+     module function get_rnn2(c,ineq)
+       class(crystal), intent(inout) :: c
+       integer, intent(in) :: ineq
+       real*8 :: get_rnn2
+     end function get_rnn2
      module function identify_atom(c,x0,icrd,lvec,dist,distmax)
        class(crystal), intent(inout) :: c
        real*8, intent(in) :: x0(3)
@@ -675,7 +681,7 @@ module subroutine amd(c,imax,res)
        real*8, intent(out) :: res(imax)
      end subroutine amd
      module subroutine struct_report(c,lcrys,lq)
-       class(crystal), intent(in) :: c
+       class(crystal), intent(inout) :: c
        logical, intent(in) :: lcrys
        logical, intent(in) :: lq
      end subroutine struct_report
@@ -689,7 +695,7 @@ module subroutine struct_report_symxyz(c,strfin,doaxes)
      end subroutine struct_report_symxyz
      module subroutine struct_write_json(c,json,p)
        use json_module, only: json_value, json_core
-       class(crystal), intent(in) :: c
+       class(crystal), intent(inout) :: c
        type(json_core), intent(inout) :: json
        type(json_value), pointer, intent(inout) :: p
      end subroutine struct_write_json

src/crystalmod@complex.f90

@@ -286,10 +286,12 @@ module function get_pack_ratio(c) result(px)
     real*8 :: px
 
     integer :: i
+    real*8 :: rnn2
 
     px = 0d0
     do i = 1, c%nneq
-       px = px + c%at(i)%mult * 4d0/3d0 * pi * c%at(i)%rnn2**3
+       rnn2 = c%get_rnn2(i)
+       px = px + c%at(i)%mult * 4d0/3d0 * pi * rnn2**3
     end do
     px = px / c%omega * 100d0
 

src/crystalmod@env.f90

@@ -555,6 +555,21 @@ module subroutine nearest_atom(c,xp,icrd,nid,dist,distmax,lvec,nid0,id0,iz0,ispc
 
   end subroutine nearest_atom
 
+  !> Returns half of the nearest neighbor distance for non-equivalent
+  !> atom ineq.
+  module function get_rnn2(c,ineq)
+    use param, only: icrd_cart
+    class(crystal), intent(inout) :: c
+    integer, intent(in) :: ineq
+    real*8 :: get_rnn2
+
+    integer :: iaux
+
+    call c%nearest_atom(c%at(ineq)%r,icrd_cart,iaux,get_rnn2,nozero=.true.)
+    get_rnn2 = 0.5d0 * get_rnn2
+
+  end function get_rnn2
+
   !> Given point x0 (with icrd input coordinates), if x0 corresponds
   !> to an atomic position (to within distmax or atomeps if distmax is
   !> not given), return the complete-list ID of the atom. Otherwise,

src/crystalmod@proc.f90

@@ -83,11 +83,6 @@ module subroutine struct_init(c)
        c%spc(i)%qat = 0d0
     end do
 
-    ! initialize atoms
-    do i = 1, mneq0
-       c%at(i)%rnn2 = 0d0
-    end do
-
     ! no molecular fragments
     c%nmol = 0
     if (allocated(c%nstar)) deallocate(c%nstar)
@@ -539,16 +534,6 @@ module subroutine struct_new(c,seed,crashfail,noenv,ti)
           end do
        end do
        deallocate(xcoord,iord)
-
-       ! Write the half nearest-neighbor distance
-       do i = 1, c%nneq
-          if (.not.c%ismolecule .or. c%ncel > 1) then
-             call c%nearest_atom(c%at(i)%r,icrd_cart,iat,dist,nozero=.true.)
-             c%at(i)%rnn2 = 0.5d0 * dist
-          else
-             c%at(i)%rnn2 = 0d0
-          end if
-       end do
     end if
 
     ! the initialization is done - this crystal is ready to use

src/crystalmod@symmetry.f90

@@ -409,7 +409,6 @@ module subroutine calcsym(c,usenneq,errmsg,ti)
           c%at(c%nneq)%r = c%atcel(idx)%r
           c%at(c%nneq)%is = c%atcel(idx)%is
           c%at(c%nneq)%mult = 1
-          c%at(c%nneq)%rnn2 = 0d0
           c%at(c%nneq)%wyc = "?"
        else
           c%at(iidx(idx))%mult = c%at(iidx(idx))%mult + 1

src/crystalmod@write.f90

@@ -30,14 +30,14 @@ module subroutine struct_report(c,lcrys,lq)
     use tools_math, only: gcd
     use tools_io, only: uout, string, ioj_center, ioj_left, ioj_right
     use param, only: bohrtoa, maxzat, pi, atmass, pcamu, bohrtocm
-    class(crystal), intent(in) :: c
+    class(crystal), intent(inout) :: c
     logical, intent(in) :: lcrys
     logical, intent(in) :: lq
 
     integer :: i, j, k, iz, is
     integer :: nelec
     real*8 :: maxdv, xcm(3), x0(3), xlen(3), xang(3), xred(3,3)
-    real*8 :: dens, mass
+    real*8 :: dens, mass, rnn2
     character(len=:), allocatable :: str1
     integer, allocatable :: nis(:)
     integer :: izp0
@@ -213,11 +213,12 @@ module subroutine struct_report(c,lcrys,lq)
           else
              str1 = ""
           end if
+          rnn2 = c%get_rnn2(c%atcel(i)%idx)
           write (uout,'("  ",99(A," "))') &
              string(i,3,ioj_center),&
              (string((c%atcel(i)%r(j)+c%molx0(j))*dunit0(iunit),'f',length=16,decimal=10,justify=5),j=1,3),&
              string(is,3,ioj_center),string(c%spc(is)%name,7,ioj_center), string(c%spc(is)%z,3,ioj_center),&
-             string(2d0*c%at(c%atcel(i)%idx)%rnn2*dunit0(iunit),'f',length=10,decimal=4,justify=4), &
+             string(2d0*rnn2*dunit0(iunit),'f',length=10,decimal=4,justify=4), &
              string(c%atcel(i)%idx,3,ioj_center), str1
        enddo
        write (uout,*)
@@ -610,14 +611,15 @@ end subroutine struct_report_symxyz
   !> json with root p.
   module subroutine struct_write_json(c,json,p)
     use json_module, only: json_value, json_core
-    class(crystal), intent(in) :: c
+    class(crystal), intent(inout) :: c
     type(json_core), intent(inout) :: json
     type(json_value), pointer, intent(inout) :: p
 
     type(json_value), pointer :: s, ap, arr
 
     integer :: i
     character(len=mlen), allocatable :: strfin(:)
+    real*8 :: rnn2
 
     if (.not.c%isinit) return
     call json%create_object(s,'structure')
@@ -658,7 +660,8 @@ module subroutine struct_write_json(c,json,p)
        call json%add(ap,'cartesian_coordinates',c%at(i)%r(:))
        call json%add(ap,'multiplicity',c%at(i)%mult)
        call json%add(ap,'wyckoff_letter',c%at(i)%wyc)
-       call json%add(ap,'half_nn_distance',c%at(i)%rnn2)
+       rnn2 = c%get_rnn2(i)
+       call json%add(ap,'half_nn_distance',rnn2)
        call json%add(arr,ap)
        nullify(ap)
     end do

src/qtree@proc.f90

@@ -916,7 +916,7 @@ subroutine qtree_setsph1(lvl,verbose)
     integer :: l2, i, j, k, l, tt, tto, idum
     integer :: nid, nid0, vin(3)
     integer(qtreeidx) :: idx
-    real*8 :: xx(3), dist, mdist, rmin, rmax, rdist, rmt
+    real*8 :: xx(3), dist, mdist, rmin, rmax, rdist, rmt, rnn2
     real*8, allocatable :: rref(:)
     logical, allocatable :: nucmask(:), nfrozen(:)
     integer :: niter, icolor
@@ -948,16 +948,17 @@ subroutine qtree_setsph1(lvl,verbose)
     ! Calculate sphere sizes (r_betagp and r_betaint)
     nfrozen = .false.
     do i = 1, nnuc
+       rnn2 = sy%c%get_rnn2(i)
        if (i<=sy%c%nneq .and. sy%f(sy%iref)%type == type_elk) then
           rmt = sy%f(sy%iref)%elk%rmt(sy%c%at(i)%is)
        elseif (i<=sy%c%nneq .and. sy%f(sy%iref)%type == type_wien) then
           rmt = sy%f(sy%iref)%wien%rmt_atom(sy%c%at(i)%x)
        else
-          rmt = sy%c%at(i)%rnn2
+          rmt = rnn2
        end if
 
        if (i <= sy%c%nneq) then
-          rref(i) = sy%c%at(i)%rnn2
+          rref(i) = rnn2
        else
           xx = sy%f(sy%iref)%cp(i)%x
           call sy%f(sy%iref)%nearest_cp(xx,idum,rref(i),type=sy%f(sy%iref)%typnuc,nozero=.true.)
@@ -1074,9 +1075,10 @@ subroutine qtree_setsph1(lvl,verbose)
     do i = 1, nnuc
        r_betaint(i) = sphintfactor(i) * r_betagp(i)
        if (verbose) then
+          rnn2 = sy%c%get_rnn2(i)
           write (uout,'("+ Sphfactor/rbeta/rnn2 of nuc ",A," (",A,") :",3("  ",A))') &
              string(i), string(sy%c%spc(sy%c%at(i)%is)%name), string(sphfactor(i),'f',decimal=6), &
-             string(r_betagp(i),'f',decimal=6), string(sy%c%at(i)%rnn2,'f',decimal=6)
+             string(r_betagp(i),'f',decimal=6), string(rnn2,'f',decimal=6)
        end if
     end do
     if (verbose) then
@@ -1114,7 +1116,7 @@ subroutine qtree_setsph2(verbose)
 
     type(minisurf) :: srf
     integer :: i, j, idum, ii
-    real*8 :: xx(3), dist, plen
+    real*8 :: xx(3), dist, plen, rnn2
     real*8 :: rref, x0(3), unit(3), rmt
     integer :: ier
     integer :: ndo, nstep
@@ -1138,18 +1140,19 @@ subroutine qtree_setsph2(verbose)
     ! Calculate sphere sizes (r_betagp and r_betaint)
     ndo = 0
     do i = 1, nnuc
+       rnn2 = sy%c%get_rnn2(i)
        if (i<=sy%c%nneq .and. sy%f(sy%iref)%type == type_elk) then
           rmt = sy%f(sy%iref)%elk%rmt(sy%c%at(i)%is)
        elseif (i<=sy%c%nneq .and. sy%f(sy%iref)%type == type_wien) then
           rmt = sy%f(sy%iref)%wien%rmt_atom(sy%c%at(i)%x)
        elseif (i<=sy%c%nneq) then
-          rmt = sy%c%at(i)%rnn2
+          rmt = rnn2
        else
           rmt = 1d30
        end if
 
        if (i <= sy%c%nneq) then
-          rref = sy%c%at(i)%rnn2
+          rref = rnn2
        else
           xx = sy%f(sy%iref)%cp(i)%x
           call sy%f(sy%iref)%nearest_cp(xx,idum,rref,type=sy%f(sy%iref)%typnuc,nozero=.true.)
@@ -1179,7 +1182,7 @@ subroutine qtree_setsph2(verbose)
        i = ido(ii)
        if (i <= sy%c%nneq) then
           x0 = sy%c%at(i)%r
-          rref = sy%c%at(i)%rnn2
+          rref = sy%c%get_rnn2(i)
        else
           x0 = sy%f(sy%iref)%cp(i)%r
           call sy%f(sy%iref)%nearest_cp(xx,idum,rref,type=sy%f(sy%iref)%typnuc,nozero=.true.)
@@ -1217,9 +1220,10 @@ subroutine qtree_setsph2(verbose)
        r_betaint(i) = sphintfactor(i) * r_betagp(i)
        if (verbose) then
           if (i <= sy%c%nneq) then
+             rnn2 = sy%c%get_rnn2(i)
              write (uout,'("+ Sphfactor/rbeta/rnn2 of nuc ",A," (",A,") :",3(A,"  "))') &
                 string(i), string(sy%c%spc(sy%c%at(i)%is)%name), string(sphfactor(i),'f',decimal=6), &
-                string(r_betagp(i),'f',decimal=6), string(sy%c%at(i)%rnn2,'f',decimal=6)
+                string(r_betagp(i),'f',decimal=6), string(rnn2,'f',decimal=6)
           else
              write (uout,'("+ Sphfactor/rbeta/rnn2 of nuc ",A," (",A,") :",3(A,"  "))') &
                 string(i), "nnm", string(sphfactor(i),'f',decimal=6), &

src/types.f90

@@ -103,7 +103,6 @@ module types
      integer :: is = 0 !< species
      integer :: mult  !< multiplicity
      character*1 :: wyc !< Wyckoff letter
-     real*8 :: rnn2   !< half the nearest neighbor distance
   end type neqatom
 
   !> Any atom in the crystal type