TABLE OF CONTENTS
- ABINIT/m_bz_mesh
- m_bz_mesh/box_len
- m_bz_mesh/bz_mesh_isirred
- m_bz_mesh/find_qmesh
- m_bz_mesh/findnq
- m_bz_mesh/findq
- m_bz_mesh/findqg0
- m_bz_mesh/get_BZ_diff
- m_bz_mesh/get_bz_item
- m_bz_mesh/get_IBZ_item
- m_bz_mesh/get_ng0sh
- m_bz_mesh/getkptnorm_bycomponent
- m_bz_mesh/has_BZ_item
- m_bz_mesh/has_IBZ_item
- m_bz_mesh/identk
- m_bz_mesh/isamek
- m_bz_mesh/isequalk
- m_bz_mesh/kmesh_free
- m_bz_mesh/kmesh_init
- m_bz_mesh/kmesh_print
- m_bz_mesh/kmesh_t
- m_bz_mesh/kpath_free
- m_bz_mesh/kpath_new
- m_bz_mesh/kpath_print
- m_bz_mesh/kpath_t
- m_bz_mesh/littlegroup_free_0D
- m_bz_mesh/littlegroup_free_1D
- m_bz_mesh/littlegroup_init
- m_bz_mesh/littlegroup_print
- m_bz_mesh/littlegroup_t
- m_bz_mesh/make_mesh
- m_bz_mesh/make_path
- m_bz_mesh/setup_k_rotation
ABINIT/m_bz_mesh [ Modules ]
NAME
m_bz_mesh
FUNCTION
This module provides the definition of the kmesh_t structure gathering information on the sampling of the Brillouin zone. It also contains useful tools to operate on k-points. and the definition of the littlegroup_t data type. The littlegroup_t structure is used to store tables and useful info on the set of k-points belonging to the irreducible wedge defined by the symmetry properties of the point group that preserve the external q-point.
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, GMR, VO, LR, RWG, MT) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
NOTES
One has to use a fixed ordering of the loops over nsym and time-reversal when the full zone is reconstructed by symmetry starting from the IBZ. This is especially important in systems with both time-reversal and spatial inversion as the effect of the two operation in reciprocal space is very similar the only difference being the possibly non-zero fractional translation associated to the spatial inversion. In the present implementation, the spatial inversion has the precedence wrt time-reversal (i.e., do itim; do isym). Note that this particular ordering should be used in any routine used to symmetrize k-dependent quantities in the full BZ zone to avoid possible errors. * This module is deprecated and should be used only in the GW/BSE part. Some of the routines will be gradually moved to m_kpts
SOURCE
36 #if defined HAVE_CONFIG_H 37 #include "config.h" 38 #endif 39 40 #include "abi_common.h" 41 42 module m_bz_mesh 43 44 use defs_basis 45 use m_errors 46 use m_abicore 47 use m_sort 48 use m_xmpi 49 50 use m_fstrings, only : ltoa, itoa, sjoin, ktoa 51 use m_numeric_tools, only : is_zero, isinteger, imin_loc, imax_loc, bisect, wrap2_pmhalf 52 use m_symtk, only : sg_multable, littlegroup_q 53 use m_geometry, only : normv 54 use m_crystal, only : crystal_t 55 use m_kpts, only : getkgrid 56 use m_symkpt, only : symkpt 57 58 implicit none 59 60 private 61 62 real(dp),parameter :: TOL_KDIFF = 0.0001_dp 63 ! Tolerance below which two points are considered equal within a RL vector: 64 ! for each reduced direction the absolute difference between the coordinates must be less that TOL_KDIFF 65 66 integer,parameter :: NONE_KPTRLATT(3,3) = RESHAPE((/0,0,0,0,0,0,0,0,0/),(/3,3/))
m_bz_mesh/box_len [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
box_len
FUNCTION
Given a direction in q-space defined by the q-point qpt, this function returns the length of the vector connecting the origin with one the faces of the cell defined by the lattice vectors gprimd.
INPUTS
qpt(3)=The reduced coordinates of the q-point defining the direction. Normalization is not mandatory. gprimd(3,3)=Cartesian coordinates of the vectors defining the lattice.
SOURCE
2715 function box_len(qpt, gprimd) 2716 2717 !Arguments ------------------------------------ 2718 !scalars 2719 real(dp) :: box_len 2720 !arrays 2721 real(dp),intent(in) :: qpt(3),gprimd(3,3) 2722 2723 !Local variables------------------------------- 2724 !scalars 2725 integer :: idir,iplane 2726 real(dp) :: x1,x2,x3 2727 !arrays 2728 real(dp) :: my_qpt(3),gmet(3,3),q0box(3) 2729 2730 ! ************************************************************************* 2731 2732 ! Compute reciprocal space metric 2733 gmet = MATMUL(TRANSPOSE(gprimd),gprimd) 2734 2735 ! Rotate the input q-point such that it is always in the first octant then normalize it. 2736 ! Bravais lattices are invariant under inversion of any of the basis vectors. 2737 my_qpt = ABS(qpt)/normv(qpt,gmet,"G") 2738 2739 ! Check whether the q is along one of the reduced directions. 2740 idir=0; if (COUNT(ABS(qpt)<tol16) == 2) idir = imax_loc(ABS(qpt)) 2741 2742 if (idir/=0) then ! easy as q is along vector idir. 2743 box_len = normv(gprimd(:,idir), gmet, "G") 2744 RETURN 2745 2746 else 2747 !iplane/=0 means that q is placed on one the planes defined by two reciprocal lattice vectors. 2748 iplane=0; if (COUNT(ABS(qpt)<tol16) == 1) iplane = imin_loc(ABS(qpt)) 2749 q0box = [-1,-1,-1] 2750 2751 if (iplane/=1) then 2752 x1=one 2753 x2=my_qpt(2)/my_qpt(1) 2754 x3=my_qpt(3)/my_qpt(1) 2755 if (x2<=one+tol16 .and. x3<=one+tol16) q0box=(/x1,x2,x3/) 2756 end if 2757 if (iplane/=2) then 2758 x1=my_qpt(1)/my_qpt(2) 2759 x2=one 2760 x3=my_qpt(3)/my_qpt(2) 2761 if (x1<=one+tol16 .and. x3<=one+tol16) q0box=(/x1,x2,x3/) 2762 end if 2763 if (iplane/=3) then 2764 x1=my_qpt(1)/my_qpt(3) 2765 x2=my_qpt(2)/my_qpt(3) 2766 x3=one 2767 if (x1<=one+tol16 .and. x2<=one+tol16) q0box=(/x1,x2,x3/) 2768 end if 2769 2770 if (ALL(q0box == [-1,-1,-1])) then 2771 ABI_BUG("Cannot found q0box") 2772 end if 2773 2774 box_len = normv(q0box,gmet,"G") 2775 RETURN 2776 end if 2777 2778 end function box_len
m_bz_mesh/bz_mesh_isirred [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
bz_mesh_isirred
FUNCTION
bz_mesh_isirred
INPUTS
ik_bz=Index of the k-point in the BZ.
OUTPUT
bz_mesh_isirred=.TRUE. if the k-point is in the IBZ (a non-zero umklapp is not allowed)
SOURCE
1230 pure logical function bz_mesh_isirred(Kmesh, ik_bz) 1231 1232 !Arguments ------------------------------------ 1233 !scalars 1234 class(kmesh_t),intent(in) :: Kmesh 1235 integer,intent(in) :: ik_bz 1236 1237 !Local variables------------------------------- 1238 !scalars 1239 integer :: isym,itim 1240 1241 ! ********************************************************************* 1242 1243 isym = Kmesh%tabo(ik_bz) 1244 itim = (3-Kmesh%tabi(ik_bz))/2 1245 1246 ! Be careful here as we assume a particular ordering of symmetries. 1247 bz_mesh_isirred = ( isym==1.and.itim==1.and.ALL(Kmesh%umklp(:,ik_bz)==(/0,0,0/)) ) 1248 1249 end function bz_mesh_isirred
m_bz_mesh/find_qmesh [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
find_qmesh
FUNCTION
Find the q-mesh defined as all the possible differences between k-points Find the irreducible q-points using a special treatment for the Gamma point. Then call setup_kmesh to initialize the Qmesh datatype
INPUTS
Cryst<crystal_t>=datatype gathering info on the unit cell and symmetries %nsym=number of symmetry operations %symrec(3,3,nsym)=symmetry operations in reciprocal space Kmesh<kmesh_t>=datatype gathering information on the k-mesh
OUTPUT
Qmesh<kmesh_t>=datatype gathering information on the q point sampling.
SOURCE
1887 subroutine find_qmesh(Qmesh,Cryst,Kmesh) 1888 1889 !Arguments ------------------------------------ 1890 !scalars 1891 type(kmesh_t),intent(in) :: Kmesh 1892 type(kmesh_t),intent(inout) :: Qmesh 1893 type(crystal_t),intent(in) :: Cryst 1894 1895 !Local variables------------------------------- 1896 !scalars 1897 integer :: nqibz,kptopt 1898 !arrays 1899 real(dp),allocatable :: qibz(:,:) 1900 1901 ! ************************************************************************* 1902 1903 !@kmesh_t 1904 ! * Find the number of q-points such that q = k1-k2. 1905 call findnq(Kmesh%nbz,Kmesh%bz,Cryst%nsym,Cryst%symrec,Cryst%symafm,nqibz,Cryst%timrev) 1906 ! 1907 ! * Find the coordinates of the q-points in the IBZ. 1908 ABI_MALLOC(qibz,(3,nqibz)) 1909 call findq(Kmesh%nbz,Kmesh%bz,Cryst%nsym,Cryst%symrec,Cryst%symafm,Cryst%gprimd,nqibz,qibz,Cryst%timrev) 1910 ! 1911 ! * Create the Qmesh object starting from the IBZ. 1912 kptopt = Kmesh%kptopt 1913 call kmesh_init(Qmesh,Cryst,nqibz,qibz,kptopt) 1914 1915 ABI_FREE(qibz) 1916 1917 end subroutine find_qmesh
m_bz_mesh/findnq [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
findnq
FUNCTION
Identify the number of q-points in the IBZ by which the k-points in BZ differ (count the q points in the k-point difference set)
INPUTS
nkbz=number of k points in Brillouin zone kbz(3,nkbz)=coordinates of k points in BZ nsym=number of symmetry operations symrec(3,3,nsym)=symmetry operations in reciprocal space symafm(nsym)=-1 if AFM symmetry. timrev=2 if time-reversal symmetry is used, 1 otherwise
OUTPUT
nqibz=number of q points
SOURCE
1943 subroutine findnq(nkbz,kbz,nsym,symrec,symafm,nqibz,timrev) 1944 1945 !Arguments ------------------------------------ 1946 !scalars 1947 integer,intent(in) :: timrev,nkbz,nsym 1948 integer,intent(out) :: nqibz 1949 !arrays 1950 integer,intent(in) :: symrec(3,3,nsym),symafm(nsym) 1951 real(dp),intent(in) :: kbz(3,nkbz) 1952 1953 !Local variables ------------------------------ 1954 !scalars 1955 integer :: ifound,ik,isym,iq,memory_exhausted,nqall,nqallm,itim,ierr 1956 !arrays 1957 integer :: g0(3) 1958 real(dp) :: qposs(3),qrot(3) 1959 real(dp),allocatable :: qall(:,:) 1960 !************************************************************************ 1961 1962 ! Infinite do-loop to be able to allocate sufficient memory 1963 nqallm=1000 1964 do 1965 memory_exhausted=0 1966 ABI_MALLOC_OR_DIE(qall,(3,nqallm), ierr) 1967 nqall=0 1968 1969 ! Loop over all k-points in BZ, forming k-k1. 1970 do ik=1,nkbz 1971 qposs(:)=kbz(:,ik)-kbz(:,1) 1972 1973 ! Check whether this q (or its equivalent) has already been found within a reciprocal lattice vector. 1974 ! Use spatial inversion instead of time reversal whenever possible. 1975 ifound=0 1976 do iq=1,nqall 1977 do itim=1,timrev 1978 do isym=1,nsym 1979 if (symafm(isym)==-1) CYCLE 1980 qrot = (3-2*itim) * MATMUL(symrec(:,:,isym),qall(:,iq)) 1981 if (isamek(qrot,qposs,g0)) ifound=ifound+1 1982 end do 1983 end do 1984 end do 1985 1986 if (ifound==0) then 1987 nqall=nqall+1 1988 1989 ! If not yet found, check that the allocation is big enough. 1990 if (nqall>nqallm) then 1991 memory_exhausted=1 1992 ABI_FREE(qall) 1993 nqallm=nqallm*2; EXIT ! Exit the do ik=1 loop 1994 end if 1995 1996 ! Add it to the list. 1997 qall(:,nqall)=qposs(:) 1998 end if 1999 end do 2000 2001 if (memory_exhausted==0) EXIT 2002 end do !infinite loop 2003 2004 ABI_FREE(qall) 2005 nqibz=nqall 2006 2007 end subroutine findnq
m_bz_mesh/findq [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
findq
FUNCTION
Identify the q-points by which the k-points in BZ differ
INPUTS
nkbz=number of k points in Brillouin zone kbz(3,nkbz)=coordinates of k points in BZ nsym=number of symmetry operations symrec(3,3,nsym)=symmetry operations in reciprocal space symafm(nsym)=-1 is symmetry is AFM, +1 otherwise. gprimd(3,3)=dimensional reciprocal space primitive translations nqibz=number of q points in the IBZ by which k points differ (computed in findnq) timrev=2 if time-reversal symmetry is used, 1 otherwise
OUTPUT
qibz(3,nqibz)=coordinates of q points by which k points differ
SOURCE
2035 subroutine findq(nkbz,kbz,nsym,symrec,symafm,gprimd,nqibz,qibz,timrev) 2036 2037 !Arguments ------------------------------------ 2038 !scalars 2039 integer,intent(in) :: nkbz,nqibz,nsym,timrev 2040 !arrays 2041 integer,intent(in) :: symrec(3,3,nsym),symafm(nsym) 2042 real(dp),intent(in) :: gprimd(3,3),kbz(3,nkbz) 2043 real(dp),intent(inout) :: qibz(3,nqibz) !vz_i num m_numeric_tools 2044 2045 !Local variables ------------------------------ 2046 !scalars 2047 integer :: ii,ik,iq,iqp,isym,itim 2048 real(dp) :: shift1,qred 2049 logical :: found 2050 character(len=500) :: msg 2051 !arrays 2052 integer :: g0(3) 2053 real(dp) :: gmet(3,3),qposs(3),qrot(3) 2054 2055 !************************************************************************ 2056 2057 ! Compute reciprocal space metrics 2058 do ii=1,3 2059 gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+& 2060 & gprimd(2,ii)*gprimd(2,:)+& 2061 & gprimd(3,ii)*gprimd(3,:) 2062 end do 2063 ! 2064 ! === Loop over k-points in BZ, form k-k1 and translate in first BZ === 2065 ! iq is the no. of q-points found, zero at the beginning 2066 iq=0 2067 do ik=1,nkbz 2068 qposs(:)=kbz(:,ik)-kbz(:,1) 2069 ! === Check whether this q (or its equivalent) has already been found === 2070 ! * Use spatial inversion instead of time reversal whenever possible. 2071 found=.FALSE. 2072 do iqp=1,iq 2073 do itim=1,timrev 2074 do isym=1,nsym 2075 if (symafm(isym)==-1) CYCLE 2076 qrot = (3-2*itim) * MATMUL(symrec(:,:,isym),qibz(:,iqp)) 2077 if (isamek(qrot,qposs,g0)) found=.TRUE. 2078 end do 2079 end do 2080 end do 2081 if (.not.found) then 2082 iq=iq+1 2083 if (iq>nqibz) then 2084 ABI_BUG(sjoin('iq > nqibz= ',itoa(nqibz))) 2085 end if 2086 qibz(:,iq)=qposs(:) 2087 end if 2088 end do 2089 2090 if (iq/=nqibz) then 2091 write(msg,'(2(a,i5))')' iq= ',iq,'/= nqibz= ',nqibz 2092 ABI_BUG(msg) 2093 end if 2094 ! 2095 ! * Translate q-points to 1st BZ in the interval [-1/2,1/2[ 2096 do iq=1,nqibz 2097 do ii=1,3 2098 call wrap2_pmhalf(qibz(ii,iq),qred,shift1) 2099 qibz(ii,iq)=qred 2100 end do 2101 end do 2102 2103 end subroutine findq
m_bz_mesh/findqg0 [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
findqg0
FUNCTION
Identify q + g0 = k - kp
INPUTS
kmkp(3)= k - kp input vector nqbz=number of q points in the BZ qbz(3,nqbz)=coordinates of q-points in the BZ mG0(3)= For each reduced direction gives the maximum G0 component to account for umklapp processes
OUTPUT
iq=index of q in qbz array. g0(3)=reciprocal space vector, to be used in igfft
SOURCE
2127 subroutine findqg0(iq, g0, kmkp, nqbz, qbz, mG0) 2128 2129 !Arguments ------------------------------------ 2130 !scalars 2131 integer,intent(in) :: nqbz 2132 integer,intent(out) :: iq 2133 !arrays 2134 integer,intent(in) :: mG0(3) 2135 integer,intent(out) :: g0(3) 2136 real(dp),intent(in) :: kmkp(3),qbz(3,nqbz) 2137 2138 !Local variables------------------------------- 2139 !FIXME if I use 1.0d-4 the jobs crash, should understand why 2140 !scalars 2141 integer :: ig,iqbz,jg01,jg02,jg03 2142 real(dp) :: tolq0=1.0D-3 2143 character(len=500) :: msg 2144 !arrays 2145 real(dp) :: glist1(2*ABS(mG0(1))+1),glist2(2*ABS(mG0(2))+1),glist3(2*ABS(mG0(3))+1) 2146 real(dp) :: qpg0(3),rg(3) 2147 2148 ! ************************************************************************* 2149 2150 iq = 0 2151 2152 if (ALL(ABS(kmkp) < EPSILON(one))) then 2153 ! Find q close to 0 2154 do iqbz=1,nqbz 2155 if (ALL(ABS(qbz(:,iqbz)) < tolq0)) then 2156 iq = iqbz 2157 end if 2158 end do 2159 2160 ABI_CHECK(iq /= 0, 'Wrong list of q-points: q=0 not present.') 2161 g0(:) = 0; RETURN 2162 2163 else 2164 ! q is not zero, find q such as k-kp=q+G0. 2165 2166 ! Try with G0 = 0 first. 2167 !do iqbz=1,nqbz 2168 ! if (ALL(ABS(qbz(:,iqbz)-kmkp)<TOL_KDIFF)) then 2169 ! iq=iqbz 2170 ! g0(:)=0; RETURN 2171 ! end if 2172 !end do 2173 2174 ! Init G0 lists to accelerate search below (small |G0| first) 2175 glist1(1) = 0 2176 ig = 2 2177 do jg01=1,mG0(1) 2178 glist1(ig) = jg01 2179 glist1(ig+1) = -jg01 2180 ig = ig + 2 2181 end do 2182 2183 glist2(1) = 0 2184 ig = 2 2185 do jg02=1,mG0(2) 2186 glist2(ig) = jg02 2187 glist2(ig+1) = -jg02 2188 ig = ig + 2 2189 end do 2190 2191 glist3(1) = 0 2192 ig = 2 2193 do jg03=1,mG0(3) 2194 glist3(ig) = jg03 2195 glist3(ig+1) = -jg03 2196 ig = ig + 2 2197 end do 2198 2199 ! Search algorithm. 2200 g1loop: do jg01=1,2*mG0(1)+1 2201 rg(1) = glist1(jg01) 2202 do jg02=1,2*mG0(2)+1 2203 rg(2) = glist2(jg02) 2204 do jg03=1,2*mG0(3)+1 2205 rg(3) = glist3(jg03) 2206 2207 ! Form q+G0 and check if it is the one. 2208 do iqbz=1,nqbz 2209 qpg0= qbz(:,iqbz) + rg 2210 if (ALL(ABS(qpg0-kmkp) < TOL_KDIFF)) then 2211 iq = iqbz 2212 g0 = NINT(rg) 2213 EXIT g1loop 2214 end if 2215 end do 2216 2217 end do 2218 end do 2219 end do g1loop 2220 2221 if (iq == 0) then 2222 write(msg,'(a, 3f9.5)')' q = k-kp+G0 not found. kmkp:',kmkp 2223 ABI_ERROR(msg) 2224 end if 2225 end if 2226 2227 end subroutine findqg0
m_bz_mesh/get_BZ_diff [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
get_BZ_diff
FUNCTION
Given two points k1 and k2 where k1 belongs to the BZ, check if the difference k1-k2 still belongs to the BZ reporting useful quantities
INPUTS
Kmesh<kmesh_t>=datatype gathering information on the k-mesh k1(3)=the first k-points (supposed to be in the BZ) k2(3)=the second point
OUTPUT
idiff_bz=the idex of k1-k2 in the BZ G0(3)=the umklapp G0 vector required to bring k1-k2 back to the BZ nfound= the number of points in the BZ that are equal to k1-k2 (should be 1 if everything is OK)
SOURCE
967 subroutine get_BZ_diff(Kmesh,k1,k2,idiff_bz,g0,nfound) 968 969 !Arguments ------------------------------------ 970 !scalars 971 class(kmesh_t),intent(in) :: Kmesh 972 integer,intent(out) :: idiff_bz,nfound 973 !arrays 974 integer,intent(out) :: g0(3) 975 real(dp),intent(in) :: k1(3),k2(3) 976 977 !Local variables------------------------------- 978 !scalars 979 integer :: ikp 980 character(len=500) :: msg 981 !arrays 982 integer :: umklp(3) 983 real(dp) :: kdiff(3),ktrial(3) 984 985 ! ********************************************************************* 986 987 if (.not.has_BZ_item(Kmesh,k1,ikp,umklp)) then 988 write(msg,'(a,3f12.6)')' first point must be in BZ: ',k1 989 ABI_ERROR(msg) 990 end if 991 992 kdiff = k1-k2 993 nfound = 0 994 idiff_bz= 0 995 ! 996 ! === Find p such k1-k2=p+g0 where p in the BZ === 997 do ikp=1,Kmesh%nbz 998 ktrial=Kmesh%bz(:,ikp) 999 if (isamek(kdiff,ktrial,umklp)) then 1000 idiff_bz=ikp 1001 g0=umklp 1002 nfound=nfound+1 1003 end if 1004 end do 1005 ! 1006 ! === Check if p has not found of found more than once === 1007 ! * For extremely dense meshes, tol1q in defs_basis might be too large! 1008 if (nfound/=1) then 1009 if (nfound==0) then 1010 ABI_WARNING(" k1-k2-G0 not found in BZ") 1011 else 1012 ABI_WARNING(sjoin(' Multiple k1-k2-G0 found in BZ, nfound= ', itoa(nfound))) 1013 end if 1014 write(msg,'(4a,3(a,3f12.6,a))')& 1015 & ' k1 = ',k1 ,ch10,& 1016 & ' k2 = ',k2 ,ch10,& 1017 & ' k1-k2 = ',kdiff,ch10 1018 ABI_WARNING(msg) 1019 end if 1020 1021 end subroutine get_BZ_diff
m_bz_mesh/get_bz_item [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
get_bz_item
FUNCTION
Given the index of a point in the full BZ, this routine returns the index of the symmetric image in the IBZ, the index of the symmetry operation symrec needed, whether time-reversal has to be used. Optionally the non-symmorphic phase and the umklapp vector is returned.
INPUTS
ikbz=The index of the required point in the BZ Kmesh<kmesh_t>=Datatype gathering information on the k point sampling.
OUTPUT
kbz(3)=The k-point in the first BZ in reduced coordinated. isym=Index of the symrec symmetry required to rotate ik_ibz onto ik_bz. itim=2 is time-reversal has to be used, 1 otherwise ik_ibz=The index of the corresponding symmetric point in the IBZ. [ph_mkbzt]=The phase factor for non-symmorphic operations e^{-i 2 \pi k_IBZ \cdot R{^-1}t}=e{-i 2\pi k_BZ cdot t} [umklp(3)]=The umklapp G0 vector such as kbz + G0 = (IS) k_ibz, where kbz is in the BZ. [isirred]=.TRUE. if the k-point belongs to IBZ.
SOURCE
864 subroutine get_bz_item(Kmesh,ik_bz,kbz,ik_ibz,isym,itim,ph_mkbzt,umklp,isirred) 865 866 !Arguments ------------------------------------ 867 !scalars 868 class(kmesh_t),intent(in) :: Kmesh 869 integer,intent(in) :: ik_bz 870 integer,intent(out) :: ik_ibz,isym,itim 871 complex(dpc),optional,intent(out) :: ph_mkbzt 872 logical,optional,intent(out) :: isirred 873 !arrays 874 integer,optional,intent(out) :: umklp(3) 875 real(dp),intent(out) :: kbz(3) 876 877 !Local variables------------------------------- 878 !scalars 879 character(len=500) :: msg 880 881 ! ********************************************************************* 882 883 if (ik_bz>Kmesh%nbz.or.ik_bz<=0) then 884 write(msg,'(a,2i3)')' Wrong value for ik_bz: ',ik_bz,Kmesh%nbz 885 ABI_BUG(msg) 886 end if 887 888 kbz = Kmesh%bz(:,ik_bz) 889 ik_ibz = Kmesh%tab(ik_bz) 890 isym = Kmesh%tabo(ik_bz) 891 itim = (3-Kmesh%tabi(ik_bz))/2 892 893 if (PRESENT(ph_mkbzt)) ph_mkbzt=Kmesh%tabp(ik_bz) 894 if (PRESENT(umklp)) umklp =Kmesh%umklp(:,ik_bz) 895 ! Be careful here as we assume a particular ordering of symmetries. 896 if (PRESENT(isirred)) isirred = (isym==1.and.itim==1.and.ALL(Kmesh%umklp(:,ik_bz)==(/0,0,0/))) 897 898 end subroutine get_bz_item
m_bz_mesh/get_IBZ_item [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
get_IBZ_item
FUNCTION
Report useful information on a k-point in the IBZ starting from its sequential index in %ibz.
INPUTS
ik_ibz=The index of the required point in the IBZ Kmesh<kmesh_t>=datatype gathering information on the k point sampling.
OUTPUT
kibz(3)=the k-point in reduced coordinated wtk=the weight
TODO
Add mapping ibz2bz, ibz2star
SOURCE
923 subroutine get_IBZ_item(Kmesh,ik_ibz,kibz,wtk) 924 925 !Arguments ------------------------------------ 926 !scalars 927 class(kmesh_t),intent(in) :: Kmesh 928 integer,intent(in) :: ik_ibz 929 real(dp),intent(out) :: wtk 930 !arrays 931 real(dp),intent(out) :: kibz(3) 932 933 ! ********************************************************************* 934 935 if (ik_ibz>Kmesh%nibz.or.ik_ibz<=0) then 936 ABI_BUG(sjoin('wrong value for ik_ibz: ',itoa(ik_ibz))) 937 end if 938 939 kibz=Kmesh%ibz(:,ik_ibz) 940 wtk =Kmesh%wt(ik_ibz) 941 942 end subroutine get_IBZ_item
m_bz_mesh/get_ng0sh [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
get_ng0sh
FUNCTION
Given two lists of k-points, kbz1 and kbz2, calculate any possible difference k1-k2. For each difference, find the umklapp g0 vector and the point k3 in the array kfold such as k1-k2 = k3 + G0. The optimal value of G0 shells is returned, namely the smallest box around Gamma which suffices to treat all possible umklapp processes. The search algorithm uses bisection to process to scale in nk1*nk2*log(nkfold)
INPUTS
nk1, nk2=Number of points in the arrays kbz1, kbz2. kbz1(3,nk1)=Reduced coordinates of the first set of points. kbz2(3,nk2)=Reduced coordinates of the second set of points. nkfold=Number of points in the array kfold. kfold(3,nkfol)=Reduced coordinated of the points in the BZ. tolq0=Tolerance below which a q-point is treated as zero.
OUTPUT
opt_ng0(3)=Minimal reduced components of the G0 vectors to account for umklapps.
SOURCE
1600 subroutine get_ng0sh(nk1,kbz1,nk2,kbz2,nkfold,kfold,tolq0,opt_ng0) 1601 1602 !Arguments ------------------------------------ 1603 !scalars 1604 integer,intent(in) :: nk1,nk2,nkfold 1605 real(dp),intent(in) :: tolq0 1606 !arrays 1607 integer,intent(out) :: opt_ng0(3) 1608 real(dp),intent(in) :: kbz1(3,nk1),kbz2(3,nk2),kfold(3,nkfold) 1609 1610 !Local variables------------------------------- 1611 !scalars 1612 integer :: i1,i2,ikf,ind,factor 1613 real(dp) :: normdiff,tempnorm,smallestlen 1614 logical :: found 1615 character(len=500) :: msg 1616 !arrays 1617 integer :: roundk(3),kbigdiff(3),kbigfold(3,nkfold),iperm(nkfold) 1618 real(dp) :: k1mk2(3),ksmalldiff(3),norm(nkfold),ksmallfold(3,nkfold) 1619 1620 !************************************************************************ 1621 1622 ! Compute smallest length of one component 1623 ! To get a sufficiently large factor to order vectors 1624 smallestlen = one 1625 1626 ! Compute integer part and fractional part, [0,1[ of kfold 1627 do ikf = 1,nkfold 1628 kbigfold(:,ikf) = FLOOR(kfold(:,ikf)+tol7) 1629 ksmallfold(:,ikf) = kfold(:,ikf)-kbigfold(:,ikf) 1630 1631 if (ABS(ksmallfold(1,ikf)) > tol7) & 1632 & smallestlen = MIN(smallestlen, ABS(ksmallfold(1,ikf))) 1633 if (ABS(ksmallfold(2,ikf)) > tol7) & 1634 & smallestlen = MIN(smallestlen, ABS(ksmallfold(2,ikf))) 1635 if (ABS(ksmallfold(3,ikf)) > tol7) & 1636 & smallestlen = MIN(smallestlen, ABS(ksmallfold(3,ikf))) 1637 1638 end do 1639 1640 ! WARNING ! This could not be sufficient if tested k1 - k2 has lower 1641 ! components than smallestlen. The factor 10 is giving us a security margin. 1642 factor = 10*(int(one/smallestlen)+1) 1643 1644 ! Loop again over kfold vectors to give each term its norm 1645 do ikf=1, nkfold 1646 iperm(ikf) = ikf 1647 1648 ! Computing a sort of norm with order of components (used for ordering) 1649 call getkptnorm_bycomponent(ksmallfold(:,ikf),factor,tempnorm) 1650 1651 norm(ikf) = tempnorm 1652 end do 1653 1654 ! Sorting list of kfold vectors 1655 call sort_dp(nkfold,norm,iperm,tol14) 1656 1657 ! Loop over all k1 - k2 1658 opt_ng0(:)=0 1659 do i2=1,nk2 1660 ! This is used in case of screening calculation. 1661 ! If q is small treat it as zero. In this case, indeed, 1662 ! we use q=0 to calculate the oscillator matrix elements. 1663 if (is_zero(kbz2(:,i2),tolq0)) CYCLE 1664 do i1=1,nk1 1665 found=.FALSE. 1666 1667 ! Separating in integer part and fractionary part 1668 k1mk2(:) = kbz1(:,i1)-kbz2(:,i2) 1669 ! Adding small tol to prevent 1 being in fractionary part 1670 kbigdiff(:) = FLOOR(k1mk2(:)+tol7) 1671 ksmalldiff(:) = k1mk2(:)-kbigdiff(:) 1672 1673 call getkptnorm_bycomponent(ksmalldiff(:),factor,normdiff) 1674 1675 ! Try to find the right smallkfold, corresponding to ksmalldiff 1676 ind = bisect(norm,normdiff) 1677 1678 if (ind > 0) then 1679 if(ABS(norm(ind) - normdiff) < TOL_KDIFF) then 1680 found = .TRUE. 1681 end if 1682 end if 1683 if(ind < nkfold) then 1684 if(ABS(norm(ind+1) - normdiff) < TOL_KDIFF) then 1685 found = .TRUE. 1686 ind = ind + 1 1687 end if 1688 end if 1689 1690 if (.not. found) then 1691 write(msg,'(a,2(2a,i4,3es16.8),a)')& 1692 & 'Not able to found umklapp G0 vector such as k1-k2 = kf+G0',ch10,& 1693 & 'point1 = ',i1,kbz1(:,i1),ch10,& 1694 & 'point2 = ',i2,kbz2(:,i2),ch10 1695 ABI_ERROR(msg) 1696 else 1697 ! We have found one k, extracting the max g0 1698 roundk(:) = ABS(kbigdiff(:) - kbigfold(:,iperm(ind))) 1699 opt_ng0(1) = MAX(opt_ng0(1),roundk(1)) 1700 opt_ng0(2) = MAX(opt_ng0(2),roundk(2)) 1701 opt_ng0(3) = MAX(opt_ng0(3),roundk(3)) 1702 end if 1703 end do 1704 end do 1705 1706 end subroutine get_ng0sh
m_bz_mesh/getkptnorm_bycomponent [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
getkptnorm_bycomponent
FUNCTION
get the norm of one vector, by order of its components
INPUTS
vect(3) = vector which should be evaluated
factor = term which multiplies the components
WARNING ! should be large enough to get unique order
OUTPUT
norm = value of the norm
SOURCE
1728 subroutine getkptnorm_bycomponent(vect,factor,norm) 1729 1730 !Arguments ------------------------------------ 1731 !scalars 1732 integer,intent(in) :: factor 1733 real(dp),intent(out):: norm 1734 !arrays 1735 real(dp),intent(in) :: vect(3) 1736 !Local variables------------------------------- 1737 !scalars 1738 character(len=500) :: msg 1739 1740 ! ************************************************************************* 1741 1742 ! Checking the factor is large enough 1743 !(skipping zero components, since in this case the product will be 0) 1744 if(ANY(vect(:)*factor < 1.0 .and. vect(:) > tol7)) then 1745 write(msg,'(a,a,a,a,a,a,a,a)') ' Not able to give unique norm to order vectors',ch10,& 1746 & 'This is likely related to a truncation error for a k-point in the input file',ch10,& 1747 & 'Always prefer fractional numbers in the input file instead of truncated ones',ch10,& 1748 & '(e.g. 1/6 instead of 0.166666667)',ch10 1749 ABI_ERROR(msg) 1750 end if 1751 1752 norm = (vect(1)*factor+vect(2))*factor+vect(3) 1753 1754 end subroutine getkptnorm_bycomponent
m_bz_mesh/has_BZ_item [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
has_BZ_item
FUNCTION
check if item belongs to the BZ within a reciprocal lattice vector and return the index number and the reciprocal vector g0.
INPUTS
Kmesh<kmesh_t>=datatype gathering information on the k-mesh item(3)=the k-point to be checked
OUTPUT
.TRUE. if item is the BZ within a RL vector ikbz=Index of the k-point in the Kmesh%bz array g0(3)=Umklapp vector. FIXME Switch to routine version. Due to side-effects the present implementation might be source of bugs in logical statements
SOURCE
1120 logical function has_BZ_item(Kmesh, item, ikbz, g0) 1121 1122 !Arguments ------------------------------------ 1123 !scalars 1124 class(kmesh_t),intent(in) :: Kmesh 1125 integer,intent(out) :: ikbz 1126 !arrays 1127 integer,intent(out) :: g0(3) 1128 real(dp),intent(in) :: item(3) 1129 1130 !Local variables------------------------------- 1131 !scalars 1132 integer :: ik_bz,yetfound 1133 !arrays 1134 integer :: g0_tmp(3) 1135 1136 ! ************************************************************************* 1137 1138 has_BZ_item=.FALSE.; ikbz=0; g0=0; yetfound=0 1139 do ik_bz=1,Kmesh%nbz 1140 if (isamek(item,Kmesh%bz(:,ik_bz),g0_tmp)) then 1141 has_BZ_item=.TRUE. 1142 ikbz=ik_bz 1143 g0 = g0_tmp 1144 yetfound=yetfound+1 1145 !EXIT 1146 end if 1147 end do 1148 1149 if (yetfound/=0.and.yetfound/=1) then 1150 ABI_ERROR('Multiple k-points found') 1151 end if 1152 1153 end function has_BZ_item
m_bz_mesh/has_IBZ_item [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
has_IBZ_item
FUNCTION
Check if item belongs to the IBZ within a reciprocal lattice vector
INPUTS
Kmesh<kmesh_t>=Datatype gathering information on the mesh in the BZ. item(3)=the k-point to be checked
OUTPUT
Return .TRUE. if item is the IBZ within a RL vector ikibz=The index of the k-point in the IBZ. g0(3)=The reciprocal lattice vector.
SOURCE
1176 logical function has_IBZ_item(Kmesh,item,ikibz,g0) 1177 1178 !Arguments ------------------------------------ 1179 !scalars 1180 class(kmesh_t),intent(in) :: Kmesh 1181 integer,intent(out) :: ikibz 1182 !arrays 1183 integer,intent(out) :: g0(3) 1184 real(dp),intent(in) :: item(3) 1185 1186 !Local variables------------------------------- 1187 !scalars 1188 integer :: ik_ibz,yetfound 1189 !character(len=500) :: msg 1190 !arrays 1191 integer :: g0_tmp(3) 1192 1193 ! ************************************************************************* 1194 1195 has_IBZ_item=.FALSE.; ikibz=0; g0=0; yetfound=0 1196 do ik_ibz=1,Kmesh%nibz 1197 if (isamek(item,Kmesh%ibz(:,ik_ibz),g0_tmp)) then 1198 has_IBZ_item=.TRUE. 1199 ikibz=ik_ibz 1200 g0 = g0_tmp 1201 yetfound=yetfound+1 1202 !EXIT 1203 end if 1204 end do 1205 1206 if (yetfound/=0.and.yetfound/=1) then 1207 ABI_BUG("multiple k-points found") 1208 end if 1209 1210 end function has_IBZ_item
m_bz_mesh/identk [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
identk
FUNCTION
Identify k-points in the whole BZ starting from the IBZ. Generate also symmetry tables relating the BZ to the IBZ.
INPUTS
kibz(3,nkibz)=Coordinates of k-points in the IBZ. nkibz=Number of k points in IBZ. nkbzmx=Maximum number of k points in BZ. nsym=Number of symmetry operations. timrev=2 if time reversal symmetry can be used; 1 otherwise. symrec(3,3,nsym)=Symmetry operation matrices in reciprocal space. symafm(nsym)=(anti)ferromagnetic part of symmetry operations. [ref_bz(:,:)]= Reference set of points in the full Brillouin zone.
OUTPUT
kbz(3,nkbzmx)= k-points in whole BZ
ktab(nkbzmx)= table giving for each k-point in the BZ (array kbz),
the corresponding irreducible point in the array (kibz)
k_BZ= (IS) kIBZ where S is one of the symrec operations and I is the inversion or the identity
where k_BZ = (IS) k_IBZ and S = \transpose R^{-1}
ktabi(nkbzmx)= for each k-point in the BZ defines whether inversion has to be
considered in the relation k_BZ=(IS) k_IBZ (1 => only S; -1 => -S)
ktabo(nkbzmx)= the symmetry operation S that takes k_IBZ to each k_BZ
nkbz= no. of k-points in the whole BZ
wtk(nkibz)= weight for each k-point in IBZ for symmetric quantities:
no. of distinct ks in whole BZ/(timrev*nsym)
NOTES
The logic of the routine relies on the assumption that kibz really represent an irreducible set. If symmetrical points are present in the input list, indeed, some the output weights will turn out to be zero. An initial check is done at the beginning of the routine to trap this possible error.
SOURCE
1416 subroutine identk(kibz,nkibz,nkbzmx,nsym,timrev,symrec,symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk,ref_bz) 1417 1418 !Arguments ------------------------------------ 1419 !scalars 1420 integer,intent(in) :: nkbzmx,nkibz,nsym,timrev 1421 integer,intent(out) :: nkbz 1422 !arrays 1423 integer,intent(in) :: symafm(nsym),symrec(3,3,nsym) 1424 integer,intent(out) :: ktab(nkbzmx),ktabi(nkbzmx),ktabo(nkbzmx) 1425 real(dp),intent(in) :: kibz(3,nkibz) 1426 real(dp),intent(out) :: kbz(3,nkbzmx),wtk(nkibz) 1427 real(dp),optional,intent(in) :: ref_bz(:,:) 1428 1429 !Local variables ------------------------------ 1430 !scalars 1431 integer :: ik1,ik2,ikbz,ikibz,iold,isym,itim 1432 integer :: ikref,nkref,isym_swp,itim_swp,ikibz_swp 1433 logical :: is_irred_set 1434 logical :: found,ltest 1435 character(len=500) :: msg 1436 !arrays 1437 integer :: g0(3) 1438 real(dp) :: knew(3),k1(3),k2(3) 1439 real(dp) :: kref(3),kbz_swp(3) 1440 1441 ! ************************************************************************* 1442 1443 DBG_ENTER("COLL") 1444 ! 1445 ! === Check whether kibz really forms an irreducible set === 1446 is_irred_set=.TRUE. 1447 do ik1=1,nkibz-1 1448 k1=kibz(:,ik1) 1449 do ik2=ik1+1,nkibz 1450 k2=kibz(:,ik2) 1451 1452 do itim=1,timrev 1453 do isym=1,nsym 1454 if (symafm(isym)==-1) CYCLE 1455 knew = (3-2*itim) * MATMUL(symrec(:,:,isym),k2) 1456 if (isamek(k1,knew,g0)) then 1457 is_irred_set=.FALSE. 1458 write(msg,'(2(a,3f8.4),2(a,i2))')& 1459 & ' k1 = ',k1,' is symmetrical of k2 = ',k2,' through sym = ',isym,' itim = ',itim 1460 ABI_WARNING(msg) 1461 end if 1462 end do 1463 end do 1464 1465 end do 1466 end do 1467 1468 !call klist_isirred(nkibz,kibz,Cryst,nimg) 1469 1470 if (.not.is_irred_set) then 1471 ABI_WARNING("Input array kibz does not constitute an irreducible set.") 1472 end if 1473 ! 1474 ! === Loop over k-points in IBZ === 1475 ! * Start with zero no. of k-points found. 1476 nkbz=0 1477 do ikibz=1,nkibz 1478 wtk(ikibz) = zero 1479 1480 ! === Loop over time-reversal I and symmetry operations S === 1481 ! * Use spatial inversion instead of time reversal whenever possible. 1482 do itim=1,timrev 1483 do isym=1,nsym 1484 if (symafm(isym)==-1) CYCLE 1485 ! 1486 ! * Form IS k 1487 knew=(3-2*itim)*MATMUL(symrec(:,:,isym),kibz(:,ikibz)) 1488 ! 1489 ! * Check whether it has already been found (to within a RL vector). 1490 iold=0 1491 do ikbz=1,nkbz 1492 if (isamek(knew,kbz(:,ikbz),g0)) then 1493 iold=iold+1 1494 exit 1495 end if 1496 end do 1497 ! 1498 ! * If not yet found add to kbz and increase the weight. 1499 if (iold==0) then 1500 nkbz=nkbz+1 1501 wtk(ikibz)=wtk(ikibz)+one 1502 if (nkbz>nkbzmx) then 1503 ABI_BUG(sjoin('nkbzmx too small, nkbzmx = ',itoa(nkbzmx),', increase nkbzmx !')) 1504 end if 1505 kbz(:,nkbz) = knew(:) 1506 ktab (nkbz) = ikibz 1507 ktabo(nkbz) = isym 1508 ktabi(nkbz) = 3-2*itim 1509 end if 1510 ! 1511 end do 1512 end do 1513 1514 end do !ikibz 1515 1516 if (PRESENT(ref_bz)) then 1517 call wrtout(std_out," Pruning the k-points not in ref_bz then reordering tables","COLL") 1518 1519 nkref = SIZE(ref_bz,DIM=2) 1520 ltest = (nkref>=nkbz.and.nkref<=nkbzmx) 1521 if (.not.ltest) then 1522 write(msg,'(3(a,i0))')" Wrong value for nkref: nkref= ",nkref," nkbz= ",nkbz," nkbzmx =",nkbzmx 1523 ABI_WARNING(msg) 1524 end if 1525 1526 do ikref=1,nkref 1527 kref = ref_bz(:,ikref) 1528 found=.FALSE. 1529 1530 do ikbz=1,nkbz ! Loop on the set of BZ points found above. 1531 if (isequalk(kref,kbz(:,ikbz))) then ! Swap indices. 1532 kbz_swp = kbz(:,ikref) 1533 ikibz_swp = ktab (ikref) 1534 isym_swp = ktabo(ikref) 1535 itim_swp = ktabi(ikref) 1536 1537 kbz(:,ikref) = kref 1538 ktab (ikref) = ktab (ikbz) 1539 ktabo(ikref) = ktabo(ikbz) 1540 ktabi(ikref) = ktabi(ikbz) 1541 1542 kbz(:,ikbz) = kbz_swp 1543 ktab (ikbz) = ikibz_swp 1544 ktabo(ikbz) = isym_swp 1545 ktabi(ikbz) = itim_swp 1546 1547 found=.TRUE.; EXIT 1548 end if 1549 end do 1550 1551 if (.not.found) then 1552 write(msg,'(a,3es16.8)')" One of the k-point in ref_bz is not a symmetrical image of the IBZ: ",kref 1553 ABI_ERROR(msg) 1554 end if 1555 end do 1556 ! 1557 ! Change nkbz to nkref, then get the new weights. 1558 nkbz=nkref; wtk=zero 1559 do ikref=1,nkref 1560 ikibz = ktab(ikref) 1561 wtk(ikibz) = wtk(ikibz) + 1 1562 end do 1563 end if ! PRESENT(ref_bz) 1564 ! 1565 ! * Weights are normalized to 1. 1566 wtk = wtk/SUM(wtk) 1567 1568 DBG_EXIT("COLL") 1569 1570 end subroutine identk
m_bz_mesh/isamek [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
isamek
FUNCTION
Test two k-points for equality. Return .TRUE. is they are equal within a reciprocal lattice vector G0.
INPUTS
k1(3),k2(3)=The two k points to be compared.
OUTPUT
Return .TRUE. if they are the same within a RL vector,
.FALSE. if they are different.
G0(3)=if .TRUE. G0(3) is the reciprocal lattice vector such as k1=k2+G0
SOURCE
1044 logical function isamek(k1, k2, g0) 1045 1046 !Arguments ------------------------------------ 1047 !arrays 1048 integer,intent(out) :: g0(3) 1049 real(dp),intent(in) :: k1(3),k2(3) 1050 1051 ! ************************************************************************* 1052 1053 isamek = isinteger(k1-k2,TOL_KDIFF) 1054 1055 if (isamek) then 1056 g0=NINT(k1-k2) 1057 else 1058 g0=HUGE(1) 1059 end if 1060 1061 end function isamek
m_bz_mesh/isequalk [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
is_equalk
FUNCTION
Return .TRUE. if two points are equal within a reciprocal lattice vector.
INPUTS
q1(3),q2(3)=The two points to be compared for equivalence.
OUTPUT
SOURCE
1080 logical function isequalk(q1, q2) 1081 1082 !Arguments ------------------------------------ 1083 real(dp),intent(in) :: q1(3),q2(3) 1084 1085 !Local variables------------------------------- 1086 !arrays 1087 integer :: g0(3) 1088 ! ************************************************************************* 1089 1090 isequalk = isamek(q1,q2,g0) 1091 1092 end function isequalk
m_bz_mesh/kmesh_free [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
kmesh_free
FUNCTION
Deallocate all dynamics entities present in a kmesh_t structure.
INPUTS
Kmesh<kmesh_t>=The datatype to be freed.
SOURCE
560 subroutine kmesh_free(Kmesh) 561 562 !Arguments ------------------------------------ 563 class(kmesh_t),intent(inout) :: Kmesh 564 565 ! ********************************************************************* 566 567 !@kmesh_t 568 !integer 569 ABI_SFREE(Kmesh%rottb) 570 ABI_SFREE(Kmesh%rottbm1) 571 ABI_SFREE(Kmesh%tab) 572 ABI_SFREE(Kmesh%tabi) 573 ABI_SFREE(Kmesh%tabo) 574 ABI_SFREE(Kmesh%umklp) 575 576 !real 577 ABI_SFREE(Kmesh%ibz) 578 ABI_SFREE(Kmesh%bz) 579 ABI_SFREE(Kmesh%shift) 580 ABI_SFREE(Kmesh%wt) 581 582 !complex 583 ABI_SFREE(Kmesh%tabp) 584 585 end subroutine kmesh_free
m_bz_mesh/kmesh_init [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
kmesh_init
FUNCTION
Initialize and construct a kmesh_t datatype gathering information on the mesh in the Brilloin zone.
INPUTS
nkibz=Number of irreducible k-points. kibz(3,nkibz)=Irreducible k-points in reduced coordinates. Cryst<crystal_t> = Info on unit cell and its symmetries %nsym=number of symmetry operations %symrec(3,3,nsym)=symmetry operations in reciprocal space %tnons(3,nsym)=fractional translations kptopt=option for the generation of k points (see input variable description) [wrap_1zone]=If .TRUE., the points are wrapped in in the first BZ. Defaults to .FALSE. to preserve GW implementation. [ref_bz(:,:)]= Reference set of points in the full Brillouin zone used to prune k-points.
OUTPUT
Kmesh<kmesh_t>=Datatype gathering information on the k point sampling.
SOURCE
402 subroutine kmesh_init(Kmesh,Cryst,nkibz,kibz,kptopt,wrap_1zone,ref_bz,break_symmetry) 403 404 !Arguments ------------------------------------ 405 !scalars 406 class(kmesh_t),intent(inout) :: Kmesh 407 integer,intent(in) :: nkibz,kptopt 408 logical,optional,intent(in) :: wrap_1zone,break_symmetry 409 type(crystal_t),intent(in) :: Cryst 410 !arrays 411 real(dp),intent(in) :: kibz(3,nkibz) 412 real(dp),optional,intent(in) :: ref_bz(:,:) 413 414 !Local variables------------------------------- 415 !scalars 416 integer :: ik_bz,ik_ibz,isym,nkbz,nkbzX,nsym,timrev,itim 417 real(dp) :: shift1,shift2,shift3 418 logical :: ltest,do_wrap,do_hack 419 !arrays 420 integer,allocatable :: ktab(:),ktabi(:),ktabo(:) 421 real(dp) :: rm1t(3),kbz_wrap(3) 422 real(dp),allocatable :: kbz(:,:),wtk(:) 423 424 ! ************************************************************************* 425 426 DBG_ENTER("COLL") 427 428 !@kmesh_t 429 ! === Initial tests on input arguments === 430 ltest=(Cryst%timrev==1.or.Cryst%timrev==2) 431 ABI_CHECK(ltest, sjoin('Wrong value for timrev= ', itoa(Cryst%timrev))) 432 433 if (ALL(kptopt/=(/1,3/))) then 434 ABI_WARNING(sjoin("Not allowed value for kptopt: ", itoa(kptopt))) 435 end if 436 437 Kmesh%kptopt = kptopt 438 439 nsym = Cryst%nsym 440 timrev = Cryst%timrev 441 442 ! 443 ! === Find BZ from IBZ and fill tables === 444 nkbzX=nkibz*nsym*timrev ! Maximum possible number 445 ABI_MALLOC(kbz,(3,nkbzX)) 446 ABI_MALLOC(wtk,(nkibz)) 447 ABI_MALLOC(ktab,(nkbzX)) 448 ABI_MALLOC(ktabi,(nkbzX)) 449 ABI_MALLOC(ktabo,(nkbzX)) 450 451 if (PRESENT(ref_bz)) then 452 call identk(kibz,nkibz,nkbzX,nsym,timrev,cryst%symrec,cryst%symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk,ref_bz=ref_bz) 453 else 454 call identk(kibz,nkibz,nkbzX,nsym,timrev,cryst%symrec,cryst%symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk) 455 end if 456 457 ! TODO: Force the k-points to be in the first Brillouin zone. 458 ! Now the GW tests seem to be OK, additional tests have to be done though. 459 do_wrap=.FALSE.; if (PRESENT(wrap_1zone)) do_wrap=wrap_1zone 460 !do_wrap=.TRUE. 461 462 if (do_wrap) then ! Wrap the BZ points in the interval ]-1/2,1/2] 463 do ik_bz=1,nkbz 464 call wrap2_pmhalf(kbz(1,ik_bz),kbz_wrap(1),shift1) 465 call wrap2_pmhalf(kbz(2,ik_bz),kbz_wrap(2),shift2) 466 call wrap2_pmhalf(kbz(3,ik_bz),kbz_wrap(3),shift3) 467 kbz(:,ik_bz) = kbz_wrap 468 end do 469 end if 470 ! 471 ! ================================================================ 472 ! ==== Create data structure to store information on k-points ==== 473 ! ================================================================ 474 ! 475 ! * Dimensions. 476 Kmesh%nbz =nkbz ! Number of points in the full BZ 477 Kmesh%nibz =nkibz ! Number of points in the IBZ 478 Kmesh%nsym =nsym ! Number of operations 479 Kmesh%timrev=timrev ! 2 if time-reversal is used, 1 otherwise 480 ! 481 ! * Arrays. 482 Kmesh%gmet = Cryst%gmet 483 Kmesh%gprimd = Cryst%gprimd 484 485 ABI_MALLOC(Kmesh%bz ,(3,nkbz)) 486 Kmesh%bz = kbz(:,1:nkbz ) ! Red. coordinates of points in full BZ. 487 ABI_MALLOC(Kmesh%ibz,(3,nkibz)) 488 Kmesh%ibz = kibz(:,1:nkibz) ! Red. coordinates of points in IBZ. 489 490 ABI_MALLOC(Kmesh%tab ,(nkbz)) 491 Kmesh%tab = ktab (1:nkbz) ! Index of the irred. point in the array IBZ. 492 ABI_MALLOC(Kmesh%tabi,(nkbz)) 493 Kmesh%tabi = ktabi(1:nkbz) !-1 if time reversal must be used to obtain this point, 494 ABI_MALLOC(Kmesh%tabo,(nkbz)) 495 Kmesh%tabo = ktabo(1:nkbz) ! Symm. operation that rotates k_IBZ onto \pm k_BZ 496 ! (depending on tabi) 497 ABI_MALLOC(Kmesh%wt,(nkibz)) 498 Kmesh%wt(:)= wtk(1:nkibz) ! Weight for each k_IBZ 499 500 ABI_MALLOC(Kmesh%rottbm1,(nkbz,timrev,nsym)) 501 ABI_MALLOC(Kmesh%rottb ,(nkbz,timrev,nsym)) 502 503 do_hack = .FALSE. 504 if (PRESENT(ref_bz) .and. PRESENT(break_symmetry)) then 505 do_hack = break_symmetry 506 end if 507 508 if (do_hack) then 509 ABI_WARNING("Hacking the rottb tables!") 510 do ik_bz=1,nkbz 511 Kmesh%rottbm1(ik_bz,:,:) = ik_bz 512 Kmesh%rottb (ik_bz,:,:) = ik_bz 513 end do 514 else 515 call setup_k_rotation(nsym,timrev,cryst%symrec,nkbz,Kmesh%bz,Cryst%gmet,Kmesh%rottb,Kmesh%rottbm1) 516 end if 517 518 ! TODO umklp can be calculated inside setup_k_rotation. 519 ABI_MALLOC(Kmesh%umklp,(3,nkbz)) 520 do ik_bz=1,nkbz 521 ik_ibz= Kmesh%tab (ik_bz) 522 isym = Kmesh%tabo(ik_bz) 523 itim = (3-Kmesh%tabi(ik_bz))/2 524 Kmesh%umklp(:,ik_bz) = NINT( -Kmesh%bz(:,ik_bz) + (3-2*itim)*MATMUL(cryst%symrec(:,:,isym),Kmesh%ibz(:,ik_ibz)) ) 525 end do 526 527 ABI_MALLOC(Kmesh%tabp,(nkbz)) 528 do ik_bz=1,nkbz 529 isym =Kmesh%tabo(ik_bz) 530 ik_ibz=Kmesh%tab (ik_bz) 531 rm1t=MATMUL(TRANSPOSE(cryst%symrec(:,:,isym)),cryst%tnons(:,isym)) 532 Kmesh%tabp(ik_bz)=EXP(-(0.,1.)*two_pi*DOT_PRODUCT(kibz(:,ik_ibz),rm1t)) 533 end do 534 535 ABI_FREE(kbz) 536 ABI_FREE(wtk) 537 ABI_FREE(ktab) 538 ABI_FREE(ktabi) 539 ABI_FREE(ktabo) 540 541 DBG_EXIT("COLL") 542 543 end subroutine kmesh_init
m_bz_mesh/kmesh_print [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
kmesh_print
FUNCTION
Print the content of a kmesh_t datatype
INPUTS
Kmesh<kmesh_t>=the datatype to be printed [header]=optional header [unit]=the unit number for output [prtvol]=verbosity level [mode_paral]=either "COLL" or "PERS"
OUTPUT
Only printing.
SOURCE
609 subroutine kmesh_print(Kmesh,header,unit,prtvol,mode_paral) 610 611 !Arguments ------------------------------------ 612 !scalars 613 class(kmesh_t),intent(in) :: Kmesh 614 integer,optional,intent(in) :: prtvol,unit 615 character(len=4),optional,intent(in) :: mode_paral 616 character(len=*),optional,intent(in) :: header 617 618 !Local variables------------------------------- 619 !scalars 620 integer,parameter :: nmaxk=50 621 integer :: ii,ik,my_unt,my_prtvol 622 character(len=100) :: fmt 623 character(len=4) :: my_mode 624 character(len=500) :: msg 625 626 ! ************************************************************************* 627 628 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 629 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 630 my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 631 632 msg=' ==== Info on the Kmesh% object ==== ' 633 if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== ' 634 call wrtout(my_unt,msg,my_mode) 635 636 write(msg,'(a,i5,3a)')& 637 ' Number of points in the irreducible wedge : ',Kmesh%nibz,ch10,& 638 ' Reduced coordinates and weights : ',ch10 639 call wrtout(my_unt,msg,my_mode) 640 641 write(fmt,*)'(1x,i5,a,2x,3es16.8,3x,f11.5)' 642 do ik=1,Kmesh%nibz ! Add tol8 for portability reasons. 643 write(msg,fmt) ik,') ',(Kmesh%ibz(ii,ik),ii=1,3),Kmesh%wt(ik)+tol8 644 call wrtout(my_unt,msg,my_mode) 645 end do 646 647 SELECT CASE (Kmesh%timrev) 648 CASE (1) 649 write(msg,'(2a,i2,3a,i5,a)')ch10,& 650 ' Together with ',Kmesh%nsym,' symmetry operations (time-reversal symmetry not used) ',ch10,& 651 ' yields ',Kmesh%nbz,' points in the full Brillouin Zone.' 652 653 CASE (2) 654 write(msg,'(2a,i2,3a,i5,a)')ch10,& 655 ' Together with ',Kmesh%nsym,' symmetry operations and time-reversal symmetry ',ch10,& 656 ' yields ',Kmesh%nbz,' points in the full Brillouin Zone.' 657 658 CASE DEFAULT 659 ABI_BUG(sjoin('Wrong value for timrev:', itoa(Kmesh%timrev))) 660 END SELECT 661 662 call wrtout(my_unt,msg,my_mode) 663 664 if (my_prtvol > 0) then 665 write(fmt,*)'(1x,i5,a,2x,3es16.8)' 666 do ik=1,Kmesh%nbz 667 if (my_prtvol==1 .and. ik>nmaxk) then 668 write(msg,'(a)')' prtvol=1, do not print more points.' 669 call wrtout(my_unt,msg,my_mode); EXIT 670 end if 671 write(msg,fmt)ik,') ',(Kmesh%bz(ii,ik),ii=1,3) 672 call wrtout(my_unt,msg,my_mode) 673 end do 674 end if 675 676 ! Additional printing 677 if (my_prtvol >= 10) then 678 write(msg,'(2a)')ch10,& 679 ' Full point -------> Irred point --> through: Symrec Time-Rev (1=No,-1=Yes) G0(1:3) ' 680 call wrtout(my_unt,msg,my_mode) 681 write(fmt,*)'(2x,i5,2x,2(3(f7.4,2x)),i3,2x,i2,3(i3))' 682 do ik=1,Kmesh%nbz 683 write(msg,fmt)ik,Kmesh%bz(:,ik),Kmesh%ibz(:,Kmesh%tab(ik)),Kmesh%tabo(ik),Kmesh%tabi(ik),Kmesh%umklp(:,ik) 684 call wrtout(my_unt,msg,my_mode) 685 end do 686 end if 687 688 write(msg,'(a)')ch10 689 call wrtout(my_unt,msg,my_mode) 690 691 end subroutine kmesh_print
m_bz_mesh/kmesh_t [ Types ]
[ Top ] [ m_bz_mesh ] [ Types ]
NAME
kmesh_t
FUNCTION
The kmesh_t structured datatype contains different information on the grid used to sample the BZ : the k-points in the full Brillouin zone BZ, the irreducible wedge IBZ as well as tables describing the symmetry relationship between the points.
SOURCE
80 type,public :: kmesh_t 81 82 !scalars 83 integer :: nshift = 0 84 85 integer :: nbz = 0 86 ! Number of points in the BZ. 87 88 integer :: nibz = 0 89 ! Number of points in the IBZ. 90 91 integer :: nsym 92 ! Number of symmetry operations. 93 94 integer :: kptopt 95 ! kptopt=option for the generation of k points (see input variable description) 96 ! 97 ! 1 if both time-reversal and point group symmetries are used. 98 ! 2 if only time-reversal symmetry is used. 99 ! 3 do not take into account any symmetry (except the identity). 100 ! 4 if time-reversal is not used (spin-orbit coupling). 101 ! <0 number of segments used to construct the k-path for NSCF calculation. 102 103 integer :: timrev 104 ! 2 if time reversal symmetry can be used, 1 otherwise. 105 106 !arrays 107 integer :: kptrlatt(3,3) = NONE_KPTRLATT 108 ! Coordinates of three vectors in real space, expressed in reduced coordinates. 109 ! They define a super-lattice in real space. The k point lattice is the reciprocal of 110 ! this super-lattice, eventually shifted by shift. 111 ! Not available if the structure is initialized from the points in the IBZ. 112 113 integer,allocatable :: rottb(:,:,:) 114 ! rottb(nbz,timrev,nsym), 115 ! Index of (IS)k in the BZ array where S is a sym operation in reciprocal space, 116 ! I is the identity or the inversion operator (1,2 resp) 117 118 integer,allocatable :: rottbm1(:,:,:) 119 ! rottbm1(nbz,timrev,nsym) 120 ! Index of IS^{-1} k in the BZ array. 121 122 integer,allocatable :: tab(:) 123 ! tab(nbz) 124 ! For each point in the BZ, it gives the index of the symmetric irreducible point in the array ibz. 125 126 integer,allocatable :: tabi(:) 127 ! tabi(nbz) 128 ! For each point in the BZ, tabi tells whether time-reversal has to be 129 ! used to obtain k_BZ starting from the corresponding point in the IBZ (1=>no, -1=>yes) 130 131 integer,allocatable :: tabo(:) 132 ! tabo(nbz) 133 ! For each point in the BZ, it gives the index in the array symrec of the 134 ! symmetry operation in reciprocal space which rotates k_IBZ onto \pm k_BZ (depending on tabi) 135 136 integer,allocatable :: umklp(:,:) 137 ! umklp(3,nbz) 138 ! The Umklapp G0-vector such as kbz + G0 = (IS) k_ibz, where kbz is in the first BZ. 139 140 real(dp) :: gmet(3,3) 141 ! Reciprocal space metric ($\textrm{bohr}^{-2}$). 142 143 real(dp) :: gprimd(3,3) 144 ! Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$) 145 146 real(dp),allocatable :: bz(:,:) 147 ! bz(3,nbz) 148 ! Points in the BZ in reduced coordinates. 149 ! TODO should packed in shells. 150 151 real(dp),allocatable :: ibz(:,:) 152 ! ibz(3,nibz) 153 ! Points in the IBZ in reduced coordinates. 154 155 real(dp),allocatable :: shift(:,:) 156 ! shift(3,nshift) 157 ! Shift for k-points, not available is nshift=0. Usually nshift=1 158 159 real(dp),allocatable :: wt(:) 160 ! wt(nibz) 161 ! Weights for each point in the IBZ. 162 163 !%real(dp),allocatable :: vbox(:) 164 ! vbox(nkbz) 165 ! Volume of the small box centered on the k-point in the full BZ. 166 ! Mainly used for inhomogeneous meshes. 167 168 complex(dpc),allocatable :: tabp(:) 169 ! tabp(nkbz) 170 ! For each point in the BZ, this table gives the phase factors associated 171 ! to non-symmorphic operations, i.e., e^{-i2\pi k_IBZ.R{^-1}t}=e^{-i2\pi k_BZ cdot t} 172 ! where \transpose R{-1}=S and (S k_IBZ)=\pm k_BZ (depending on tabi) 173 174 contains 175 176 ! Methods 177 procedure :: init => kmesh_init ! Main creation method. 178 procedure :: free => kmesh_free ! Free memory 179 procedure :: print => kmesh_print ! Printout of basic info on the object. 180 procedure :: get_bz_item => get_bz_item ! Get point in the BZ and other useful quantities. 181 procedure :: get_ibz_item => get_IBZ_item ! Get point in the IBZ and other useful quantities. 182 procedure :: get_bz_diff => get_BZ_diff ! Get the difference k1-k2 in the BZ (if any). 183 184 185 procedure :: has_bz_item => has_BZ_item ! Check if a point belongs to the BZ mesh. 186 procedure :: has_ibz_item => has_IBZ_item ! Check if a point is in the IBZ 187 procedure :: isirred => bz_mesh_isirred ! TRUE if ik_bz is in the IBZ (a non-zero umklapp is not allowed) 188 189 end type kmesh_t 190 191 public :: make_mesh ! Initialize the mesh starting from kptrlatt and shift. 192 public :: find_qmesh ! Find the Q-mesh defined as the set of all possible k1-k2 differences. 193 194 public :: isamek ! Check whether two points are equal within an umklapp G0. 195 public :: isequalk ! Check whether two points are equal within an umklapp G0 (does not report G0) 196 public :: findqg0 ! Identify q + G0 = k1-k2. 197 public :: findq ! Helper routine returning the list of q-points. 198 public :: findnq ! Helper routine returning the number of q-points. 199 public :: identk ! Find the BZ starting from the irreducible k-points. 200 public :: get_ng0sh ! Calculate the smallest box in RSpace able to treat all possible umklapp processes. 201 public :: box_len ! Return the length of the vector connecting the origin with one the faces of the unit cell.
m_bz_mesh/kpath_free [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
kpath_free
FUNCTION
Free memory allocated in the object
SOURCE
2864 subroutine kpath_free(Kpath) 2865 2866 !Arguments ------------------------------------ 2867 !scalars 2868 class(kpath_t),intent(inout) :: Kpath 2869 2870 ! ************************************************************************* 2871 2872 ABI_SFREE(Kpath%ndivs) 2873 ABI_SFREE(Kpath%bounds2kpt) 2874 ABI_SFREE(Kpath%bounds) 2875 ABI_SFREE(Kpath%points) 2876 ABI_SFREE(Kpath%dl) 2877 2878 end subroutine kpath_free
m_bz_mesh/kpath_new [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
kpath_new
FUNCTION
Create a normalized path given the extrema.
INPUTS
bounds(3,nbounds)=The points defining the path in reduced coordinates. gprimd(3,3)=Reciprocal lattice vectors ndivsm=Number of divisions to be used for the smallest segment. A negative value activates a specialized mode if with bounds is suppose to supply the full list of k-points.
SOURCE
2798 type(kpath_t) function kpath_new(bounds, gprimd, ndivsm) result(kpath) 2799 2800 !Arguments ------------------------------------ 2801 !scalars 2802 integer,intent(in) :: ndivsm 2803 !!arrays 2804 real(dp),intent(in) :: bounds(:,:),gprimd(3,3) 2805 2806 !Local variables------------------------------- 2807 integer :: ii 2808 !arrays 2809 real(dp) :: dk(3) 2810 2811 ! ************************************************************************* 2812 2813 ABI_CHECK(size(bounds, dim=1) == 3, "Wrong dim1 in bounds") 2814 !ABI_CHECK(ndivsm > 0, sjoin("ndivsm:", itoa(ndivsm))) 2815 Kpath%nbounds = size(bounds, dim=2) 2816 Kpath%ndivsm = ndivsm 2817 2818 ! Compute reciprocal space metric. 2819 Kpath%gprimd = gprimd; Kpath%gmet = matmul(transpose(gprimd), gprimd) 2820 2821 ABI_MALLOC(Kpath%ndivs, (Kpath%nbounds-1)) 2822 2823 if (kpath%ndivsm > 0) then 2824 call make_path(Kpath%nbounds, bounds, Kpath%gmet, "G", ndivsm, Kpath%ndivs, Kpath%npts, kpath%points, unit=dev_null) 2825 else 2826 ! Get list of points directly. 2827 kpath%ndivs = 1 2828 kpath%npts = kpath%nbounds 2829 ABI_MALLOC(Kpath%points, (3, Kpath%npts)) 2830 kpath%points = bounds 2831 end if 2832 2833 ABI_MALLOC(Kpath%bounds, (3, Kpath%nbounds)) 2834 Kpath%bounds = bounds 2835 2836 ! Compute distance between point i-1 and i 2837 ABI_CALLOC(kpath%dl, (kpath%npts)) 2838 do ii=2,kpath%npts 2839 dk = kpath%points(:, ii-1) - kpath%points(:,ii) 2840 kpath%dl(ii) = normv(dk, kpath%gmet, "G") 2841 end do 2842 2843 ! Mapping bounds --> points 2844 ABI_MALLOC(kpath%bounds2kpt, (kpath%nbounds)) 2845 kpath%bounds2kpt(1) = 1 2846 do ii=1,kpath%nbounds-1 2847 kpath%bounds2kpt(ii+1) = sum(kpath%ndivs(:ii)) + 1 2848 end do 2849 2850 end function kpath_new
m_bz_mesh/kpath_print [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
kpath_print
FUNCTION
Print info on the path.
INPUTS
[unit]=Unit number for output. Defaults to std_out [prtvol]=Verbosity level. [header]=String to be printed as header for additional info. [pre]=Optional string prepended to output e.g. #. Default: " "
OUTPUT
Only printing
SOURCE
2901 subroutine kpath_print(kpath, header, unit, prtvol, pre) 2902 2903 !Arguments ------------------------------------ 2904 !scalars 2905 integer,optional,intent(in) :: unit,prtvol 2906 character(len=*),optional,intent(in) :: header,pre 2907 class(kpath_t),intent(in) :: kpath 2908 2909 !Local variables------------------------------- 2910 integer :: unt,my_prtvol,ii 2911 character(len=500) :: my_pre 2912 2913 ! ************************************************************************* 2914 2915 unt = std_out; if (present(unit)) unt = unit 2916 my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol 2917 my_pre = " "; if (present(pre)) my_pre = pre 2918 if (unt <= 0) return 2919 2920 if (present(header)) write(unt,"(a)") sjoin(my_pre, '==== '//trim(adjustl(header))//' ==== ') 2921 write(unt, "(a)") sjoin(my_pre, " Number of points:", itoa(kpath%npts), ", ndivsmall:", itoa(kpath%ndivsm)) 2922 write(unt, "(a)") sjoin(my_pre, " Boundaries and corresponding index in the k-points array:") 2923 do ii=1,kpath%nbounds 2924 write(unt, "(a)") sjoin(my_pre, itoa(kpath%bounds2kpt(ii)), ktoa(kpath%bounds(:,ii))) 2925 end do 2926 write(unt, "(a)") sjoin(my_pre, " ") 2927 2928 if (my_prtvol > 10) then 2929 do ii=1,kpath%npts 2930 write(unt, "(a)") sjoin(my_pre, ktoa(kpath%points(:,ii))) 2931 end do 2932 end if 2933 2934 end subroutine kpath_print
m_bz_mesh/kpath_t [ Types ]
[ Top ] [ m_bz_mesh ] [ Types ]
NAME
path_t
FUNCTION
A (normalized) path in reciprocal space
SOURCE
215 type,public :: kpath_t 216 217 integer :: nbounds = 0 218 ! Number of extrema defining the path. 219 220 integer :: ndivsm = 0 221 ! Number of divisions used to sample the smallest segment. 222 223 integer :: npts = 0 224 ! Total number of points in the path. 225 226 real(dp) :: gprimd(3,3) 227 ! Reciprocal lattice vectors. 228 229 real(dp) :: gmet(3,3) 230 ! Metric matrix in G-space. 231 232 integer,allocatable :: ndivs(:) 233 ! ndivs(nbounds-1) 234 ! Number of divisions for each segment. 235 236 integer,allocatable :: bounds2kpt(:) 237 ! bounds2kpt(nbounds) 238 ! bounds2kpt(i): Index of the i-th extrema in the pts(:) array. 239 240 real(dp),allocatable :: bounds(:,:) 241 ! bounds(3,nbounds) 242 ! The points defining the path in reduced coordinates. 243 244 real(dp),allocatable :: points(:,:) 245 ! points(3,npts) 246 ! The points of the path in reduced coordinates. 247 248 real(dp),allocatable :: dl(:) 249 ! dl(npts) 250 ! dl(i) = Distance between the (i-1)-th and the i-th k-point. dl(1) = zero 251 252 contains 253 254 procedure :: free => kpath_free 255 ! Free memory 256 257 procedure :: print => kpath_print 258 ! Print the path. 259 260 end type kpath_t 261 262 public :: kpath_new ! Construct a new path 263 !public :: kpath_from_cryst ! High-level API to construct a path from a crystal_t 264 public :: make_path ! Construct a normalized path. TODO: Remove
m_bz_mesh/littlegroup_free_0D [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
littlegroup_free_0D
FUNCTION
Deallocate all associated pointers defined in the littlegroup_t data type.
INPUTS
Ltg=datatype to be freed
OUTPUT
SOURCE
2571 subroutine littlegroup_free_0D(Ltg) 2572 2573 !Arguments ------------------------------------ 2574 class(littlegroup_t),intent(inout) :: Ltg 2575 2576 ! ********************************************************************* 2577 2578 !@littlegroup_t 2579 ABI_SFREE(Ltg%g0) 2580 ABI_SFREE(Ltg%ibzq) 2581 ABI_SFREE(Ltg%bz2ibz) 2582 ABI_SFREE(Ltg%ibz2bz) 2583 ABI_SFREE(Ltg%igmG0) 2584 ABI_SFREE(Ltg%flag_umklp) 2585 ABI_SFREE(Ltg%preserve) 2586 ABI_SFREE(Ltg%tab) 2587 ABI_SFREE(Ltg%tabo) 2588 ABI_SFREE(Ltg%tabi) 2589 ABI_SFREE(Ltg%wtksym) 2590 2591 end subroutine littlegroup_free_0D
m_bz_mesh/littlegroup_free_1D [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
littlegroup_free_1D
FUNCTION
Deallocate all the associated pointers defined in the littlegroup_t data type.
INPUTS
Ltg=datatype to be freed
OUTPUT
SOURCE
2610 subroutine littlegroup_free_1D(Ltg) 2611 2612 !Arguments ------------------------------------ 2613 class(littlegroup_t),intent(inout) :: Ltg(:) 2614 2615 !Local variables------------------------------- 2616 integer :: ipt 2617 2618 ! ********************************************************************* 2619 2620 do ipt=1,SIZE(Ltg) 2621 call littlegroup_free_0D(Ltg(ipt)) 2622 end do 2623 2624 end subroutine littlegroup_free_1D
m_bz_mesh/littlegroup_init [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
littlegroup_init
FUNCTION
Finds symmetry operations belonging to the little group associated to an external point ext_pt and fills symmetry tables.
INPUTS
ext_pt(3)= External point in the Brillouin zone in reduce coordinated nbz=number of points in the full BZ. bz(3,nbz)=points in the full BZ. Cryst<crystal_t>= Info on symmetries and unit cell. use_umklp=flag to include umklapp G0 vectors in the definition of the little group (0:n0,1:yes) npwe=If greater than 0, the index of G-Go in the gvec(:,1:npwvec) array will be calculated and stored in %igmG0 for each symmetry preserving the external q. Note that G is one of the npwe vectors. gvec(3,npwe) coordinates of G vectors [timrev]=Optional argument to change the value of time-reversal. If not given, the value from cryst is used.
SOURCE
2253 subroutine littlegroup_init(Ltg, ext_pt, nbz, bz, Cryst, use_umklp, npwe, gvec, timrev) 2254 2255 !Arguments ------------------------------------ 2256 !scalars 2257 class(littlegroup_t),intent(inout) :: Ltg 2258 integer,intent(in) :: nbz, npwe, use_umklp 2259 type(crystal_t),target,intent(in) :: Cryst 2260 real(dp),intent(in) :: bz(3, nbz) 2261 !arrays 2262 integer,optional,intent(in) :: gvec(:,:) ! (3,npwe) 2263 real(dp),intent(in) :: ext_pt(3) 2264 integer,optional,intent(in) :: timrev 2265 2266 !Local variables------------------------------- 2267 !scalars 2268 integer :: dummy_timrev,enough,idx,ige,igpw,ik,ind,iold,iout,isym,itest,itim 2269 integer :: nkibzq,nsym,nsym_Ltg,ntest,my_timrev,ierr,npwvec 2270 real(dp) :: G0len,kin,mG0len,max_kin 2271 logical :: found,found_identity,use_antiferro 2272 character(len=500) :: msg 2273 !arrays 2274 integer :: g0(3),gg(3),gmG0(3),identity(3,3),nop(Cryst%timrev),nopg0(2) 2275 integer :: symxpt(4,2,Cryst%nsym) 2276 integer,allocatable :: indkpt1(:),symafm_ltg(:),symrec_Ltg(:,:,:),bz2ibz_smap(:,:) 2277 integer,pointer :: symafm(:),symrec(:,:,:) 2278 real(dp) :: knew(3) 2279 real(dp),allocatable :: ktest(:,:),wtk(:),wtk_folded(:) 2280 2281 !************************************************************************ 2282 2283 ABI_CHECK(any(cryst%timrev == [1, 2]), sjoin("Wrong value for cryst%timrev:", itoa(cryst%timrev))) 2284 2285 ! Destroy structure if it already exists 2286 call Ltg%free() 2287 2288 ! Copy useful data. 2289 nsym = Cryst%nsym 2290 my_timrev = Cryst%timrev 2291 if (present(timrev)) then 2292 my_timrev = timrev 2293 ABI_CHECK_ILEQ(my_timrev, cryst%timrev, "my_timrev cannot be greater that cryst%timrev") 2294 end if 2295 2296 symrec => Cryst%symrec 2297 symafm => Cryst%symafm 2298 use_antiferro = Cryst%use_antiferro 2299 2300 ! Store dimensions and useful info. 2301 Ltg%nsym_sg = nsym 2302 Ltg%timrev = my_timrev 2303 Ltg%nbz = nbz 2304 !Ltg%use_umklp=use_umklp ! 1 if umklapp processes are used 2305 Ltg%ext_pt(:)=ext_pt(:) 2306 2307 ABI_MALLOC(Ltg%G0, (3,2,nsym)) 2308 ABI_MALLOC(Ltg%ibzq, (nbz)) 2309 ABI_MALLOC(Ltg%bz2ibz, (nbz)) 2310 ABI_MALLOC(Ltg%preserve, (2, nsym)) 2311 ABI_MALLOC(Ltg%wtksym, (2, nsym, nbz)) 2312 ABI_MALLOC(Ltg%tab, (nbz)) 2313 ABI_MALLOC(Ltg%tabi, (nbz)) 2314 ABI_MALLOC(Ltg%tabo, (nbz)) 2315 ABI_MALLOC(Ltg%flag_umklp, (2, nsym)) 2316 2317 ! In the old GW implementation we were removing symmetries related by time-reversal and 2318 ! sometimes it happened that only the inversion was reported in the KSS file (see outkss.F90). 2319 identity(:,:)=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/)) ; found_identity=.FALSE. 2320 do isym=1,nsym 2321 if (ALL(symrec(:,:,isym) == identity)) then 2322 found_identity=.TRUE.; EXIT 2323 end if 2324 end do 2325 if (.not. found_identity) then 2326 write(msg,'(5a)')& 2327 'Only the inversion was found in the set of symmetries read from the KSS file ',ch10,& 2328 'Likely you are using a KSS file generated with an old version of Abinit, ',ch10,& 2329 'To run a GW calculation with an old KSS file, use version < 5.5 ' 2330 ABI_ERROR(msg) 2331 end if 2332 2333 ! Find operations in the little group as well as umklapp vectors G0 2334 call littlegroup_q(nsym, ext_pt, symxpt, symrec, symafm, dummy_timrev, prtvol=0) 2335 2336 Ltg%preserve(:,:)=0; Ltg%g0(:,:,:)=0; Ltg%flag_umklp(:,:)=0; mG0len=zero 2337 2338 do itim=1,my_timrev 2339 do isym=1,nsym 2340 if (symafm(isym)==-1) CYCLE 2341 2342 if (symxpt(4, itim, isym) == 1) then !\pm Sq = q+g0 2343 if (ANY(symxpt(1:3, itim, isym) /= 0) .and. use_umklp == 0) CYCLE ! Exclude non zero G0 vectors 2344 Ltg%preserve(itim, isym) = 1 2345 g0(:)=symxpt(1:3, itim, isym); Ltg%g0(:, itim, isym) = g0(:) 2346 if (ANY(Ltg%g0(:, itim, isym) /= 0)) Ltg%flag_umklp(itim, isym) = 1 2347 ! Max radius to be considered to include all G0s 2348 G0len = normv(g0,Cryst%gmet,'G') 2349 mG0len = MAX(mG0len, G0len) 2350 end if 2351 end do 2352 end do 2353 2354 nop(:) = 0; nopg0(:) = 0 2355 do itim=1,my_timrev 2356 nop (itim) = SUM(Ltg%preserve (itim,:)) 2357 nopg0(itim) = SUM(Ltg%flag_umklp(itim,:)) 2358 end do 2359 nsym_Ltg = SUM(nop(:)) 2360 2361 ! Store little group operations, include time-reversal if present. 2362 Ltg%nsym_Ltg = nsym_Ltg 2363 ABI_MALLOC(symrec_Ltg, (3, 3, Ltg%nsym_Ltg)) 2364 2365 ind = 1 2366 do itim=1,my_timrev 2367 do isym=1,nsym 2368 if (Ltg%preserve(itim,isym)==1) then 2369 if (itim==1) symrec_Ltg(:,:,ind) = symrec(:,:,isym) 2370 if (itim==2) symrec_Ltg(:,:,ind) =-symrec(:,:,isym) 2371 ind = ind+1 2372 end if 2373 end do 2374 end do 2375 2376 ! Check the closure of the (ferromagnetic) little group 2377 ABI_MALLOC(symafm_ltg,(Ltg%nsym_Ltg)) 2378 symafm_ltg(:) = 1 2379 call sg_multable(Ltg%nsym_Ltg,symafm_ltg,symrec_Ltg,ierr) 2380 ABI_CHECK(ierr == 0, "Error in group closure") 2381 2382 ABI_FREE(symafm_ltg) 2383 2384 ! Find the irreducible zone associated to ext_pt 2385 ! Do not use time-reversal since it has been manually introduced previously 2386 ABI_MALLOC(indkpt1, (nbz)) 2387 ABI_MALLOC(wtk_folded, (nbz)) 2388 ABI_MALLOC(wtk, (nbz)) 2389 ABI_MALLOC(bz2ibz_smap, (6, nbz)) 2390 wtk=one; iout=0; dummy_timrev=0 2391 2392 call symkpt(0,Cryst%gmet,indkpt1,iout, bz, nbz, nkibzq, Ltg%nsym_Ltg, symrec_Ltg, dummy_timrev, wtk, wtk_folded, & 2393 bz2ibz_smap, xmpi_comm_self) 2394 2395 ABI_FREE(bz2ibz_smap) 2396 ABI_FREE(indkpt1) 2397 ABI_FREE(wtk) 2398 2399 Ltg%nibz_Ltg = nkibzq 2400 ! 2401 ! === Set up table in the BZ === 2402 ! * 0 if the point does not belong to IBZ_xpt, 1 otherwise 2403 ABI_MALLOC(Ltg%ibz2bz,(nkibzq)) 2404 Ltg%ibzq(:)=0; Ltg%bz2ibz(:)=0; Ltg%ibz2bz(:)=0 2405 2406 ind=0; enough=0 2407 do ik=1,nbz 2408 if (wtk_folded(ik)>tol8) then 2409 ind = ind + 1 2410 Ltg%ibzq(ik) = 1 2411 Ltg%bz2ibz(ik) = ind 2412 Ltg%ibz2bz(ind)= ik 2413 end if 2414 end do 2415 ABI_CHECK_IEQ(ind, Ltg%nibz_Ltg, "BUG: ind /= Ltg%nibz_Ltg") 2416 2417 ! Reconstruct full BZ starting from IBZ_q. 2418 ! Calculate appropriate weight for each item (point,symmetry operation,time-reversal) 2419 Ltg%tab=0; Ltg%tabo=0; Ltg%tabi=0; Ltg%wtksym(:,:,:)=0 2420 2421 ! Start with zero no. of k-points found 2422 ntest = 0 2423 ABI_MALLOC(ktest, (3, nbz)) 2424 ktest = zero 2425 2426 do ik=1,nbz 2427 if (Ltg%ibzq(ik) /= 1) CYCLE 2428 ! * Loop over symmetry operations S and time-reversal. 2429 ! * Use spatial inversion instead of time reversal whenever possible. 2430 do itim=1,my_timrev 2431 do isym=1,nsym 2432 2433 ! Form IS k only for (IS) pairs in the (ferromagnetic) little group. 2434 if (symafm(isym)==-1) CYCLE 2435 if (Ltg%preserve(itim,isym)==0) CYCLE 2436 knew(:)=(3-2*itim)*MATMUL(symrec(:,:,isym), bz(:,ik)) 2437 ! 2438 ! Check whether it has already been found (to within a RL vector) 2439 iold=0 2440 do itest=1,ntest 2441 if (isamek(knew(:),ktest(:,itest),gg)) iold=iold+1 2442 end do 2443 2444 if (iold==0) then 2445 ! Found new BZ point 2446 ! For this point the operation (isym,itim) must be considered to reconstruct the full BZ 2447 Ltg%wtksym(itim,isym,ik)=1 2448 ntest=ntest+1 2449 ktest(:,ntest)=knew(:) 2450 ! 2451 ! Now find knew in the BZ array 2452 found=.FALSE. 2453 do idx=1,nbz 2454 if (isamek(knew(:), bz(:,idx), gg)) then ! They are the same within a RL vector 2455 Ltg%tab (idx)=ik 2456 Ltg%tabo(idx)=isym 2457 Ltg%tabi(idx)=3-2*itim 2458 found=.TRUE.; EXIT 2459 end if 2460 end do 2461 if (.not.found) then 2462 write(msg,'(a,3f12.6,a)')'Not able to find the ',knew(:),' in the array BZ ' 2463 ABI_ERROR(msg) 2464 end if 2465 end if 2466 2467 end do !isym 2468 end do !itim 2469 2470 end do !nbz 2471 2472 ABI_FREE(ktest) 2473 2474 if (ntest/=nbz) then 2475 ABI_BUG(sjoin('ntest - nbz = ',itoa(ntest-nbz))) 2476 end if 2477 2478 if (sum(Ltg%wtksym) /= nbz) then 2479 ABI_BUG(sjoin('sum(Ltg%wtksym)-nbz = ', itoa(SUM(Ltg%wtksym)-nbz))) 2480 end if 2481 2482 Ltg%max_kin_gmG0=zero 2483 2484 if (npwe > 0.and. PRESENT(gvec)) then 2485 npwvec = SIZE(gvec,DIM=2) 2486 ! This correspond to the case in which we need to know the index of G-Go in the gvec array 2487 ! where G is one of the npwe vectors. This is required in screening but not in sigma. 2488 ! The drawback is that the effective G sphere used to calculate the oscillators must be smaller 2489 ! that gvec if we want to avoid possible aliasing effects. Lifting this constraint would require 2490 ! a lot of boring coding. (no need to do this if ext_pt=zero, but oh well) 2491 ABI_MALLOC(Ltg%igmG0,(npwe,2,nsym)) 2492 Ltg%igmG0(:,:,:)=0 2493 max_kin=zero 2494 2495 ! Loop over symmetry operations S and time-reversal 2496 do itim=1,my_timrev 2497 do isym=1,nsym 2498 ! Form IS k only for (IS) pairs in the little group 2499 if (symafm(isym)==-1) CYCLE 2500 if (Ltg%preserve(itim,isym)/=0) then 2501 g0(:)=Ltg%g0(:,itim,isym) 2502 do ige=1,npwe 2503 gmG0(:)=gvec(:,ige)-g0(:) 2504 kin=half*normv(gmG0,Cryst%gmet,'G')**2 2505 max_kin=MAX(max_kin,kin) 2506 2507 found=.FALSE. 2508 do igpw=1,npwvec 2509 if (ALL(gvec(:,igpw)-gmG0(:)==0)) then 2510 Ltg%igmG0(ige,itim,isym)=igpw 2511 found=.TRUE.; EXIT 2512 end if 2513 end do 2514 if (.not. found) then 2515 write(msg,'(5a,f8.3,2a,3i5)')& 2516 'Not able to found G-G0 in the largest G-spere ',ch10,& 2517 'Decrease the size of epsilon or, if possible, increase ecutwfn (>ecuteps) ',ch10,& 2518 'Minimum required cutoff energy for G-G0 sphere= ',kin,ch10,& 2519 'G0 = ',g0(:) 2520 ABI_ERROR(msg) 2521 end if 2522 end do 2523 end if 2524 end do 2525 end do 2526 Ltg%max_kin_gmG0=max_kin 2527 end if 2528 ABI_FREE(symrec_Ltg) 2529 2530 ! DEBUG SECTION 2531 !do ik=1,nbz 2532 ! if (ABS(SUM(Ltg%wtksym(1,:,ik)+Ltg%wtksym(2,:,ik))-wtk_folded(ik))>tol6) then 2533 ! write(std_out,*)' sum(Ltg%wtksym,ik)-wtk_folded(ik) = ',sum(Ltg%wtksym(1,:,ik)+Ltg%wtksym(2,:,ik))-wtk_folded(ik) 2534 ! write(std_out,*)Ltg%wtksym(1,:,ik),Ltg%wtksym(2,:,ik),wtk_folded(ik) 2535 ! write(std_out,*)ik, bz(:,ik) 2536 ! ABI_BUG("Wrong weight") 2537 ! end if 2538 !end do 2539 !do ik=1,nbz 2540 ! knew = Ltg%tabi(ik) * MATMUL(symrec(:,:,Ltg%tabo(ik)), bz(:,Ltg%tab(ik))) 2541 ! if (.not.isamek(knew, bz(:,ik),gg)) then 2542 ! write(std_out,*)knew, bz(:,ik) 2543 ! write(std_out,*)Ltg%tabo(ik),Ltg%tabi(ik),Ltg%tab(ik) 2544 ! ABI_BUG("Wrong tables") 2545 ! end if 2546 !end do 2547 2548 ABI_FREE(wtk_folded) 2549 2550 DBG_EXIT("COLL") 2551 2552 end subroutine littlegroup_init
m_bz_mesh/littlegroup_print [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
littlegroup_print
FUNCTION
Print info on the littlegroup_t data type.
INPUTS
Ltg=the datatype to be printed [unit]=the unit number for output [prtvol]=verbosity level [mode_paral]=either "COLL" or "PERS"
OUTPUT
Only printing
SOURCE
2647 subroutine littlegroup_print(Ltg, unit, prtvol, mode_paral) 2648 2649 !Arguments ------------------------------------ 2650 class(littlegroup_t),intent(in) :: Ltg 2651 integer,optional,intent(in) :: prtvol,unit 2652 character(len=4),optional,intent(in) :: mode_paral 2653 2654 !Local variables------------------------------- 2655 !scalars 2656 integer :: itim,my_unt,my_prtvol 2657 character(len=4) :: my_mode 2658 character(len=500) :: msg 2659 !arrays 2660 integer :: nop(Ltg%timrev),nopg0(Ltg%timrev) 2661 2662 ! ********************************************************************* 2663 2664 my_unt =std_out; if (PRESENT(unit )) my_unt =unit 2665 my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol 2666 my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral 2667 2668 write(msg,'(7a,i0,a,i0,2a,i0,a,i0)')ch10, & 2669 ' ==== Little Group Info ==== ', ch10, & 2670 ' External point: ',trim(ktoa(Ltg%ext_pt)), ch10, & 2671 ' Number of points in the IBZ defined by little group: ', Ltg%nibz_Ltg, '/', Ltg%nbz,ch10, & 2672 ' Number of operations in the little group: ',Ltg%nsym_Ltg,'/',Ltg%nsym_sg 2673 call wrtout(my_unt,msg,my_mode) 2674 2675 nop=0 ; nopg0=0 2676 do itim=1,Ltg%timrev 2677 nop (itim)=SUM(Ltg%preserve (itim,:)) 2678 nopg0(itim)=SUM(Ltg%flag_umklp(itim,:)) 2679 end do 2680 2681 do itim=1,Ltg%timrev 2682 if (itim==1) then 2683 write(msg,'(2(a,i2,a))') & 2684 ' No time-reversal symmetry with zero umklapp: ',nop(1)-nopg0(1),ch10,& 2685 ' No time-reversal symmetry with non-zero umklapp: ',nopg0(1),ch10 2686 call wrtout(my_unt,msg,my_mode) 2687 else if (itim==2) then 2688 write(msg,'(2(a,i2,a))') & 2689 ' time-reversal symmetry with zero umklapp: ',nop(2)-nopg0(2),ch10,& 2690 ' time-reversal symmetry with non-zero umklapp: ',nopg0(2),ch10 2691 call wrtout(my_unt,msg,my_mode) 2692 end if 2693 end do 2694 2695 end subroutine littlegroup_print
m_bz_mesh/littlegroup_t [ Types ]
[ Top ] [ m_bz_mesh ] [ Types ]
NAME
littlegroup_t
FUNCTION
For the GW part of ABINIT. The littlegroup_t structured datatype gather information on the little group associated to an external vector q. The little group associated to q is defined as the subset of the space group that preserves q, modulo a G0 vector (also called umklapp vector). Namely Sq = q +G0, where S is an operation in reciprocal space. If time reversal symmetry holds true, it is possible to enlarge the little group by including the operations such as -Sq = q+ G0. The operations belongin to the little group define an irriducible wedge in the Brillouin zone that is, usually, larger than the irredubile zone defined by the space group. The two zone coincide when q=0
TODO
Rationalize most of the arrays, in particular the tables This structure shoud be rewritten almost from scratch, thus avoid using it for your developments.
SOURCE
296 type,public :: littlegroup_t 297 298 integer :: npw ! No. of planewaves used to describe the wavefuntion, used to dimension igmG0 299 integer :: nsym_sg ! No. of operations in the space group (*NOT* the little group) 300 integer :: nsym_ltg ! No. of symmetry operations in the little group (time-reversal is included, if can be used) 301 integer :: timrev ! 2 if time-reversal is considered, 1 otherwise 302 integer :: nbz ! No. of kpoints in the full BZ 303 integer :: nibz_ltg ! No. of points in the irreducible wedge defined by the little group 304 !integer :: use_umklp ! 1 if umklapp processes are included 305 306 real(dp) :: max_kin_gmG0 307 ! Max kinetic energy of G-G0 in case of umklapp. 308 309 integer,allocatable :: G0(:,:,:) 310 ! (3,2,nsym_sg) 311 ! Reduced coordinates of the umklapp G0 vector. 312 313 integer,allocatable :: ibzq(:) 314 ! ibzq(nbz) 315 ! 1 if the point belongs to the IBZ_q defined by ext_pt, 0 otherwise. 316 317 integer,allocatable :: bz2ibz(:) 318 ! bz2ibz(nbz) 319 ! Index of the point in the irreducible wedge defined by the little group, 0 otherwise. 320 321 integer,allocatable :: ibz2bz(:) 322 ! ibz2bz(nibz_ltg) 323 ! The correspondind index in the BZ array 324 325 integer,allocatable :: igmG0(:,:,:) 326 ! iumklp(npw,2,nsym_sg) 327 ! Index of G-G0 in the FFT array for each operations IS (I=\pm 1). 328 329 integer,allocatable :: flag_umklp(:,:) 330 ! flag_umklp(2,nsym_sg) 331 ! 1 if the operation IS requires a non null G0 vector to preserve q, 0 otherwise. 332 333 integer,allocatable :: preserve(:,:) 334 ! preserve(2, nsym_sg) 335 ! (1,S) is 1 if the operation S in rec space preserves the external q-point i.e Sq=q+G0 336 ! (2,S) is 1 if -Sq=q+G0. G0 is a reciprocal lattice vector also called "umklapp vector". 337 338 integer,allocatable :: tab(:) 339 ! tab(nbz) 340 ! For each point in BZ, the index of the irreducible point (kIBZ_q) in the irreducible 341 ! wedge defined by the little group of q. kBZ= (IS) kIBZ where I is the inversion or the identity. 342 343 integer,allocatable :: tabo(:) 344 ! tabo(nbz) 345 ! The index of the operation S in the little group that rotates kIBZ_q into \pm kBZ. 346 347 integer,allocatable :: tabi(:) 348 ! tabi(nbz) 349 ! for each k-point in the BZ defines whether inversion has to be 350 ! considered in the relation kBZ= IS kIBZ_q (1 => only S; -1 => -S). 351 352 integer,allocatable :: wtksym(:,:,:) 353 ! (2, nsym_sg, kbz) 354 ! 1 if IS belongs to the little group, 0 otherwise TODO (should invert the first two dimensions) 355 356 real(dp) :: ext_pt(3) 357 ! The external point defining the little group. 358 359 contains 360 361 procedure :: init => littlegroup_init 362 procedure :: print => littlegroup_print 363 procedure :: free => littlegroup_free_0D 364 365 end type littlegroup_t 366 367 public :: littlegroup_free
m_bz_mesh/make_mesh [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
make_mesh
FUNCTION
Initialize the kmesh_t starting from kptrlatt and shiftk
INPUTS
Cryst<crystal_t>=Info on the crystalline structure. nshiftk=Number of shifts for the mesh. kptrlatt(3,3)= Coordinates of three vectors in real space, expressed in reduced coordinates. They define a super-lattice in real space. The k point lattice is the reciprocal of this super-lattice, eventually shifted by shift. shiftk(3,nshiftk)=Shifts for the k-mesh. [vacuum(3)]=For each direction, 0 if no vacuum, 1 if vacuum
OUTPUT
Kmesh<kmesh_t>=Object gathering info on the sampling of the Brillouin zone.
SOURCE
1275 subroutine make_mesh(Kmesh,Cryst,kptopt,kptrlatt,nshiftk,shiftk,& 1276 & vacuum,break_symmetry) ! Optional 1277 1278 !Arguments ------------------------------- 1279 !scalars 1280 integer,intent(in) :: nshiftk,kptopt 1281 logical,optional,intent(in) :: break_symmetry 1282 type(kmesh_t),intent(inout) :: Kmesh 1283 type(crystal_t),intent(in) :: Cryst 1284 !arrays 1285 integer,intent(inout) :: kptrlatt(3,3) 1286 integer,optional,intent(in) :: vacuum(3) 1287 real(dp),intent(in) :: shiftk(3,nshiftk) 1288 1289 !Local variables ------------------------- 1290 !scalars 1291 integer,parameter :: chksymbreak0=0 1292 integer :: iscf,nkbz,nkibz,nkpt_computed,my_nshiftk 1293 real(dp) :: kptrlen 1294 logical :: my_break_symmetry 1295 !arrays 1296 integer :: my_vacuum(3) 1297 real(dp),allocatable :: kibz(:,:),wtk(:),my_shiftk(:,:),ref_kbz(:,:) 1298 1299 ! ************************************************************************* 1300 1301 DBG_ENTER("COLL") 1302 1303 !@kmesh_t 1304 !if (ALL(kptopt/=(/1,2,3,4/))) then 1305 if (ALL(kptopt/=(/1,3/))) then 1306 ABI_WARNING(sjoin(" Not allowed value for kptopt: ", itoa(kptopt))) 1307 end if 1308 ! 1309 ! ====================================================================== 1310 ! ==== First call to getkgrid to obtain nkibz as well as the BZ set ==== 1311 ! ====================================================================== 1312 iscf=7 ! use for the Weights in NSCF calculation. check it more carefully. 1313 nkibz=0 ! Compute number of k-points in the BZ and IBZ 1314 1315 my_vacuum = (/0,0,0/); if (PRESENT(vacuum)) my_vacuum=vacuum 1316 1317 my_nshiftk = nshiftk 1318 ABI_CHECK(my_nshiftk>0.and.my_nshiftk<=MAX_NSHIFTK, sjoin("Wrong nshiftk must be between 1 and ", itoa(MAX_NSHIFTK))) 1319 ABI_MALLOC(my_shiftk, (3, MAX_NSHIFTK)) 1320 my_shiftk=zero; my_shiftk(:,1:nshiftk) = shiftk(:,:) 1321 1322 !write(std_out,*)" In make_mesh" 1323 !write(std_out,*)" kptopt ",kptopt," kptrlatt ",kptrlatt 1324 !write(std_out,*)" nshiftk ",nshiftk," shiftk ",shiftk 1325 1326 ABI_MALLOC(kibz,(3,nkibz)) 1327 ABI_MALLOC(wtk,(nkibz)) 1328 1329 call getkgrid(chksymbreak0,0,iscf,kibz,kptopt,kptrlatt,kptrlen,Cryst%nsym,0,nkibz,my_nshiftk,& 1330 & Cryst%nsym,Cryst%rprimd,my_shiftk,Cryst%symafm,Cryst%symrel,my_vacuum,wtk,fullbz=ref_kbz) 1331 1332 nkbz = SIZE(ref_kbz,DIM=2) 1333 1334 ABI_FREE(kibz) 1335 ABI_FREE(wtk) 1336 1337 !write(std_out,*)" after getkgrid1: nkbz = ",nkbz," nkibz=",nkibz 1338 !write(std_out,*)" ref_kbz = ",ref_kbz 1339 1340 ! ============================================================== 1341 ! ==== Recall getkgrid to get kibz(3,nkibz) and wtk(nkibz) ===== 1342 ! ============================================================== 1343 1344 ABI_MALLOC(kibz,(3,nkibz)) 1345 ABI_MALLOC(wtk,(nkibz)) 1346 1347 call getkgrid(chksymbreak0,0,iscf,kibz,kptopt,kptrlatt,kptrlen,Cryst%nsym,nkibz,nkpt_computed,my_nshiftk,& 1348 & Cryst%nsym,Cryst%rprimd,my_shiftk,Cryst%symafm,Cryst%symrel,my_vacuum,wtk) 1349 1350 ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ. 1351 Kmesh%nshift = my_nshiftk 1352 Kmesh%kptrlatt = kptrlatt 1353 1354 ! Call the main creation method to get the tables tabo, tabi, tabp, umklp... 1355 ! kmesh_init will reconstruct the BZ from kibz but pruning the k-points not in ref_bz 1356 ! TODO: solve problem with timrev 1357 my_break_symmetry=.FALSE.; if (PRESENT(break_symmetry)) my_break_symmetry=break_symmetry 1358 call kmesh_init(Kmesh,Cryst,nkibz,kibz,kptopt,ref_bz=ref_kbz,break_symmetry=my_break_symmetry) 1359 1360 ABI_FREE(ref_kbz) 1361 ABI_FREE(kibz) 1362 ABI_FREE(wtk) 1363 1364 ABI_MALLOC(Kmesh%shift,(3,my_nshiftk)) 1365 Kmesh%shift=my_shiftk(:,1:my_nshiftk) 1366 ! Init Kmesh is breaking nshiftk 1367 Kmesh%nshift=my_nshiftk 1368 1369 ABI_FREE(my_shiftk) 1370 1371 DBG_EXIT("COLL") 1372 1373 end subroutine make_mesh
m_bz_mesh/make_path [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
make_path
FUNCTION
Generate a normalized path given the extrema. See also kpath_t and kpath_new (recommended API).
INPUTS
nbounds=Number of extrema defining the path. bounds(3,nbounds)=The points defining the path in reduced coordinates. met(3,3)=Metric matrix. space='R' for real space, G for reciprocal space. ndivsm=Number of divisions to be used for the smallest segment. [unit]=Fortran unit for formatted output. Default: dev_null
OUTPUT
npts=Total number of points in the normalized circuit.
ndivs(nbounds-1)=Number of division for each segment
path: allocated inside the routine. When the subroutine returns, path(3,npts) will
contain the path in reduced coordinates.
SOURCE
1783 subroutine make_path(nbounds, bounds, met, space, ndivsm, ndivs, npts, path, unit) 1784 1785 !Arguments ------------------------------------ 1786 !scalars 1787 integer,intent(in) :: nbounds,ndivsm 1788 integer,optional,intent(in) :: unit 1789 integer,intent(out) :: npts 1790 character(len=1),intent(in) :: space 1791 !arrays 1792 integer,intent(out) :: ndivs(nbounds-1) 1793 real(dp),intent(in) :: bounds(3,nbounds),met(3,3) 1794 real(dp),allocatable,intent(out) :: path(:,:) 1795 1796 !Local variables------------------------------- 1797 !scalars 1798 integer,parameter :: prtvol=0 1799 integer :: idx,ii,jp,ount 1800 real(dp) :: nfact 1801 character(len=500) :: msg 1802 !arrays 1803 real(dp) :: diff(3),lng(nbounds-1) 1804 1805 ! ************************************************************************* 1806 1807 ABI_CHECK(ndivsm > 0, sjoin('ndivsm', itoa(ndivsm))) 1808 1809 ount = dev_null; if (present(unit)) ount = unit 1810 1811 do ii=1,nbounds-1 1812 diff(:)=bounds(:,ii+1)-bounds(:,ii) 1813 lng(ii) = normv(diff,met,space) 1814 end do 1815 1816 ! Avoid division by zero if any k(:,i+1)=k(:,i). 1817 nfact=MINVAL(lng) 1818 if (ABS(nfact)<tol6) then 1819 write(msg,'(3a)')& 1820 'Found two equivalent consecutive points in the path ',ch10,& 1821 'This is not allowed, modify the path in your input file' 1822 ABI_ERROR(msg) 1823 end if 1824 1825 nfact=nfact/ndivsm 1826 ndivs(:)=NINT(lng(:)/nfact) 1827 npts=SUM(ndivs)+1 !1 for the first point 1828 1829 write(msg,'(2a,i0,2a)')ch10,& 1830 ' Total number of points in the path: ',npts,ch10,& 1831 ' Number of divisions for each segment of the normalized path: ' 1832 call wrtout(ount,msg) 1833 1834 do ii=1,nbounds-1 1835 write(msg,'(2(3f8.5,a),i0,a)')bounds(:,ii),' ==> ',bounds(:,ii+1),' ( ndivs : ',ndivs(ii),' )' 1836 call wrtout(ount,msg) 1837 end do 1838 call wrtout(ount,ch10) 1839 1840 ! Allocate and construct the path. 1841 ABI_MALLOC(path,(3,npts)) 1842 1843 if (prtvol > 0) call wrtout(ount,' Normalized Path: ') 1844 idx=0 1845 do ii=1,nbounds-1 1846 do jp=1,ndivs(ii) 1847 idx=idx+1 1848 path(:,idx)=bounds(:,ii)+(jp-1)*(bounds(:,ii+1)-bounds(:,ii))/ndivs(ii) 1849 if (prtvol > 0) then 1850 write(msg,'(i4,4x,3(f8.5,1x))')idx,path(:,idx) 1851 call wrtout(ount,msg) 1852 end if 1853 end do 1854 end do 1855 path(:,npts)=bounds(:,nbounds) 1856 1857 if (prtvol > 0) then 1858 write(msg,'(i0,4x,3(f8.5,1x))')npts,path(:,npts) 1859 call wrtout(ount,msg) 1860 end if 1861 1862 end subroutine make_path
m_bz_mesh/setup_k_rotation [ Functions ]
[ Top ] [ m_bz_mesh ] [ Functions ]
NAME
setup_k_rotation
FUNCTION
Set up tables giving the correspondence btw a k-point and its rotated image.
INPUTS
timrev=2 if time-reversal can be used, 1 otherwise. nsym=Number of symmetry operations symrec(3,3,nsym)=Symmetry operations in reciprocal space in reduced coordinates. nbz=Number of k-points kbz(3,nbz)=k-points in reduced coordinates. gmet(3,3)=Metric in reciprocal space.
OUTPUT
krottb(k,I,S)=Index of (IS) k in the array bz
krottbm1(k,I,S)=Index of IS^{-1} k
SOURCE
717 subroutine setup_k_rotation(nsym,timrev,symrec,nbz,kbz,gmet,krottb,krottbm1) 718 719 !Arguments ------------------------------------ 720 !scalars 721 integer,intent(in) :: nbz,nsym,timrev 722 !arrays 723 integer,intent(in) :: symrec(3,3,nsym) 724 integer,intent(out) :: krottb(nbz,timrev,nsym),krottbm1(nbz,timrev,nsym) 725 real(dp),intent(in) :: kbz(3,nbz),gmet(3,3) 726 727 !Local variables ------------------------------ 728 !scalars 729 integer :: ik,ikp,isym,itim,nsh,ik_st !,sh_start 730 real(dp),parameter :: KTOL=tol6 731 real(dp) :: norm_old,norm !,norm_rot !norm_base 732 logical :: found,isok 733 character(len=500) :: msg 734 !arrays 735 integer :: g0(3) 736 integer :: iperm(nbz),shlim(nbz+1) 737 real(dp) :: shift(3),kbase(3),krot(3),knorm(nbz),kwrap(3),shlen(nbz+1) 738 739 !************************************************************************ 740 741 DBG_ENTER("COLL") 742 743 ! Sort the k-points according to their norm to speed up the search below. 744 do ik=1,nbz 745 call wrap2_pmhalf(kbz(1,ik),kwrap(1),shift(1)) 746 call wrap2_pmhalf(kbz(2,ik),kwrap(2),shift(2)) 747 call wrap2_pmhalf(kbz(3,ik),kwrap(3),shift(3)) 748 knorm(ik) = normv(kwrap,gmet,"G") 749 iperm(ik)= ik 750 end do 751 752 call sort_dp(nbz,knorm,iperm,KTOL) 753 ! 754 ! The index of the initial (sorted) k-point in each shell 755 nsh=1; norm_old=knorm(1) 756 757 shlim(1)=1; shlen(1) = norm_old 758 do ik_st=2,nbz 759 norm = knorm(ik_st) 760 if (ABS(norm-norm_old) > KTOL) then 761 norm_old = norm 762 nsh = nsh+1 763 shlim(nsh) = ik_st 764 shlen(nsh) = norm 765 end if 766 end do 767 shlim(nsh+1)=nbz+1 768 shlen(nsh+1)=HUGE(one) 769 ! 770 ! === Set up k-rotation tables === 771 ! * Use spatial inversion instead of time reversal whenever possible. 772 !call wrtout(std_out," Begin sorting ","COLL") 773 774 isok=.TRUE. 775 do ik=1,nbz 776 kbase(:)=kbz(:,ik) 777 ! 778 do itim=1,timrev 779 do isym=1,nsym 780 krot(:)=(3-2*itim)*MATMUL(symrec(:,:,isym),kbase) 781 782 found=.FALSE. 783 #if 1 784 ! Old code 785 do ikp=1,nbz 786 if (isamek(krot,kbz(:,ikp),g0)) then 787 found=.TRUE. 788 krottb (ik ,itim,isym)=ikp 789 krottbm1(ikp,itim,isym)=ik 790 EXIT 791 end if 792 end do 793 #else 794 ! 795 ! Locate the shell index with bisection. 796 call wrap2_pmhalf(krot(1),kwrap(1),shift(1)) 797 call wrap2_pmhalf(krot(2),kwrap(2),shift(2)) 798 call wrap2_pmhalf(krot(3),kwrap(3),shift(3)) 799 norm_rot = normv(kwrap,gmet,"G") 800 sh_start = bisect(shlen(1:nsh+1),norm_rot) 801 802 do ik_st=shlim(sh_start),nbz 803 ikp = iperm(ik_st) 804 if (isamek(krot,kbz(:,ikp),g0)) then 805 found=.TRUE. 806 krottb (ik ,itim,isym)=ikp 807 krottbm1(ikp,itim,isym)=ik 808 !write(std_out,*)ik_st,shlim(sh_start),nbz 809 EXIT 810 end if 811 end do 812 #endif 813 if (.not.found) then 814 isok=.FALSE. 815 !write(std_out,*)" norm_base,norm_rot ",norm_base,norm_rot 816 !write(std_out,*)normv(kbase,gmet,"G"),normv(krot,gmet,"G") 817 write(msg,'(2(a,i4),2x,2(3f12.6,2a),i3,a,i2)')& 818 & 'Initial k-point ',ik,'/',nbz,kbase(:),ch10,& 819 & 'Rotated k-point (not found) ',krot(:),ch10,& 820 & 'Through symmetry operation ',isym,' and itim ',itim 821 ABI_ERROR(msg) 822 end if 823 824 end do 825 end do 826 end do 827 828 if (.not.isok) then 829 ABI_ERROR('k-mesh not closed') 830 end if 831 832 DBG_EXIT("COLL") 833 834 end subroutine setup_k_rotation