From 2841f92b6abdd5fbaa5dd6c49b20eb9af5b2fe02 Mon Sep 17 00:00:00 2001 From: Alberto Otero de la Roza Date: Sun, 4 Feb 2024 10:32:14 +0100 Subject: [PATCH] initial implementation of simpler atomic environments --- src/CMakeLists.txt | 3 +- src/crystalmod.f90 | 34 ++++++ src/crystalmod@env.f90 | 238 +++++++++++++++++++++++++++++++++++++++++ src/global@proc.F90 | 30 ++++++ src/types.f90 | 1 + 5 files changed, 305 insertions(+), 1 deletion(-) create mode 100644 src/crystalmod@env.f90 diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 7e708e96..7db2a8be 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -2,7 +2,8 @@ set(SOURCES_Fortran abinit_private.f90 abinit_private@proc.f90 autocp.f90 autocp@proc.f90 bader.f90 bader@proc.f90 bisect.f90 bisect@proc.f90 cfftnd.f90 c_interface_module.f90 config.f90 - crystalmod.f90 crystalmod@complex.f90 crystalmod@edit.f90 crystalmod@mols.f90 + crystalmod.f90 crystalmod@complex.f90 crystalmod@edit.f90 crystalmod@env.f90 + crystalmod@mols.f90 crystalmod@proc.f90 crystalmod@powderproc.f90 crystalmod@symmetry.f90 crystalmod@write.f90 crystalseedmod.f90 crystalseedmod@proc.f90 dftb_private.f90 diff --git a/src/crystalmod.f90 b/src/crystalmod.f90 index 39893888..c85dfae9 100644 --- a/src/crystalmod.f90 +++ b/src/crystalmod.f90 @@ -132,6 +132,11 @@ module crystalmod !! Initialization level: isenv !! ! atomic environment of the cell type(environ) :: env + integer :: nblock(3) ! number of environemt blocks + integer,allocatable :: iblock0(:,:,:) ! starting atomic index for each block + real*8 :: blockrmax ! radius of the largest sphere contained in a block + real*8 :: blockomega ! volume of a block + real*8 :: blockcv(3) ! cross products of the block lattice vectors !! Initialization level: isast !! ! asterisms @@ -176,6 +181,10 @@ module crystalmod procedure :: identify_atom !< Identify an atom in the unit cell procedure :: identify_spc !< Identify species in a structure from a string + ! atomic environments and distance calculations + procedure :: build_env + procedure :: list_near_atoms + ! molecular environments and neighbors (mols) procedure :: identify_fragment !< Build an atomic fragment of the crystal procedure :: identify_fragment_from_xyz !< Build a crystal fragment from an xyz file @@ -382,6 +391,31 @@ module function identify_spc(c,str) result(res) character*(*), intent(in) :: str integer :: res end function identify_spc + module subroutine build_env(c) + class(crystal), intent(inout) :: c + end subroutine build_env + module subroutine list_near_atoms(c,xp,icrd,sorted,nat,eid,dist,lvec,ishell0,up2d,& + up2dsp,up2dcidx,up2sh,up2n,nid0,id0,iz0,ispc0,nozero) + class(crystal), intent(inout) :: c + real*8, intent(in) :: xp(3) + integer, intent(in) :: icrd + logical, intent(in) :: sorted + integer, intent(out) :: nat + integer, allocatable, intent(inout), optional :: eid(:) + real*8, allocatable, intent(inout), optional :: dist(:) + integer, allocatable, intent(inout), optional :: lvec(:,:) + integer, allocatable, intent(inout), optional :: ishell0(:) + real*8, intent(in), optional :: up2d + real*8, intent(in), optional :: up2dsp(:,:) + real*8, intent(in), optional :: up2dcidx(:) + integer, intent(in), optional :: up2sh + integer, intent(in), optional :: up2n + integer, intent(in), optional :: nid0 + integer, intent(in), optional :: id0 + integer, intent(in), optional :: iz0 + integer, intent(in), optional :: ispc0 + logical, intent(in), optional :: nozero + end subroutine list_near_atoms 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 new file mode 100644 index 00000000..14d0df0b --- /dev/null +++ b/src/crystalmod@env.f90 @@ -0,0 +1,238 @@ +! Copyright (c) 2007-2022 Alberto Otero de la Roza , +! Ángel Martín Pendás and Víctor Luaña +! . +! +! critic2 is free software: you can redistribute it and/or modify +! it under the terms of the GNU General Public License as published by +! the Free Software Foundation, either version 3 of the License, or (at +! your option) any later version. +! +! critic2 is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU General Public License for more details. +! +! You should have received a copy of the GNU General Public License +! along with this program. If not, see . + +! Routines for atomic environments and distance calculations. +submodule (crystalmod) env + implicit none + +contains + + !> Fill the environment variables (nblock, iblock0, blockrmax) as + !> well as the %inext field in the complete atom list + !> atoms. Requires the reduced lattice vectors as well as the + !> complete atom list. + module subroutine build_env(c) + use tools_math, only: det3, cross + use param, only: pi + class(crystal), intent(inout) :: c + + integer :: i, ix(3) + real*8 :: xred(3,3), xlen(3), aux, xr(3) + + real*8, parameter :: lenmax = 15d0 + + ! calculate the number of cells in each direction and the lattice vectors + do i = 1, 3 + c%nblock(i) = ceiling(norm2(c%m_xr2c(:,i))/lenmax) + xred(:,i) = c%m_xr2c(:,i) / c%nblock(i) + end do + + ! block volume and block radius + c%blockomega = det3(xred) + c%blockrmax = 0.5d0 * c%blockomega / max(norm2(cross(xred(:,1),xred(:,2))),& + norm2(cross(xred(:,1),xred(:,3))),norm2(cross(xred(:,2),xred(:,3)))) + do i = 1, 3 + xlen(i) = norm2(xred(:,i)) + end do + c%blockcv(1) = norm2(cross(xred(:,2),xred(:,3))) + c%blockcv(2) = norm2(cross(xred(:,1),xred(:,3))) + c%blockcv(3) = norm2(cross(xred(:,1),xred(:,2))) + + ! assign atomic indices + if (allocated(c%iblock0)) deallocate(c%iblock0) + allocate(c%iblock0(0:c%nblock(1)-1,0:c%nblock(2)-1,0:c%nblock(3)-1)) + c%iblock0 = 0 + c%atcel(:)%inext = 0 + do i = 1, c%ncel + xr = c%x2xr(c%atcel(i)%x) + xr = xr - floor(xr) + ix = min(max(floor(xr * c%nblock),0),c%nblock-1) + c%atcel(i)%inext = c%iblock0(ix(1),ix(2),ix(3)) + c%iblock0(ix(1),ix(2),ix(3)) = i + end do + + end subroutine build_env + + !> Given the point xp (in icrd coordinates), calculate the list of + !> nearest atoms. If sorted is true, the list is sorted by distance + !> and shell. One of these optional criteria must be given: + !> - up2d = list up to a distance up2d. + !> - up2dsp = between a minimum and a maximum species-dependent distance. + !> - up2dcidx = up to a distance to each atom in the complete list. + !> - up2sh = up to a number of shells. + !> - up2n = up to a number of atoms. + !> If up2n or up2sh is used, the list is always sorted. + !> + !> Optional input: + !> - nid0 = consider only atoms with index nid0 from the nneq list. + !> - id0 = consider only atoms with index id0 from the complete list. + !> - iz0 = consider only atoms with atomic number iz0. + !> - ispc0 = consider only species with species ID ispc0. + !> - nozero = disregard zero-distance atoms. + !> + !> Output: + !> - nat = the list contains nat atoms. + !> + !> Optional output: + !> - eid(nid) = atom ID from the complete atom list. + !> - dist(nid) = distance from xp. + !> - lvec(3,nid) = lattice vectors for the atoms. + !> - ishell0(nid) = shell ID. + !> + !> Sorted and the up2sh and up2n options work only if the dist(:) + !> array is present. + module subroutine list_near_atoms(c,xp,icrd,sorted,nat,eid,dist,lvec,ishell0,up2d,& + up2dsp,up2dcidx,up2sh,up2n,nid0,id0,iz0,ispc0,nozero) + use tools, only: mergesort + use tools_io, only: ferror, faterr + use tools_math, only: cross, det3 + use types, only: realloc + use param, only: icrd_cart, icrd_crys, icrd_rcrys + class(crystal), intent(inout) :: c + real*8, intent(in) :: xp(3) + integer, intent(in) :: icrd + logical, intent(in) :: sorted + integer, intent(out) :: nat + integer, allocatable, intent(inout), optional :: eid(:) + real*8, allocatable, intent(inout), optional :: dist(:) + integer, allocatable, intent(inout), optional :: lvec(:,:) + integer, allocatable, intent(inout), optional :: ishell0(:) + real*8, intent(in), optional :: up2d + real*8, intent(in), optional :: up2dsp(:,:) + real*8, intent(in), optional :: up2dcidx(:) + integer, intent(in), optional :: up2sh + integer, intent(in), optional :: up2n + integer, intent(in), optional :: nid0 + integer, intent(in), optional :: id0 + integer, intent(in), optional :: iz0 + integer, intent(in), optional :: ispc0 + logical, intent(in), optional :: nozero + + logical :: doshell + real*8 :: x(3), xorigc(3), dmax, dd, lvecx(3), xr(3) + integer :: i, j, k, ix(3), nx(3), i0(3), i1(3), idx + integer :: ib(3) + integer, allocatable :: at_id(:), at_lvec(:,:), iord(:) + real*8, allocatable :: at_dist(:) + + ! process input + if (.not.present(up2d).and..not.present(up2dsp).and..not.present(up2dcidx).and.& + .not.present(up2sh).and..not.present(up2n)) & + call ferror("list_near_atoms","must give one of up2d, up2dsp, up2dcidx, up2sh, or up2n",faterr) + doshell = present(ishell0) .or. present(up2sh) + + ! locate the block containing the point + if (icrd == icrd_crys) then + x = c%x2xr(xp) + elseif (icrd == icrd_cart) then + x = c%c2xr(xp) + else + x = xp + end if + xorigc = c%xr2c(x) + ix = floor(x * c%nblock) + + ! up to a distance + if (present(up2d)) then + dmax = up2d + else + write (*,*) "fixme!" + stop 1 + end if + + ! 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) + nx = ceiling(c%blockcv / (0.5d0 * c%blockomega / max(dmax,1d-40))) + + ! define the search space + do i = 1, 3 + if (mod(nx(i),2) == 0) then + i0(i) = ix(i) - nx(i)/2 + i1(i) = ix(i) + nx(i)/2 + else + i0(i) = ix(i) - (nx(i)/2+1) + i1(i) = ix(i) + (nx(i)/2+1) + end if + end do + + ! allocate space for atoms + nat = 0 + allocate(at_id(20),at_dist(20),at_lvec(3,20)) + + ! run over the atoms and compile the list + do i = i0(1), i1(1) + do j = i0(2), i1(2) + do k = i0(3), i1(3) + ib = (/modulo(i,c%nblock(1)), modulo(j,c%nblock(2)), modulo(k,c%nblock(3))/) + lvecx = real(((/i,j,k/) - ib) / c%nblock,8) + + idx = c%iblock0(ib(1),ib(2),ib(3)) + do while (idx /= 0) + xr = c%x2xr(c%atcel(idx)%x) + xr = xr - floor(xr) + lvecx + x = c%xr2c(xr) + dd = norm2(x - xorigc) + if (dd <= dmax) then + call add_atom_to_output_list() + end if + idx = c%atcel(idx)%inext + end do + end do + end do + end do + + ! sort if necessary + if (sorted .and. nat > 0) then + ! permutation vector + allocate(iord(nat)) + do i = 1, nat + iord(i) = i + end do + + ! sort by distance + call mergesort(at_dist,iord,1,nat) + + ! reorder and clean up + at_id = at_id(iord) + at_dist = at_dist(iord) + at_lvec = at_lvec(:,iord) + deallocate(iord) + end if + + ! assign optional outputs + if (present(eid)) eid = at_id + if (present(dist)) dist = at_dist + if (present(lvec)) lvec = at_lvec + + contains + subroutine add_atom_to_output_list() + nat = nat + 1 + if (nat > size(at_id,1)) then + call realloc(at_id,2*nat) + call realloc(at_dist,2*nat) + call realloc(at_lvec,3,2*nat) + end if + at_id(nat) = idx + at_dist(nat) = dd + at_lvec(:,nat) = nint(c%xr2x(xr) - c%atcel(idx)%x) + end subroutine add_atom_to_output_list + + end subroutine list_near_atoms + + +end submodule env diff --git a/src/global@proc.F90 b/src/global@proc.F90 index ef0dfd25..d91c5b9c 100644 --- a/src/global@proc.F90 +++ b/src/global@proc.F90 @@ -870,6 +870,8 @@ module subroutine critic_setvariables(line,lp) use tools_io, only: lgetword, getword, equal, isinteger, isreal, ferror, & faterr, string, uout, isassignment, getword, zatguess use param, only: maxzat0, atmcov, atmvdw + use systemmod, only: sy ! xxxx + use param, only: icrd_crys, icrd_rcrys ! xxxx character*(*), intent(in) :: line integer, intent(inout) :: lp @@ -879,6 +881,10 @@ module subroutine critic_setvariables(line,lp) integer :: lp2, iz logical :: iscov character(len=:), allocatable :: errmsg + real*8 :: x(3), xx(3), dmax ! xxxx + integer :: j, nat, ierr, lvec(3), nat1, nat2 ! xxxx + integer, allocatable :: eid(:), lvec2(:,:) ! xxxx + real*8, allocatable :: dist(:) ! xxxx word = lgetword(line,lp) if (equal(word,'bondfactor')) then @@ -1064,6 +1070,30 @@ module subroutine critic_setvariables(line,lp) atmvdw(iz) = rdum end if end do + elseif (equal(word,'temp')) then ! xxxx + call sy%c%build_env() + + dmax = 20d0 + x = (/-10.0d0,20.2d0,30.3d0/) + ! x = 0.5d0 + ! x = (/0.9d0,0.8d0,0.3d0/) + call sy%c%env%list_near_atoms(x,icrd_rcrys,.true.,nat,ierr,eid,dist,lvec,up2d=dmax) + do j = 1, nat + xx = sy%c%env%at(eid(j))%x + sy%c%env%x2xr(real(lvec,8)) + write (*,*) eid(j), xx, dist(j) + end do + ! write (*,*) "xx1 ", nat + write (*,*) + nat1 = nat + + call sy%c%list_near_atoms(x,icrd_rcrys,.true.,nat,eid,dist,lvec2,up2d=dmax) + do j = 1, nat + write (*,*) eid(j), sy%c%x2xr(sy%c%atcel(eid(j))%x + lvec2(:,j)), dist(j) + end do + ! write (*,*) "xx2 ", nat + nat2 = nat + ! write (*,*) "xx ", nat1, nat2 + elseif (isassignment(var,word,line)) then rdum = eval(word,errmsg) if (len_trim(errmsg) > 0) then diff --git a/src/types.f90 b/src/types.f90 index 742c1ae7..62a12c91 100644 --- a/src/types.f90 +++ b/src/types.f90 @@ -125,6 +125,7 @@ module types integer :: lvec(3) !< lattice vector to the atom in the complete atom list integer :: ir !< rotation matrix to the representative equivalent atom integer :: ic !< translation vector to the representative equivalent atom + integer :: inext = 0 !< next atom in the atomic environment block end type celatom !> Result of the evaluation of a scalar field