From 3879efed86b5dd78fa5672b9b24c296db63b5776 Mon Sep 17 00:00:00 2001 From: Alberto Otero de la Roza Date: Thu, 14 Mar 2024 14:30:59 +0100 Subject: [PATCH] remove calculating the rnn2 from the crystal input --- src/bisect@proc.f90 | 9 +++++---- src/crystalmod.f90 | 10 ++++++++-- src/crystalmod@complex.f90 | 4 +++- src/crystalmod@env.f90 | 15 +++++++++++++++ src/crystalmod@proc.f90 | 15 --------------- src/crystalmod@symmetry.f90 | 1 - src/crystalmod@write.f90 | 13 ++++++++----- src/qtree@proc.f90 | 22 +++++++++++++--------- src/types.f90 | 1 - 9 files changed, 52 insertions(+), 38 deletions(-) diff --git a/src/bisect@proc.f90 b/src/bisect@proc.f90 index 23e4c11b1..94718e874 100644 --- a/src/bisect@proc.f90 +++ b/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 diff --git a/src/crystalmod.f90 b/src/crystalmod.f90 index fca77991b..1fc117575 100644 --- a/src/crystalmod.f90 +++ b/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 diff --git a/src/crystalmod@complex.f90 b/src/crystalmod@complex.f90 index 4200ce27c..0a8a7dfd9 100644 --- a/src/crystalmod@complex.f90 +++ b/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 diff --git a/src/crystalmod@env.f90 b/src/crystalmod@env.f90 index a8a7ad585..299c340c9 100644 --- a/src/crystalmod@env.f90 +++ b/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, diff --git a/src/crystalmod@proc.f90 b/src/crystalmod@proc.f90 index 94ce06b91..27ae1503d 100644 --- a/src/crystalmod@proc.f90 +++ b/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 diff --git a/src/crystalmod@symmetry.f90 b/src/crystalmod@symmetry.f90 index 3fa1dab02..04d1f85d4 100644 --- a/src/crystalmod@symmetry.f90 +++ b/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 diff --git a/src/crystalmod@write.f90 b/src/crystalmod@write.f90 index f9c7c4046..65bf4f2f7 100644 --- a/src/crystalmod@write.f90 +++ b/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,7 +611,7 @@ 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 @@ -618,6 +619,7 @@ module subroutine struct_write_json(c,json,p) 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 diff --git a/src/qtree@proc.f90 b/src/qtree@proc.f90 index 9f831c748..1a70ea2a7 100644 --- a/src/qtree@proc.f90 +++ b/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), & diff --git a/src/types.f90 b/src/types.f90 index 942c7d287..973ef4920 100644 --- a/src/types.f90 +++ b/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