TABLE OF CONTENTS


ABINIT/m_bz_mesh [ Modules ]

[ Top ] [ 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-2018 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

PARENTS

CHILDREN

SOURCE

41 #if defined HAVE_CONFIG_H
42 #include "config.h"
43 #endif
44 
45 #include "abi_common.h"
46 
47 MODULE m_bz_mesh
48 
49  use defs_basis
50  use m_errors
51  use m_abicore
52  use m_sort
53 
54  use m_fstrings,       only : ltoa, itoa, sjoin, ktoa
55  use m_numeric_tools,  only : is_zero, isinteger, imin_loc, imax_loc, bisect, wrap2_pmhalf
56  use m_symtk,          only : chkgrp, littlegroup_q
57  use m_geometry,       only : normv
58  use m_crystal,        only : crystal_t
59  use m_kpts,           only : getkgrid
60  use m_symkpt,     only : symkpt
61 
62  implicit none
63 
64  private
65 
66  real(dp),parameter :: TOL_KDIFF=0.0001_dp
67  ! Tolerance below which two points are considered equal within a RL vector:
68  ! for each reduced direction the absolute difference between the coordinates must be less that TOL_KDIFF
69 
70  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.

PARENTS

CHILDREN

SOURCE

3156 function box_len(qpt,gprimd)
3157 
3158 
3159 !This section has been created automatically by the script Abilint (TD).
3160 !Do not modify the following lines by hand.
3161 #undef ABI_FUNC
3162 #define ABI_FUNC 'box_len'
3163 !End of the abilint section
3164 
3165  implicit none
3166 
3167 !Arguments ------------------------------------
3168 !scalars
3169  real(dp) :: box_len
3170 !arrays
3171  real(dp),intent(in) :: qpt(3),gprimd(3,3)
3172 
3173 !Local variables-------------------------------
3174 !scalars
3175  integer :: idir,iplane
3176  real(dp) :: x1,x2,x3
3177 !arrays
3178  real(dp) :: my_qpt(3),gmet(3,3),q0box(3)
3179 
3180 ! *************************************************************************
3181 
3182 !Compute reciprocal space metric
3183  gmet = MATMUL(TRANSPOSE(gprimd),gprimd)
3184 
3185 !Rotate the input q-point such that it is always in the first octant then normalize it.
3186 !Bravais lattices are invariant under inversion of any of the basis vectors.
3187  my_qpt = ABS(qpt)/normv(qpt,gmet,"G")
3188 
3189 !Check whether the q is along one of the reduced directions.
3190  idir=0; if (COUNT(ABS(qpt)<tol16) == 2) idir = imax_loc(ABS(qpt))
3191 
3192  if (idir/=0) then ! easy as q is along vector idir.
3193    box_len =  normv(gprimd(:,idir),gmet,"G")
3194    RETURN
3195 
3196  else
3197    !iplane/=0 means that q is placed on one the planes defined by two reciprocal lattice vectors.
3198    iplane=0; if (COUNT(ABS(qpt)<tol16) == 1) iplane = imin_loc(ABS(qpt))
3199    q0box = (/-1,-1,-1/)
3200 
3201    if (iplane/=1) then
3202      x1=one
3203      x2=my_qpt(2)/my_qpt(1)
3204      x3=my_qpt(3)/my_qpt(1)
3205      if (x2<=one+tol16 .and. x3<=one+tol16) q0box=(/x1,x2,x3/)
3206    end if
3207    if (iplane/=2) then
3208      x1=my_qpt(1)/my_qpt(2)
3209      x2=one
3210      x3=my_qpt(3)/my_qpt(2)
3211      if (x1<=one+tol16 .and. x3<=one+tol16) q0box=(/x1,x2,x3/)
3212    end if
3213    if (iplane/=3) then
3214      x1=my_qpt(1)/my_qpt(3)
3215      x2=my_qpt(2)/my_qpt(3)
3216      x3=one
3217      if (x1<=one+tol16 .and. x2<=one+tol16) q0box=(/x1,x2,x3/)
3218    end if
3219 
3220    if (ALL(q0box==(/-1,-1,-1/) )) then
3221      MSG_BUG("Cannot found q0box")
3222    end if
3223 
3224    box_len = normv(q0box,gmet,"G")
3225    RETURN
3226  end if
3227 
3228 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)

PARENTS

CHILDREN

SOURCE

1400 pure logical function bz_mesh_isirred(Kmesh,ik_bz)
1401 
1402 
1403 !This section has been created automatically by the script Abilint (TD).
1404 !Do not modify the following lines by hand.
1405 #undef ABI_FUNC
1406 #define ABI_FUNC 'bz_mesh_isirred'
1407 !End of the abilint section
1408 
1409  implicit none
1410 
1411 !Arguments ------------------------------------
1412 !scalars
1413  integer,intent(in) :: ik_bz
1414  type(kmesh_t),intent(in) :: Kmesh
1415 !arrays
1416 
1417 !Local variables-------------------------------
1418 !scalars
1419  integer :: isym,itim
1420 
1421 ! *********************************************************************
1422 
1423  isym = Kmesh%tabo(ik_bz)
1424  itim = (3-Kmesh%tabi(ik_bz))/2
1425 
1426  ! Be careful here as we assume a particular ordering of symmetries.
1427  bz_mesh_isirred = ( isym==1.and.itim==1.and.ALL(Kmesh%umklp(:,ik_bz)==(/0,0,0/)) )
1428 
1429 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.

PARENTS

      gwls_hamiltonian,mrgscr,setup_bse,setup_bse_interp,setup_screening
      setup_sigma

CHILDREN

SOURCE

2142 subroutine find_qmesh(Qmesh,Cryst,Kmesh)
2143 
2144 
2145 !This section has been created automatically by the script Abilint (TD).
2146 !Do not modify the following lines by hand.
2147 #undef ABI_FUNC
2148 #define ABI_FUNC 'find_qmesh'
2149 !End of the abilint section
2150 
2151  implicit none
2152 
2153 !Arguments ------------------------------------
2154 !scalars
2155  type(kmesh_t),intent(in) :: Kmesh
2156  type(kmesh_t),intent(inout) :: Qmesh
2157  type(crystal_t),intent(in) :: Cryst
2158 
2159 !Local variables-------------------------------
2160 !scalars
2161  integer :: nqibz,kptopt
2162 !arrays
2163  real(dp),allocatable :: qibz(:,:)
2164 
2165 ! *************************************************************************
2166 
2167  !@kmesh_t
2168  ! * Find the number of q-points such that q = k1-k2.
2169  call findnq(Kmesh%nbz,Kmesh%bz,Cryst%nsym,Cryst%symrec,Cryst%symafm,nqibz,Cryst%timrev)
2170  !
2171  ! * Find the coordinates of the q-points in the IBZ.
2172  ABI_MALLOC(qibz,(3,nqibz))
2173  call findq(Kmesh%nbz,Kmesh%bz,Cryst%nsym,Cryst%symrec,Cryst%symafm,Cryst%gprimd,nqibz,qibz,Cryst%timrev)
2174  !
2175  ! * Create the Qmesh object starting from the IBZ.
2176  kptopt = Kmesh%kptopt
2177  call kmesh_init(Qmesh,Cryst,nqibz,qibz,kptopt)
2178 
2179  ABI_FREE(qibz)
2180 
2181 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

PARENTS

      m_bz_mesh

CHILDREN

SOURCE

2212 subroutine findnq(nkbz,kbz,nsym,symrec,symafm,nqibz,timrev)
2213 
2214 
2215 !This section has been created automatically by the script Abilint (TD).
2216 !Do not modify the following lines by hand.
2217 #undef ABI_FUNC
2218 #define ABI_FUNC 'findnq'
2219 !End of the abilint section
2220 
2221  implicit none
2222 
2223 !Arguments ------------------------------------
2224 !scalars
2225  integer,intent(in) :: timrev,nkbz,nsym
2226  integer,intent(out) :: nqibz
2227 !arrays
2228  integer,intent(in) :: symrec(3,3,nsym),symafm(nsym)
2229  real(dp),intent(in) :: kbz(3,nkbz)
2230 
2231 !Local variables ------------------------------
2232 !scalars
2233  integer :: ifound,ik,isym,iq,memory_exhausted,nqall,nqallm,itim,ierr
2234 !arrays
2235  integer :: g0(3)
2236  real(dp) :: qposs(3),qrot(3)
2237  real(dp),allocatable :: qall(:,:)
2238 !************************************************************************
2239 
2240  ! Infinite do-loop to be able to allocate sufficient memory
2241  nqallm=1000
2242  do
2243    memory_exhausted=0
2244    ABI_STAT_MALLOC(qall,(3,nqallm), ierr)
2245    ABI_CHECK(ierr==0, 'out-of-memory qall')
2246    nqall=0
2247 
2248    ! Loop over all k-points in BZ, forming k-k1.
2249    do ik=1,nkbz
2250      qposs(:)=kbz(:,ik)-kbz(:,1)
2251 
2252      ! Check whether this q (or its equivalent) has already been found within a reciprocal lattice vector.
2253      ! Use spatial inversion instead of time reversal whenever possible.
2254      ifound=0
2255      do iq=1,nqall
2256        do itim=1,timrev
2257          do isym=1,nsym
2258            if (symafm(isym)==-1) CYCLE
2259            qrot = (3-2*itim) * MATMUL(symrec(:,:,isym),qall(:,iq))
2260            if (isamek(qrot,qposs,g0)) ifound=ifound+1
2261          end do
2262        end do
2263      end do
2264 
2265      if (ifound==0) then
2266        nqall=nqall+1
2267 
2268        ! If not yet found, check that the allocation is big enough.
2269        if (nqall>nqallm) then
2270          memory_exhausted=1
2271          ABI_FREE(qall)
2272          nqallm=nqallm*2; EXIT ! Exit the do ik=1 loop
2273        end if
2274 
2275        !  Add it to the list.
2276        qall(:,nqall)=qposs(:)
2277      end if
2278    end do
2279 
2280    if (memory_exhausted==0) EXIT
2281  end do !infinite loop
2282 
2283  ABI_FREE(qall)
2284  nqibz=nqall
2285 
2286 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

PARENTS

      m_bz_mesh

CHILDREN

SOURCE

2319 subroutine findq(nkbz,kbz,nsym,symrec,symafm,gprimd,nqibz,qibz,timrev)
2320 
2321 
2322 !This section has been created automatically by the script Abilint (TD).
2323 !Do not modify the following lines by hand.
2324 #undef ABI_FUNC
2325 #define ABI_FUNC 'findq'
2326 !End of the abilint section
2327 
2328  implicit none
2329 
2330 !Arguments ------------------------------------
2331 !scalars
2332  integer,intent(in) :: nkbz,nqibz,nsym,timrev
2333 !arrays
2334  integer,intent(in) :: symrec(3,3,nsym),symafm(nsym)
2335  real(dp),intent(in) :: gprimd(3,3),kbz(3,nkbz)
2336  real(dp),intent(inout) :: qibz(3,nqibz) !vz_i num m_numeric_tools
2337 
2338 !Local variables ------------------------------
2339 !scalars
2340  integer :: ii,ik,iq,iqp,isym,itim
2341  real(dp) :: shift1,qred
2342  logical :: found
2343  character(len=500) :: msg
2344 !arrays
2345  integer :: g0(3)
2346  real(dp) :: gmet(3,3),qposs(3),qrot(3)
2347 
2348 !************************************************************************
2349 
2350  ! Compute reciprocal space metrics
2351  do ii=1,3
2352    gmet(ii,:)=gprimd(1,ii)*gprimd(1,:)+&
2353 &             gprimd(2,ii)*gprimd(2,:)+&
2354 &             gprimd(3,ii)*gprimd(3,:)
2355  end do
2356  !
2357  ! === Loop over k-points in BZ, form k-k1 and translate in first BZ ===
2358  ! iq is the no. of q-points found, zero at the beginning
2359  iq=0
2360  do ik=1,nkbz
2361    qposs(:)=kbz(:,ik)-kbz(:,1)
2362    ! === Check whether this q (or its equivalent) has already been found ===
2363    ! * Use spatial inversion instead of time reversal whenever possible.
2364    found=.FALSE.
2365    do iqp=1,iq
2366      do itim=1,timrev
2367        do isym=1,nsym
2368         if (symafm(isym)==-1) CYCLE
2369         qrot = (3-2*itim) * MATMUL(symrec(:,:,isym),qibz(:,iqp))
2370         if (isamek(qrot,qposs,g0)) found=.TRUE.
2371        end do
2372      end do
2373    end do
2374    if (.not.found) then
2375      iq=iq+1
2376      if (iq>nqibz) then
2377        MSG_BUG(sjoin('iq > nqibz= ',itoa(nqibz)))
2378      end if
2379      qibz(:,iq)=qposs(:)
2380    end if
2381  end do
2382 
2383  if (iq/=nqibz) then
2384    write(msg,'(2(a,i5))')' iq= ',iq,'/= nqibz= ',nqibz
2385    MSG_BUG(msg)
2386  end if
2387  !
2388  ! * Translate q-points to 1st BZ in the interval [-1/2,1/2[
2389  do iq=1,nqibz
2390    do ii=1,3
2391      call wrap2_pmhalf(qibz(ii,iq),qred,shift1)
2392      qibz(ii,iq)=qred
2393    end do
2394  end do
2395 
2396 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
  qbz(3,nqbz)=coordinates of q points
  mG0(3)= For each reduced direction gives the max G0 component to account for umklapp processes

OUTPUT

  iq=index of q vector in BZ table
  g0(3)=reciprocal space vector, to be used in igfft

PARENTS

      calc_sigc_me,calc_sigx_me,cohsex_me,exc_build_block,m_gkk,m_phpi
      m_sigma,prep_calc_ucrpa

CHILDREN

SOURCE

2427 subroutine findqg0(iq,g0,kmkp,nqbz,qbz,mG0)
2428 
2429 
2430 !This section has been created automatically by the script Abilint (TD).
2431 !Do not modify the following lines by hand.
2432 #undef ABI_FUNC
2433 #define ABI_FUNC 'findqg0'
2434 !End of the abilint section
2435 
2436  implicit none
2437 
2438 !Arguments ------------------------------------
2439 !scalars
2440  integer,intent(in) :: nqbz
2441  integer,intent(out) :: iq
2442 !arrays
2443  integer,intent(in) :: mG0(3)
2444  integer,intent(out) :: g0(3)
2445  real(dp),intent(in) :: kmkp(3),qbz(3,nqbz)
2446 
2447 !Local variables-------------------------------
2448 !FIXME if I use 1.0d-4 the jobs crash, should understand why
2449 !scalars
2450  integer :: ig,iqbz,jg01,jg02,jg03
2451  real(dp) :: tolq0=1.0D-3
2452  character(len=500) :: msg
2453 !arrays
2454  real(dp) :: glist1(2*ABS(mG0(1))+1),glist2(2*ABS(mG0(2))+1),glist3(2*ABS(mG0(3))+1)
2455  real(dp) :: qpg0(3),rg(3)
2456 
2457 ! *************************************************************************
2458 
2459  iq=0
2460 
2461  if (ALL(ABS(kmkp(:))<EPSILON(one))) then ! Find q close to 0 ===
2462    do iqbz=1,nqbz
2463      if (ALL(ABS(qbz(:,iqbz))<tolq0)) then
2464        iq=iqbz
2465      end if
2466    end do
2467 
2468    if (iq==0) then
2469      MSG_BUG('Wrong list of q-points: q=0 not present.')
2470    end if
2471    g0(:)=0; RETURN
2472 
2473  else ! q is not zero, find q such as k-kp=q+G0.
2474 
2475    ! Try with G0=0 first.
2476    !do iqbz=1,nqbz
2477    !  if (ALL(ABS(qbz(:,iqbz)-kmkp)<TOL_KDIFF)) then
2478    !    iq=iqbz
2479    !    g0(:)=0; RETURN
2480    !  end if
2481    !end do
2482 
2483 #if 1
2484    glist1(1) = 0
2485    ig = 2
2486    do jg01=1,mG0(1)
2487      glist1(ig)   =  jg01
2488      glist1(ig+1) = -jg01
2489      ig = ig + 2
2490    end do
2491 
2492    glist2(1) = 0
2493    ig = 2
2494    do jg02=1,mG0(2)
2495      glist2(ig)   =  jg02
2496      glist2(ig+1) = -jg02
2497      ig = ig + 2
2498    end do
2499 
2500    glist3(1) = 0
2501    ig = 2
2502    do jg03=1,mG0(3)
2503      glist3(ig)   =  jg03
2504      glist3(ig+1) = -jg03
2505      ig = ig + 2
2506    end do
2507 
2508   g1loop: do jg01=1,2*mG0(1)+1
2509     rg(1) = glist1(jg01)
2510     do jg02=1,2*mG0(2)+1
2511       rg(2) = glist2(jg02)
2512       do jg03=1,2*mG0(3)+1
2513          rg(3) = glist3(jg03)
2514          !
2515          ! * Form q+G0 and check if it is the one.
2516          do iqbz=1,nqbz
2517           qpg0= qbz(:,iqbz) + rg
2518           if (ALL(ABS(qpg0-kmkp)<TOL_KDIFF)) then
2519             iq = iqbz
2520             g0 = NINT(rg)
2521             EXIT g1loop
2522           end if
2523         end do
2524 
2525       end do
2526     end do
2527   end do g1loop
2528 
2529 #else
2530    g1loop: do jg01=-mG0(1),mG0(1)
2531      rg(1) = DBLE(jg01)
2532      do jg02=-mG0(2),mG0(2)
2533        rg(2) = DBLE(jg02)
2534        do jg03=-mG0(3),mG0(3)
2535           rg(3) = DBLE(jg03)
2536           !
2537           ! * Form q+G0 and check if it is the one.
2538           do iqbz=1,nqbz
2539            qpg0= qbz(:,iqbz) + rg
2540            if (ALL(ABS(qpg0-kmkp)<TOL_KDIFF)) then
2541              iq = iqbz
2542              g0 = (/jg01,jg02,jg03/)
2543              EXIT g1loop
2544            end if
2545          end do
2546 
2547        end do
2548      end do
2549    end do g1loop
2550 #endif
2551 
2552    if (iq==0) then
2553      write(msg,'(a,3f9.5)')' q = k-kp+G0 not found. kmkp = ',kmkp
2554      MSG_ERROR(msg)
2555    end if
2556 
2557  end if
2558 
2559 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)

PARENTS

      cchi0

CHILDREN

SOURCE

1072 subroutine get_BZ_diff(Kmesh,k1,k2,idiff_bz,g0,nfound)
1073 
1074 
1075 !This section has been created automatically by the script Abilint (TD).
1076 !Do not modify the following lines by hand.
1077 #undef ABI_FUNC
1078 #define ABI_FUNC 'get_BZ_diff'
1079 !End of the abilint section
1080 
1081  implicit none
1082 
1083 !Arguments ------------------------------------
1084 !scalars
1085  integer,intent(out) :: idiff_bz,nfound
1086  type(kmesh_t),intent(in) :: Kmesh
1087 !arrays
1088  integer,intent(out) :: g0(3)
1089  real(dp),intent(in) :: k1(3),k2(3)
1090 
1091 !Local variables-------------------------------
1092 !scalars
1093  integer :: ikp
1094  character(len=500) :: msg
1095 !arrays
1096  integer :: umklp(3)
1097  real(dp) :: kdiff(3),ktrial(3)
1098 
1099 ! *********************************************************************
1100 
1101  if (.not.has_BZ_item(Kmesh,k1,ikp,umklp)) then
1102    write(msg,'(a,3f12.6)')' first point must be in BZ: ',k1
1103    MSG_ERROR(msg)
1104  end if
1105 
1106  kdiff   = k1-k2
1107  nfound  = 0
1108  idiff_bz= 0
1109  !
1110  ! === Find p such k1-k2=p+g0 where p in the BZ ===
1111  do ikp=1,Kmesh%nbz
1112    ktrial=Kmesh%bz(:,ikp)
1113    if (isamek(kdiff,ktrial,umklp)) then
1114      idiff_bz=ikp
1115      g0=umklp
1116      nfound=nfound+1
1117    end if
1118  end do
1119  !
1120  ! === Check if p has not found of found more than once ===
1121  ! * For extremely dense meshes, tol1q in defs_basis might be too large!
1122  if (nfound/=1) then
1123    if (nfound==0) then
1124      MSG_WARNING(" k1-k2-G0 not found in BZ")
1125    else
1126      MSG_WARNING(sjoin(' Multiple k1-k2-G0 found in BZ, nfound= ', itoa(nfound)))
1127    end if
1128    write(msg,'(4a,3(a,3f12.6,a))')&
1129 &    ' k1    = ',k1   ,ch10,&
1130 &    ' k2    = ',k2   ,ch10,&
1131 &    ' k1-k2 = ',kdiff,ch10
1132    MSG_WARNING(msg)
1133  end if
1134 
1135 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.

PARENTS

      calc_optical_mels,calc_sigc_me,calc_sigx_me,calc_ucrpa,cchi0,cchi0q0
      cchi0q0_intraband,cohsex_me,debug_tools,exc_build_block,exc_build_ham
      m_dyson_solver,m_ppmodel,m_screen,m_screening,m_shirley,m_vcoul,m_wfd
      paw_symcprj,prep_calc_ucrpa,random_stopping_power,read_plowannier
      setup_bse,setup_screening,setup_sigma

CHILDREN

SOURCE

941 subroutine get_bz_item(Kmesh,ik_bz,kbz,ik_ibz,isym,itim,ph_mkbzt,umklp,isirred)
942 
943 
944 !This section has been created automatically by the script Abilint (TD).
945 !Do not modify the following lines by hand.
946 #undef ABI_FUNC
947 #define ABI_FUNC 'get_bz_item'
948 !End of the abilint section
949 
950  implicit none
951 
952 !Arguments ------------------------------------
953 !scalars
954  integer,intent(in) :: ik_bz
955  integer,intent(out) :: ik_ibz,isym,itim
956  complex(dpc),optional,intent(out) :: ph_mkbzt
957  logical,optional,intent(out) :: isirred
958  type(kmesh_t),intent(in) :: Kmesh
959 !arrays
960  integer,optional,intent(out) :: umklp(3)
961  real(dp),intent(out) :: kbz(3)
962 
963 !Local variables-------------------------------
964 !scalars
965  character(len=500) :: msg
966 
967 ! *********************************************************************
968 
969  if (ik_bz>Kmesh%nbz.or.ik_bz<=0) then
970    write(msg,'(a,2i3)')' Wrong value for ik_bz: ',ik_bz,Kmesh%nbz
971    MSG_BUG(msg)
972  end if
973 
974  kbz    = Kmesh%bz(:,ik_bz)
975  ik_ibz = Kmesh%tab(ik_bz)
976  isym   = Kmesh%tabo(ik_bz)
977  itim   = (3-Kmesh%tabi(ik_bz))/2
978 
979  if (PRESENT(ph_mkbzt)) ph_mkbzt=Kmesh%tabp(ik_bz)
980  if (PRESENT(umklp))    umklp   =Kmesh%umklp(:,ik_bz)
981  ! Be careful here as we assume a particular ordering of symmetries.
982  if (PRESENT(isirred))  isirred = (isym==1.and.itim==1.and.ALL(Kmesh%umklp(:,ik_bz)==(/0,0,0/)))
983 
984 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

PARENTS

      paw_symcprj

CHILDREN

SOURCE

1014 subroutine get_IBZ_item(Kmesh,ik_ibz,kibz,wtk)
1015 
1016 
1017 !This section has been created automatically by the script Abilint (TD).
1018 !Do not modify the following lines by hand.
1019 #undef ABI_FUNC
1020 #define ABI_FUNC 'get_IBZ_item'
1021 !End of the abilint section
1022 
1023  implicit none
1024 
1025 !Arguments ------------------------------------
1026 !scalars
1027  integer,intent(in) :: ik_ibz
1028  real(dp),intent(out) :: wtk
1029  type(kmesh_t),intent(in) :: Kmesh
1030 !arrays
1031  real(dp),intent(out) :: kibz(3)
1032 
1033 ! *********************************************************************
1034 
1035  if (ik_ibz>Kmesh%nibz.or.ik_ibz<=0) then
1036    MSG_BUG(sjoin('wrong value for ik_ibz: ',itoa(ik_ibz)))
1037  end if
1038 
1039  kibz=Kmesh%ibz(:,ik_ibz)
1040  wtk =Kmesh%wt(ik_ibz)
1041 
1042 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.

PARENTS

      setup_bse,setup_screening,setup_sigma

CHILDREN

SOURCE

1814 subroutine get_ng0sh(nk1,kbz1,nk2,kbz2,nkfold,kfold,tolq0,opt_ng0)
1815 
1816 
1817 !This section has been created automatically by the script Abilint (TD).
1818 !Do not modify the following lines by hand.
1819 #undef ABI_FUNC
1820 #define ABI_FUNC 'get_ng0sh'
1821 !End of the abilint section
1822 
1823  implicit none
1824 
1825 !Arguments ------------------------------------
1826 !scalars
1827  integer,intent(in) :: nk1,nk2,nkfold
1828  real(dp),intent(in) :: tolq0
1829 !arrays
1830  integer,intent(out) :: opt_ng0(3)
1831  real(dp),intent(in) :: kbz1(3,nk1),kbz2(3,nk2),kfold(3,nkfold)
1832 
1833 !Local variables-------------------------------
1834 !scalars
1835  integer :: i1,i2,ikf,ind,factor
1836  real(dp) :: normdiff,tempnorm,smallestlen
1837  logical :: found
1838  character(len=500) :: msg
1839 !arrays
1840  integer :: roundk(3),kbigdiff(3),kbigfold(3,nkfold),iperm(nkfold)
1841  real(dp) :: k1mk2(3),ksmalldiff(3),norm(nkfold),ksmallfold(3,nkfold)
1842 
1843 !************************************************************************
1844 
1845  ! Compute smallest length of one component
1846  ! To get a sufficiently large factor to order vectors
1847  smallestlen = one
1848 
1849  ! Compute integer part and fractional part, [0,1[ of kfold
1850  do ikf = 1,nkfold
1851    kbigfold(:,ikf) = FLOOR(kfold(:,ikf)+tol7)
1852    ksmallfold(:,ikf) = kfold(:,ikf)-kbigfold(:,ikf)
1853 
1854    if (ABS(ksmallfold(1,ikf)) > tol7) &
1855 &     smallestlen = MIN(smallestlen, ABS(ksmallfold(1,ikf)))
1856    if (ABS(ksmallfold(2,ikf)) > tol7) &
1857 &     smallestlen = MIN(smallestlen, ABS(ksmallfold(2,ikf)))
1858    if (ABS(ksmallfold(3,ikf)) > tol7) &
1859 &     smallestlen = MIN(smallestlen, ABS(ksmallfold(3,ikf)))
1860 
1861  end do
1862 
1863  ! WARNING ! This could not be sufficient if tested k1 - k2 has lower
1864  ! components than smallestlen. The factor 10 is giving us a security margin.
1865  factor = 10*(int(one/smallestlen)+1)
1866 
1867  ! Loop again over kfold vectors to give each term its norm
1868  do ikf=1, nkfold
1869    iperm(ikf) = ikf
1870 
1871    ! Computing a sort of norm with order of components (used for ordering)
1872    call getkptnorm_bycomponent(ksmallfold(:,ikf),factor,tempnorm)
1873 
1874    norm(ikf) = tempnorm
1875  end do
1876 
1877  ! Sorting list of kfold vectors
1878  call sort_dp(nkfold,norm,iperm,tol14)
1879 
1880  ! Loop over all k1 - k2
1881  opt_ng0(:)=0
1882  do i2=1,nk2
1883    ! This is used in case of screening calculation.
1884    ! If q is small treat it as zero. In this case, indeed,
1885    ! we use q=0 to calculate the oscillator matrix elements.
1886    if (is_zero(kbz2(:,i2),tolq0)) CYCLE
1887    do i1=1,nk1
1888      found=.FALSE.
1889 
1890      ! Separating in integer part and fractionary part
1891      k1mk2(:) = kbz1(:,i1)-kbz2(:,i2)
1892      ! Adding small tol to prevent 1 being in fractionary part
1893      kbigdiff(:) = FLOOR(k1mk2(:)+tol7)
1894      ksmalldiff(:) = k1mk2(:)-kbigdiff(:)
1895 
1896      call getkptnorm_bycomponent(ksmalldiff(:),factor,normdiff)
1897 
1898      ! Try to find the right smallkfold, corresponding to ksmalldiff
1899      ind = bisect(norm,normdiff)
1900 
1901      if (ind > 0) then
1902        if(ABS(norm(ind) - normdiff) < TOL_KDIFF) then
1903           found = .TRUE.
1904        end if
1905      end if
1906      if(ind < nkfold) then
1907        if(ABS(norm(ind+1) - normdiff) < TOL_KDIFF) then
1908           found = .TRUE.
1909           ind = ind + 1
1910        end if
1911      end if
1912 
1913      if (.not. found) then
1914        write(msg,'(a,2(2a,i4,3es16.8),a)')&
1915 &        'Not able to found umklapp G0 vector such as k1-k2 = kf+G0',ch10,&
1916 &        'point1 = ',i1,kbz1(:,i1),ch10,&
1917 &        'point2 = ',i2,kbz2(:,i2),ch10
1918        MSG_ERROR(msg)
1919      else
1920        ! We have found one k, extracting the max g0
1921        roundk(:) = ABS(kbigdiff(:) - kbigfold(:,iperm(ind)))
1922        opt_ng0(1) = MAX(opt_ng0(1),roundk(1))
1923        opt_ng0(2) = MAX(opt_ng0(2),roundk(2))
1924        opt_ng0(3) = MAX(opt_ng0(3),roundk(3))
1925      end if
1926    end do
1927  end do
1928 
1929 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

PARENTS

      m_bz_mesh

CHILDREN

SOURCE

1956 subroutine getkptnorm_bycomponent(vect,factor,norm)
1957 
1958 !Arguments ------------------------------------
1959 !scalars
1960 
1961 !This section has been created automatically by the script Abilint (TD).
1962 !Do not modify the following lines by hand.
1963 #undef ABI_FUNC
1964 #define ABI_FUNC 'getkptnorm_bycomponent'
1965 !End of the abilint section
1966 
1967  integer,intent(in) :: factor
1968  real(dp),intent(out):: norm
1969 !arrays
1970  real(dp),intent(in) :: vect(3)
1971 !Local variables-------------------------------
1972 !scalars
1973  character(len=500) :: msg
1974 
1975 ! *************************************************************************
1976 
1977  ! Checking the factor is large enough
1978  !(skipping zero components, since in this case the product will be 0)
1979  if(ANY(vect(:)*factor < 1.0 .and. vect(:) > tol7)) then
1980     write(msg,'(a,a,a,a,a,a,a,a)') ' Not able to give unique norm to order vectors',ch10,&
1981 &       'This is likely related to a truncation error for a k-point in the input file',ch10,&
1982 &       'Always prefer fractional numbers in the input file instead of truncated ones',ch10,&
1983 &       '(e.g. 1/6 instead of 0.166666667)',ch10
1984     MSG_ERROR(msg)
1985  end if
1986 
1987  norm = (vect(1)*factor+vect(2))*factor+vect(3)
1988 
1989 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.

PARENTS

 FIXME
  Switch to routine version. Due to side-effects the present implementation
  might be source of bugs in logical statements

CHILDREN

SOURCE

1264 logical function has_BZ_item(Kmesh,item,ikbz,g0)
1265 
1266 
1267 !This section has been created automatically by the script Abilint (TD).
1268 !Do not modify the following lines by hand.
1269 #undef ABI_FUNC
1270 #define ABI_FUNC 'has_BZ_item'
1271 !End of the abilint section
1272 
1273  implicit none
1274 
1275 !Arguments ------------------------------------
1276 !scalars
1277  integer,intent(out) :: ikbz
1278  type(kmesh_t),intent(in) :: Kmesh
1279 !arrays
1280  integer,intent(out) :: g0(3)
1281  real(dp),intent(in) :: item(3)
1282 
1283 !Local variables-------------------------------
1284 !scalars
1285  integer :: ik_bz,yetfound
1286 !arrays
1287  integer :: g0_tmp(3)
1288 
1289 ! *************************************************************************
1290 
1291  has_BZ_item=.FALSE.; ikbz=0; g0=0; yetfound=0
1292  do ik_bz=1,Kmesh%nbz
1293    if (isamek(item,Kmesh%bz(:,ik_bz),g0_tmp)) then
1294      has_BZ_item=.TRUE.
1295      ikbz=ik_bz
1296      g0 = g0_tmp
1297      yetfound=yetfound+1
1298      !EXIT
1299    end if
1300  end do
1301 
1302  if (yetfound/=0.and.yetfound/=1) then
1303    MSG_ERROR('Multiple k-points found')
1304  end if
1305 
1306 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.

PARENTS

CHILDREN

SOURCE

1333 logical function has_IBZ_item(Kmesh,item,ikibz,g0)
1334 
1335 
1336 !This section has been created automatically by the script Abilint (TD).
1337 !Do not modify the following lines by hand.
1338 #undef ABI_FUNC
1339 #define ABI_FUNC 'has_IBZ_item'
1340 !End of the abilint section
1341 
1342  implicit none
1343 
1344 !Arguments ------------------------------------
1345 !scalars
1346  integer,intent(out) :: ikibz
1347  type(kmesh_t),intent(in) :: Kmesh
1348 !arrays
1349  integer,intent(out) :: g0(3)
1350  real(dp),intent(in) :: item(3)
1351 
1352 !Local variables-------------------------------
1353 !scalars
1354  integer :: ik_ibz,yetfound
1355  !character(len=500) :: msg
1356 !arrays
1357  integer :: g0_tmp(3)
1358 
1359 ! *************************************************************************
1360 
1361  has_IBZ_item=.FALSE.; ikibz=0; g0=0; yetfound=0
1362  do ik_ibz=1,Kmesh%nibz
1363    if (isamek(item,Kmesh%ibz(:,ik_ibz),g0_tmp)) then
1364      has_IBZ_item=.TRUE.
1365      ikibz=ik_ibz
1366      g0 = g0_tmp
1367      yetfound=yetfound+1
1368      !EXIT
1369    end if
1370  end do
1371 
1372  if (yetfound/=0.and.yetfound/=1) then
1373    MSG_BUG("multiple k-points found")
1374  end if
1375 
1376 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.

PARENTS

      m_bz_mesh

CHILDREN

SOURCE

1615 subroutine identk(kibz,nkibz,nkbzmx,nsym,timrev,symrec,symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk,ref_bz)
1616 
1617 
1618 !This section has been created automatically by the script Abilint (TD).
1619 !Do not modify the following lines by hand.
1620 #undef ABI_FUNC
1621 #define ABI_FUNC 'identk'
1622 !End of the abilint section
1623 
1624  implicit none
1625 
1626 !Arguments ------------------------------------
1627 !scalars
1628  integer,intent(in) :: nkbzmx,nkibz,nsym,timrev
1629  integer,intent(out) :: nkbz
1630 !arrays
1631  integer,intent(in) :: symafm(nsym),symrec(3,3,nsym)
1632  integer,intent(out) :: ktab(nkbzmx),ktabi(nkbzmx),ktabo(nkbzmx)
1633  real(dp),intent(in) :: kibz(3,nkibz)
1634  real(dp),intent(out) :: kbz(3,nkbzmx),wtk(nkibz)
1635  real(dp),optional,intent(in) :: ref_bz(:,:)
1636 
1637 !Local variables ------------------------------
1638 !scalars
1639  integer :: ik1,ik2,ikbz,ikibz,iold,isym,itim
1640  integer :: ikref,nkref,isym_swp,itim_swp,ikibz_swp
1641  logical :: is_irred_set
1642  logical :: found,ltest
1643  character(len=500) :: msg
1644 !arrays
1645  integer :: g0(3)
1646  real(dp) :: knew(3),k1(3),k2(3)
1647  real(dp) :: kref(3),kbz_swp(3)
1648 
1649 ! *************************************************************************
1650 
1651  DBG_ENTER("COLL")
1652  !
1653  ! === Check whether kibz really forms an irreducible set ===
1654  is_irred_set=.TRUE.
1655  do ik1=1,nkibz-1
1656    k1=kibz(:,ik1)
1657    do ik2=ik1+1,nkibz
1658      k2=kibz(:,ik2)
1659 
1660      do itim=1,timrev
1661        do isym=1,nsym
1662          if (symafm(isym)==-1) CYCLE
1663          knew = (3-2*itim) * MATMUL(symrec(:,:,isym),k2)
1664          if (isamek(k1,knew,g0)) then
1665            is_irred_set=.FALSE.
1666            write(msg,'(2(a,3f8.4),2(a,i2))')&
1667 &            ' k1 = ',k1,' is symmetrical of k2 = ',k2,' through sym = ',isym,' itim = ',itim
1668            MSG_WARNING(msg)
1669          end if
1670        end do
1671      end do
1672 
1673    end do
1674  end do
1675 
1676  !call klist_isirred(nkibz,kibz,Cryst,nimg)
1677 
1678  if (.not.is_irred_set) then
1679    msg = "Input array kibz does not constitute an irreducible set."
1680    MSG_WARNING(msg)
1681  end if
1682  !
1683  ! === Loop over k-points in IBZ ===
1684  ! * Start with zero no. of k-points found.
1685  nkbz=0
1686  do ikibz=1,nkibz
1687    wtk(ikibz) = zero
1688 
1689    ! === Loop over time-reversal I and symmetry operations S  ===
1690    ! * Use spatial inversion instead of time reversal whenever possible.
1691    do itim=1,timrev
1692      do isym=1,nsym
1693        if (symafm(isym)==-1) CYCLE
1694        !
1695        ! * Form IS k
1696        knew=(3-2*itim)*MATMUL(symrec(:,:,isym),kibz(:,ikibz))
1697        !
1698        ! * Check whether it has already been found (to within a RL vector).
1699        iold=0
1700        do ikbz=1,nkbz
1701          if (isamek(knew,kbz(:,ikbz),g0)) then
1702            iold=iold+1
1703            exit
1704          end if
1705        end do
1706        !
1707        ! * If not yet found add to kbz and increase the weight.
1708        if (iold==0) then
1709          nkbz=nkbz+1
1710          wtk(ikibz)=wtk(ikibz)+one
1711          if (nkbz>nkbzmx) then
1712            MSG_BUG(sjoin('nkbzmx too small, nkbzmx = ',itoa(nkbzmx),', increase nkbzmx !'))
1713          end if
1714          kbz(:,nkbz) = knew(:)
1715          ktab (nkbz) = ikibz
1716          ktabo(nkbz) = isym
1717          ktabi(nkbz) = 3-2*itim
1718        end if
1719        !
1720       end do
1721    end do
1722 
1723  end do !ikibz
1724 
1725  if (PRESENT(ref_bz)) then
1726    call wrtout(std_out," Pruning the k-points not in ref_bz then reordering tables","COLL")
1727 
1728    nkref = SIZE(ref_bz,DIM=2)
1729    ltest = (nkref>=nkbz.and.nkref<=nkbzmx)
1730    if (.not.ltest) then
1731      write(msg,'(3(a,i0))')" Wrong value for nkref: nkref= ",nkref," nkbz= ",nkbz," nkbzmx =",nkbzmx
1732      MSG_WARNING(msg)
1733    end if
1734 
1735    do ikref=1,nkref
1736      kref = ref_bz(:,ikref)
1737      found=.FALSE.
1738 
1739      do ikbz=1,nkbz ! Loop on the set of BZ points found above.
1740        if (isequalk(kref,kbz(:,ikbz))) then ! Swap indeces.
1741          kbz_swp   = kbz(:,ikref)
1742          ikibz_swp = ktab (ikref)
1743          isym_swp  = ktabo(ikref)
1744          itim_swp  = ktabi(ikref)
1745 
1746          kbz(:,ikref) = kref
1747          ktab (ikref) = ktab (ikbz)
1748          ktabo(ikref) = ktabo(ikbz)
1749          ktabi(ikref) = ktabi(ikbz)
1750 
1751          kbz(:,ikbz) = kbz_swp
1752          ktab (ikbz) = ikibz_swp
1753          ktabo(ikbz) = isym_swp
1754          ktabi(ikbz) = itim_swp
1755 
1756          found=.TRUE.; EXIT
1757        end if
1758      end do
1759 
1760      if (.not.found) then
1761        write(msg,'(a,3es16.8)')" One of the k-point in ref_bz is not a symmetrical image of the IBZ: ",kref
1762        MSG_ERROR(msg)
1763      end if
1764    end do
1765    !
1766    ! Change nkbz to nkref, then get the new weights.
1767    nkbz=nkref; wtk=zero
1768    do ikref=1,nkref
1769      ikibz = ktab(ikref)
1770      wtk(ikibz) = wtk(ikibz) + 1
1771    end do
1772  end if ! PRESENT(ref_bz)
1773  !
1774  ! * Weights are normalized to 1.
1775  wtk = wtk/SUM(wtk)
1776 
1777  DBG_EXIT("COLL")
1778 
1779 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

PARENTS

CHILDREN

SOURCE

1162 logical function isamek(k1,k2,g0)
1163 
1164 
1165 !This section has been created automatically by the script Abilint (TD).
1166 !Do not modify the following lines by hand.
1167 #undef ABI_FUNC
1168 #define ABI_FUNC 'isamek'
1169 !End of the abilint section
1170 
1171  implicit none
1172 
1173 !Arguments ------------------------------------
1174 !arrays
1175  integer,intent(out) :: g0(3)
1176  real(dp),intent(in) :: k1(3),k2(3)
1177 
1178 ! *************************************************************************
1179 
1180  isamek = isinteger(k1-k2,TOL_KDIFF)
1181 
1182  if (isamek) then
1183    g0=NINT(k1-k2)
1184  else
1185    g0=HUGE(1)
1186  end if
1187 
1188 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

PARENTS

CHILDREN

SOURCE

1211 logical function isequalk(q1,q2)
1212 
1213 
1214 !This section has been created automatically by the script Abilint (TD).
1215 !Do not modify the following lines by hand.
1216 #undef ABI_FUNC
1217 #define ABI_FUNC 'isequalk'
1218 !End of the abilint section
1219 
1220  implicit none
1221 
1222 !Arguments ------------------------------------
1223  real(dp),intent(in) :: q1(3),q2(3)
1224 
1225 !Local variables-------------------------------
1226 !arrays
1227  integer :: g0(3)
1228 ! *************************************************************************
1229 
1230  isequalk = isamek(q1,q2,g0)
1231 
1232 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.

PARENTS

      bethe_salpeter,cchi0q0_intraband,gwls_hamiltonian,m_shirley,mlwfovlp_qp
      mrgscr,screening,sigma

CHILDREN

SOURCE

567 subroutine kmesh_free(Kmesh)
568 
569 
570 !This section has been created automatically by the script Abilint (TD).
571 !Do not modify the following lines by hand.
572 #undef ABI_FUNC
573 #define ABI_FUNC 'kmesh_free'
574 !End of the abilint section
575 
576  implicit none
577 
578 !Arguments ------------------------------------
579 !scalars
580  type(kmesh_t),intent(inout) :: Kmesh
581 
582 ! *********************************************************************
583 
584  !@kmesh_t
585 !integer
586  if (allocated(Kmesh%rottb)) then
587    ABI_FREE(Kmesh%rottb)
588  end if
589  if (allocated(Kmesh%rottbm1)) then
590    ABI_FREE(Kmesh%rottbm1)
591  end if
592  if (allocated(Kmesh%tab)) then
593    ABI_FREE(Kmesh%tab)
594  end if
595  if (allocated(Kmesh%tabi)) then
596    ABI_FREE(Kmesh%tabi)
597  end if
598  if (allocated(Kmesh%tabo)) then
599    ABI_FREE(Kmesh%tabo)
600  end if
601  if (allocated(Kmesh%umklp)) then
602    ABI_FREE(Kmesh%umklp)
603  end if
604 
605 !real
606  if (allocated(Kmesh%ibz)) then
607    ABI_FREE(Kmesh%ibz)
608  end if
609  if (allocated(Kmesh%bz)) then
610    ABI_FREE(Kmesh%bz)
611  end if
612  if (allocated(Kmesh%shift)) then
613    ABI_FREE(Kmesh%shift)
614  end if
615  if (allocated(Kmesh%wt)) then
616    ABI_FREE(Kmesh%wt)
617  end if
618 
619 !complex
620  if (allocated(Kmesh%tabp)) then
621    ABI_FREE(Kmesh%tabp)
622  end if
623 
624 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.

PARENTS

      cchi0q0_intraband,gwls_hamiltonian,m_bz_mesh,mlwfovlp_qp,mrgscr
      setup_bse,setup_bse_interp,setup_screening,setup_sigma

CHILDREN

SOURCE

394 subroutine kmesh_init(Kmesh,Cryst,nkibz,kibz,kptopt,wrap_1zone,ref_bz,break_symmetry)
395 
396 
397 !This section has been created automatically by the script Abilint (TD).
398 !Do not modify the following lines by hand.
399 #undef ABI_FUNC
400 #define ABI_FUNC 'kmesh_init'
401 !End of the abilint section
402 
403  implicit none
404 
405 !Arguments ------------------------------------
406 !scalars
407  integer,intent(in) :: nkibz,kptopt
408  logical,optional,intent(in) :: wrap_1zone,break_symmetry
409  type(kmesh_t),intent(inout) :: Kmesh
410  type(crystal_t),intent(in) :: Cryst
411 !arrays
412  real(dp),intent(in) :: kibz(3,nkibz)
413  real(dp),optional,intent(in) :: ref_bz(:,:)
414 
415 !Local variables-------------------------------
416 !scalars
417  integer :: ik_bz,ik_ibz,isym,nkbz,nkbzX,nsym,timrev,itim
418  real(dp) :: shift1,shift2,shift3
419  logical :: ltest,do_wrap,do_hack
420 !arrays
421  integer,allocatable :: ktab(:),ktabi(:),ktabo(:)
422  real(dp) :: rm1t(3),kbz_wrap(3)
423  real(dp),allocatable :: kbz(:,:),wtk(:)
424 
425 ! *************************************************************************
426 
427  DBG_ENTER("COLL")
428 
429  !@kmesh_t
430  ! === Initial tests on input arguments ===
431  ltest=(Cryst%timrev==1.or.Cryst%timrev==2)
432  ABI_CHECK(ltest, sjoin('Wrong value for timrev= ', itoa(Cryst%timrev)))
433 
434  if (ALL(kptopt/=(/1,3/))) then
435    MSG_WARNING(sjoin("Not allowed value for kptopt: ", itoa(kptopt)))
436  end if
437 
438  Kmesh%kptopt = kptopt
439 
440  nsym   = Cryst%nsym
441  timrev = Cryst%timrev
442 
443  !
444  ! === Find BZ from IBZ and fill tables ===
445  nkbzX=nkibz*nsym*timrev ! Maximum possible number
446  ABI_MALLOC(kbz,(3,nkbzX))
447  ABI_MALLOC(wtk,(nkibz))
448  ABI_MALLOC(ktab,(nkbzX))
449  ABI_MALLOC(ktabi,(nkbzX))
450  ABI_MALLOC(ktabo,(nkbzX))
451 
452  if (PRESENT(ref_bz)) then
453    call identk(kibz,nkibz,nkbzX,nsym,timrev,cryst%symrec,cryst%symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk,ref_bz=ref_bz)
454  else
455    call identk(kibz,nkibz,nkbzX,nsym,timrev,cryst%symrec,cryst%symafm,kbz,ktab,ktabi,ktabo,nkbz,wtk)
456  end if
457 
458  ! TODO: Force the k-points to be in the first Brillouin zone.
459  !  Now the GW tests seem to be OK, additional tests have to be done though.
460  do_wrap=.FALSE.; if (PRESENT(wrap_1zone)) do_wrap=wrap_1zone
461  !do_wrap=.TRUE.
462 
463  if (do_wrap) then ! Wrap the BZ points in the interval ]-1/2,1/2]
464    do ik_bz=1,nkbz
465      call wrap2_pmhalf(kbz(1,ik_bz),kbz_wrap(1),shift1)
466      call wrap2_pmhalf(kbz(2,ik_bz),kbz_wrap(2),shift2)
467      call wrap2_pmhalf(kbz(3,ik_bz),kbz_wrap(3),shift3)
468      kbz(:,ik_bz) = kbz_wrap
469    end do
470  end if
471  !
472  ! ================================================================
473  ! ==== Create data structure to store information on k-points ====
474  ! ================================================================
475  !
476  ! * Dimensions.
477  Kmesh%nbz   =nkbz      ! Number of points in the full BZ
478  Kmesh%nibz  =nkibz     ! Number of points in the IBZ
479  Kmesh%nsym  =nsym      ! Number of operations
480  Kmesh%timrev=timrev    ! 2 if time-reversal is used, 1 otherwise
481  !
482  ! * Arrays.
483  Kmesh%gmet   = Cryst%gmet
484  Kmesh%gprimd = Cryst%gprimd
485 
486  ABI_MALLOC(Kmesh%bz ,(3,nkbz))
487  Kmesh%bz   =  kbz(:,1:nkbz )  ! Red. coordinates of points in full BZ.
488  ABI_MALLOC(Kmesh%ibz,(3,nkibz))
489  Kmesh%ibz  = kibz(:,1:nkibz)  ! Red. coordinates of points in IBZ.
490 
491  ABI_MALLOC(Kmesh%tab ,(nkbz))
492  Kmesh%tab  = ktab (1:nkbz)    ! Index of the irred. point in the array IBZ.
493  ABI_MALLOC(Kmesh%tabi,(nkbz))
494  Kmesh%tabi = ktabi(1:nkbz)    !-1 if time reversal must be used to obtain this point,
495  ABI_MALLOC(Kmesh%tabo,(nkbz))
496  Kmesh%tabo = ktabo(1:nkbz)    ! Symm. operation that rotates k_IBZ onto \pm k_BZ
497                                                              ! (depending on tabi)
498  ABI_MALLOC(Kmesh%wt,(nkibz))
499  Kmesh%wt(:)= wtk(1:nkibz)     ! Weight for each k_IBZ
500 
501  ABI_MALLOC(Kmesh%rottbm1,(nkbz,timrev,nsym))
502  ABI_MALLOC(Kmesh%rottb  ,(nkbz,timrev,nsym))
503 
504  do_hack = .FALSE.
505  if (PRESENT(ref_bz) .and. PRESENT(break_symmetry)) then
506    do_hack = break_symmetry
507  end if
508 
509  if (do_hack) then
510    MSG_WARNING("Hacking the rottb tables!")
511    do ik_bz=1,nkbz
512      Kmesh%rottbm1(ik_bz,:,:) = ik_bz
513      Kmesh%rottb  (ik_bz,:,:) = ik_bz
514    end do
515  else
516    call setup_k_rotation(nsym,timrev,cryst%symrec,nkbz,Kmesh%bz,Cryst%gmet,Kmesh%rottb,Kmesh%rottbm1)
517  end if
518 
519  ! TODO umklp can be calculated inside setup_k_rotation.
520  ABI_MALLOC(Kmesh%umklp,(3,nkbz))
521  do ik_bz=1,nkbz
522    ik_ibz= Kmesh%tab (ik_bz)
523    isym  = Kmesh%tabo(ik_bz)
524    itim  = (3-Kmesh%tabi(ik_bz))/2
525    Kmesh%umklp(:,ik_bz) = NINT( -Kmesh%bz(:,ik_bz) + (3-2*itim)*MATMUL(cryst%symrec(:,:,isym),Kmesh%ibz(:,ik_ibz)) )
526  end do
527 
528  ABI_MALLOC(Kmesh%tabp,(nkbz))
529  do ik_bz=1,nkbz
530    isym  =Kmesh%tabo(ik_bz)
531    ik_ibz=Kmesh%tab (ik_bz)
532    rm1t=MATMUL(TRANSPOSE(cryst%symrec(:,:,isym)),cryst%tnons(:,isym))
533    Kmesh%tabp(ik_bz)=EXP(-(0.,1.)*two_pi*DOT_PRODUCT(kibz(:,ik_ibz),rm1t))
534  end do
535 
536  ABI_FREE(kbz)
537  ABI_FREE(wtk)
538  ABI_FREE(ktab)
539  ABI_FREE(ktabi)
540  ABI_FREE(ktabo)
541 
542  DBG_EXIT("COLL")
543 
544 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.

PARENTS

      gwls_hamiltonian,mrgscr,setup_bse,setup_bse_interp,setup_screening
      setup_sigma

CHILDREN

SOURCE

654 subroutine kmesh_print(Kmesh,header,unit,prtvol,mode_paral)
655 
656 
657 !This section has been created automatically by the script Abilint (TD).
658 !Do not modify the following lines by hand.
659 #undef ABI_FUNC
660 #define ABI_FUNC 'kmesh_print'
661 !End of the abilint section
662 
663  implicit none
664 
665 !Arguments ------------------------------------
666 !scalars
667  integer,optional,intent(in) :: prtvol,unit
668  character(len=4),optional,intent(in) :: mode_paral
669  character(len=*),optional,intent(in) :: header
670  type(kmesh_t),intent(in) :: Kmesh
671 
672 !Local variables-------------------------------
673 !scalars
674  integer,parameter :: nmaxk=50
675  integer :: ii,ik,my_unt,my_prtvol
676  character(len=100) :: fmt
677  character(len=4) :: my_mode
678  character(len=500) :: msg
679 
680 ! *************************************************************************
681 
682  my_unt =std_out; if (PRESENT(unit      )) my_unt   =unit
683  my_prtvol=0    ; if (PRESENT(prtvol    )) my_prtvol=prtvol
684  my_mode='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
685 
686  msg=' ==== Info on the Kmesh% object ==== '
687  if (PRESENT(header)) msg=' ==== '//TRIM(ADJUSTL(header))//' ==== '
688  call wrtout(my_unt,msg,my_mode)
689 
690  write(msg,'(a,i5,3a)')&
691 &  ' Number of points in the irreducible wedge : ',Kmesh%nibz,ch10,&
692 &  ' Reduced coordinates and weights : ',ch10
693  call wrtout(my_unt,msg,my_mode)
694 
695  write(fmt,*)'(1x,i5,a,2x,3es16.8,3x,f11.5)'
696  do ik=1,Kmesh%nibz ! Add tol8 for portability reasons.
697    write(msg,fmt) ik,') ',(Kmesh%ibz(ii,ik),ii=1,3),Kmesh%wt(ik)+tol8
698    call wrtout(my_unt,msg,my_mode)
699  end do
700 
701  SELECT CASE (Kmesh%timrev)
702  CASE (1)
703    write(msg,'(2a,i2,3a,i5,a)')ch10,&
704 &    ' Together with ',Kmesh%nsym,' symmetry operations (time-reversal symmetry not used) ',ch10,&
705 &    ' yields ',Kmesh%nbz,' points in the full Brillouin Zone.'
706 
707  CASE (2)
708    write(msg,'(2a,i2,3a,i5,a)')ch10,&
709 &    ' Together with ',Kmesh%nsym,' symmetry operations and time-reversal symmetry ',ch10,&
710 &    ' yields ',Kmesh%nbz,' points in the full Brillouin Zone.'
711 
712  CASE DEFAULT
713    MSG_BUG(sjoin('Wrong value for timrev:', itoa(Kmesh%timrev)))
714  END SELECT
715 
716  call wrtout(my_unt,msg,my_mode)
717 
718  if (my_prtvol>0) then
719    write(fmt,*)'(1x,i5,a,2x,3es16.8)'
720    do ik=1,Kmesh%nbz
721      if (my_prtvol==1 .and. ik>nmaxk) then
722        write(msg,'(a)')' prtvol=1, do not print more points.'
723        call wrtout(my_unt,msg,my_mode) ; EXIT
724      end if
725      write(msg,fmt)ik,') ',(Kmesh%bz(ii,ik),ii=1,3)
726      call wrtout(my_unt,msg,my_mode)
727    end do
728  end if
729  !
730  ! === Additional printing ===
731  if (my_prtvol>=10) then
732    write(msg,'(2a)')ch10,&
733 &   '                  Full point  ------->    Irred point -->            through:  Symrec  Time-Rev (1=No,-1=Yes) G0(1:3) '
734    call wrtout(my_unt,msg,my_mode)
735    write(fmt,*)'(2x,i5,2x,2(3(f7.4,2x)),i3,2x,i2,3(i3))'
736    do ik=1,Kmesh%nbz
737      write(msg,fmt)ik,Kmesh%bz(:,ik),Kmesh%ibz(:,Kmesh%tab(ik)),Kmesh%tabo(ik),Kmesh%tabi(ik),Kmesh%umklp(:,ik)
738      call wrtout(my_unt,msg,my_mode)
739    end do
740  end if
741 
742  write(msg,'(a)')ch10
743  call wrtout(my_unt,msg,my_mode)
744 
745 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

 84  type,public :: kmesh_t
 85 
 86   !scalars
 87   integer :: nshift=0
 88 
 89   integer :: nbz=0
 90   ! Number of points in the BZ.
 91 
 92   integer :: nibz=0
 93   ! Number of points in the IBZ.
 94 
 95   integer :: nsym
 96   ! Number of symmetry operations.
 97 
 98   integer :: kptopt
 99   ! kptopt=option for the generation of k points (see input variable description)
100   !
101   ! 1  if both time-reversal and point group symmetries are used.
102   ! 2  if only time-reversal symmetry is used.
103   ! 3  do not take into account any symmetry (except the identity).
104   ! 4  if time-reversal is not used (spin-orbit coupling).
105   ! <0 number of segments used to construct the k-path for NSCF calculation.
106 
107   integer :: timrev
108   ! 2 if time reversal symmetry can be used, 1 otherwise.
109 
110  !arrays
111   integer :: kptrlatt(3,3) = NONE_KPTRLATT
112    ! Coordinates of three vectors in real space, expressed in reduced coordinates.
113    ! They define a super-lattice in real space. The k point lattice is the reciprocal of
114    ! this super-lattice, eventually shifted by shift.
115    ! Not available if the structure is initialized from the points in the IBZ.
116 
117   integer,allocatable :: rottb(:,:,:)
118   ! rottb(nbz,timrev,nsym),
119   ! Index of (IS)k in the BZ array where S is a sym operation in reciprocal space,
120   ! I is the identity or the inversion operator (1,2 resp)
121 
122   integer,allocatable :: rottbm1(:,:,:)
123   ! rottbm1(nbz,timrev,nsym)
124   ! Index of IS^{-1} k in the BZ array.
125 
126   integer,allocatable :: tab(:)
127   ! tab(nbz)
128   ! For each point in the BZ, it gives the index of the symmetric irreducible point in the array ibz.
129 
130   integer,allocatable :: tabi(:)
131   ! tabi(nbz)
132   ! For each point in the BZ, tabi tells whether time-reversal has to be
133   ! used to obtain k_BZ starting from the corresponding point in the IBZ  (1=>no, -1=>yes)
134 
135   integer,allocatable :: tabo(:)
136   ! tabo(nbz)
137   ! For each point in the BZ, it gives the index in the array symrec of the
138   ! symmetry operation in reciprocal space which rotates k_IBZ onto \pm k_BZ (depending on tabi)
139 
140   integer,allocatable :: umklp(:,:)
141    ! umklp(3,nbz)
142    ! The Umklapp G0-vector such as kbz + G0 = (IS) k_ibz, where kbz is in the first BZ.
143 
144   real(dp) :: gmet(3,3)
145   ! Reciprocal space metric ($\textrm{bohr}^{-2}$).
146 
147   real(dp) :: gprimd(3,3)
148   ! Dimensional primitive translations for reciprocal space ($\textrm{bohr}^{-1}$)
149 
150   real(dp),allocatable :: bz(:,:)
151   ! bz(3,nbz)
152   ! Points in the BZ in reduced coordinates.
153   ! TODO should packed in shells.
154 
155   real(dp),allocatable :: ibz(:,:)
156   ! ibz(3,nibz)
157   ! Points in the IBZ in reduced coordinates.
158 
159   real(dp),allocatable :: shift(:,:)
160   !  shift(3,nshift)
161   !  Shift for k-points, not available is nshift=0. Usually nshift=1
162 
163   real(dp),allocatable :: wt(:)
164   ! wt(nibz)
165   ! Weights for each point in the IBZ.
166 
167   !%real(dp),allocatable :: vbox(:)
168   ! vbox(nkbz)
169   ! Volume of the small box centered on the k-point in the full BZ.
170   ! Mainly used for inhomogeneous meshes.
171 
172   complex(dpc),allocatable :: tabp(:)
173   ! tabp(nkbz)
174   ! For each point in the BZ, this table gives the phase factors associated
175   ! to non-symmorphic operations, i.e., e^{-i2\pi k_IBZ.R{^-1}t}=e^{-i2\pi k_BZ cdot t}
176   ! where \transpose R{-1}=S and  (S k_IBZ)=\pm k_BZ (depending on tabi)
177 
178  end type kmesh_t
179 
180  public :: kmesh_init            ! Main creation method.
181  public :: kmesh_free            ! Free memory
182  public :: kmesh_print           ! Printout of basic info on the object.
183  public :: get_bz_item           ! Get point in the  BZ as well as useful quantities.
184  public :: get_IBZ_item          ! Get point in the IBZ as well as useful quantities.
185  public :: get_BZ_diff           ! Get the difference k1-k2 in the BZ (if any).
186  public :: isamek                ! Check whether two points are equal within an umklapp G0.
187  public :: isequalk              ! Check whether two points are equal within an umklapp G0 (does not report G0)
188  public :: has_BZ_item           ! Check if a point belongs to the BZ mesh.
189  public :: has_IBZ_item          ! Check if a point is in the IBZ
190  public :: bz_mesh_isirred       ! TRUE. if ikbz is in the IBZ (a non-zero umklapp is not allowed)
191  public :: make_mesh             ! Initialize the mesh starting from kptrlatt and shift.
192  public :: identk                ! Find the BZ starting from the irreducible k-points.
193  public :: get_ng0sh             ! Calculate the smallest box in RSpace able to treat all possible umklapp processes.
194  public :: find_qmesh            ! Find the Q-mesh defined as the set of all possible k1-k2 differences.
195  public :: findnq                ! Helper routine returning the number of q-points.
196  public :: findq                 ! Helper routine returning the list of q-points.
197  public :: findqg0               ! Identify q + G0 = k1-k2.
198  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

PARENTS

      m_ebands,m_gruneisen,m_ifc,m_phgamma,m_phonons,wfk_analyze

CHILDREN

SOURCE

3331 subroutine kpath_free(Kpath)
3332 
3333 
3334 !This section has been created automatically by the script Abilint (TD).
3335 !Do not modify the following lines by hand.
3336 #undef ABI_FUNC
3337 #define ABI_FUNC 'kpath_free'
3338 !End of the abilint section
3339 
3340  implicit none
3341 
3342 !Arguments ------------------------------------
3343 !scalars
3344  type(kpath_t),intent(inout) :: Kpath
3345 
3346 ! *************************************************************************
3347 
3348  if (allocated(Kpath%ndivs)) then
3349    ABI_FREE(Kpath%ndivs)
3350  end if
3351 
3352  if (allocated(Kpath%bounds2kpt)) then
3353    ABI_FREE(Kpath%bounds2kpt)
3354  end if
3355 
3356  if (allocated(Kpath%bounds)) then
3357    ABI_FREE(Kpath%bounds)
3358  end if
3359 
3360  if (allocated(Kpath%points)) then
3361    ABI_FREE(Kpath%points)
3362  end if
3363 
3364  if (allocated(Kpath%dl)) then
3365    ABI_FREE(Kpath%dl)
3366  end if
3367 
3368 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.

OUTPUT

  Kpath<type(kpath_t)>=Object with the normalized path.

PARENTS

      wfk_analyze

CHILDREN

SOURCE

3255 type(kpath_t) function kpath_new(bounds, gprimd, ndivsm) result(kpath)
3256 
3257 
3258 !This section has been created automatically by the script Abilint (TD).
3259 !Do not modify the following lines by hand.
3260 #undef ABI_FUNC
3261 #define ABI_FUNC 'kpath_new'
3262 !End of the abilint section
3263 
3264  implicit none
3265 
3266 !Arguments ------------------------------------
3267 !scalars
3268  integer,intent(in) :: ndivsm
3269 !!arrays
3270  real(dp),intent(in) :: bounds(:,:),gprimd(3,3)
3271 
3272 !Local variables-------------------------------
3273  integer :: ii
3274 !arrays
3275  real(dp) :: dk(3)
3276  real(dp),allocatable :: pts(:,:)
3277 
3278 ! *************************************************************************
3279 
3280  ABI_CHECK(size(bounds, dim=1) == 3, "Wrong dim1 in bounds")
3281  ABI_CHECK(ndivsm > 0, sjoin("ndivsm:", itoa(ndivsm)))
3282  Kpath%nbounds = size(bounds, dim=2)
3283  Kpath%ndivsm = ndivsm
3284 
3285  ! Compute reciprocal space metric.
3286  Kpath%gprimd = gprimd; Kpath%gmet = matmul(transpose(gprimd), gprimd)
3287 
3288  ABI_MALLOC(Kpath%ndivs, (Kpath%nbounds-1))
3289  call make_path(Kpath%nbounds,bounds,Kpath%gmet,"G",ndivsm,Kpath%ndivs,Kpath%npts,pts,unit=dev_null)
3290 
3291  ABI_MALLOC(Kpath%bounds, (3,Kpath%nbounds))
3292  Kpath%bounds = bounds
3293 
3294  ABI_MALLOC(Kpath%points, (3,Kpath%npts))
3295  Kpath%points = pts
3296  ABI_FREE(pts)
3297 
3298  ! Compute distance between point i-1 and i
3299  ABI_CALLOC(kpath%dl, (kpath%npts))
3300  do ii=2,kpath%npts
3301    dk = kpath%points(:, ii-1) - kpath%points(:,ii)
3302    kpath%dl(ii) = normv(dk, kpath%gmet, "G")
3303  end do
3304 
3305  ! Mapping bounds --> points
3306  ABI_MALLOC(kpath%bounds2kpt, (kpath%nbounds))
3307  kpath%bounds2kpt(1) = 1
3308  do ii=1,kpath%nbounds-1
3309    kpath%bounds2kpt(ii+1) = sum(kpath%ndivs(:ii)) + 1
3310  end do
3311 
3312 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

PARENTS

      m_ebands,m_gruneisen,m_phgamma

CHILDREN

SOURCE

3396 subroutine kpath_print(kpath, header, unit, prtvol, pre)
3397 
3398 
3399 !This section has been created automatically by the script Abilint (TD).
3400 !Do not modify the following lines by hand.
3401 #undef ABI_FUNC
3402 #define ABI_FUNC 'kpath_print'
3403 !End of the abilint section
3404 
3405  implicit none
3406 
3407 !Arguments ------------------------------------
3408 !scalars
3409  integer,optional,intent(in) :: unit,prtvol
3410  character(len=*),optional,intent(in) :: header,pre
3411  type(kpath_t),intent(in) :: kpath
3412 
3413 !Local variables-------------------------------
3414  integer :: unt,my_prtvol,ii
3415  character(len=500) :: my_pre
3416 
3417 ! *************************************************************************
3418 
3419  unt = std_out; if (present(unit)) unt = unit
3420  my_prtvol = 0; if (present(prtvol)) my_prtvol = prtvol
3421  my_pre = " "; if (present(pre)) my_pre = pre
3422  if (unt <= 0) return
3423 
3424  if (present(header)) write(unt,"(a)") sjoin(my_pre, '==== '//trim(adjustl(header))//' ==== ')
3425  write(unt, "(a)") sjoin(my_pre, "Number of points:", itoa(kpath%npts), ", ndivsmall:", itoa(kpath%ndivsm))
3426  write(unt, "(a)") sjoin(my_pre, "Boundaries and corresponding index in the k-points array:")
3427  do ii=1,kpath%nbounds
3428    write(unt, "(a)") sjoin(my_pre, itoa(kpath%bounds2kpt(ii)), ktoa(kpath%bounds(:,ii)))
3429  end do
3430  write(unt, "(a)") sjoin(my_pre, " ")
3431 
3432  if (my_prtvol > 10) then
3433    do ii=1,kpath%npts
3434      write(unt, "(a)") sjoin(my_pre, ktoa(kpath%points(:,ii)))
3435    end do
3436  end if
3437 
3438 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

212 type,public :: kpath_t
213 
214   integer :: nbounds=0
215     ! Number of extrema defining the path.
216 
217   integer :: ndivsm=0
218     ! Number of divisions used to sample the smallest segment.
219 
220   integer :: npts=0
221     ! Total number of points in the path.
222 
223   real(dp) :: gprimd(3,3)
224    ! Reciprocal lattice vectors.
225 
226   real(dp) :: gmet(3,3)
227    ! Metric matrix in G space.
228 
229   integer,allocatable :: ndivs(:)
230    ! ndivs(nbounds-1)
231    ! Number of division for each segment.
232 
233   integer,allocatable :: bounds2kpt(:)
234    ! bounds2kpt(nbounds)
235    ! bounds2kpt(i): Index of the i-th extrema in the pts(:) array.
236 
237   real(dp),allocatable :: bounds(:,:)
238     ! bounds(3,nbounds)
239     ! The points defining the path in reduced coordinates.
240 
241   real(dp),allocatable :: points(:,:)
242     ! points(3,npts)
243     ! The points of the path in reduced coordinates.
244 
245   real(dp),allocatable :: dl(:)
246     ! dl(npts)
247     ! dl(i) = Distance between the (i-1)-th and the i-th k-point. dl(1) = zero
248 
249  end type kpath_t
250 
251  public :: kpath_new        ! Construct a new path
252  public :: kpath_free       ! Free memory
253  public :: make_path        ! Construct a normalized path. TODO: Remove
254  public :: kpath_print      ! Print the path.

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

PARENTS

      m_bz_mesh

CHILDREN

SOURCE

2947 subroutine littlegroup_free_0D(Ltg)
2948 
2949 
2950 !This section has been created automatically by the script Abilint (TD).
2951 !Do not modify the following lines by hand.
2952 #undef ABI_FUNC
2953 #define ABI_FUNC 'littlegroup_free_0D'
2954 !End of the abilint section
2955 
2956  implicit none
2957 
2958 !Arguments ------------------------------------
2959 !scalars
2960  type(littlegroup_t),intent(inout) :: Ltg
2961 
2962 ! *********************************************************************
2963 
2964  !@littlegroup_t
2965  if (allocated(Ltg%g0)) then
2966    ABI_FREE(Ltg%g0)
2967  end if
2968  if (allocated(Ltg%ibzq)) then
2969    ABI_FREE(Ltg%ibzq)
2970  end if
2971  if (allocated(Ltg%bz2ibz)) then
2972    ABI_FREE(Ltg%bz2ibz)
2973  end if
2974  if (allocated(Ltg%ibz2bz)) then
2975    ABI_FREE(Ltg%ibz2bz)
2976  end if
2977  if (allocated(Ltg%igmG0)) then
2978    ABI_FREE(Ltg%igmG0)
2979  end if
2980  if (allocated(Ltg%flag_umklp)) then
2981    ABI_FREE(Ltg%flag_umklp)
2982  end if
2983  if (allocated(Ltg%preserve)) then
2984    ABI_FREE(Ltg%preserve)
2985  end if
2986  if (allocated(Ltg%tab)) then
2987    ABI_FREE(Ltg%tab)
2988  end if
2989  if (allocated(Ltg%tabo)) then
2990    ABI_FREE(Ltg%tabo)
2991  end if
2992  if (allocated(Ltg%tabi)) then
2993    ABI_FREE(Ltg%tabi)
2994  end if
2995  if (allocated(Ltg%wtksym)) then
2996    ABI_FREE(Ltg%wtksym)
2997  end if
2998 
2999 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

PARENTS

CHILDREN

SOURCE

3022 subroutine littlegroup_free_1D(Ltg)
3023 
3024 
3025 !This section has been created automatically by the script Abilint (TD).
3026 !Do not modify the following lines by hand.
3027 #undef ABI_FUNC
3028 #define ABI_FUNC 'littlegroup_free_1D'
3029 !End of the abilint section
3030 
3031  implicit none
3032 
3033 !Arguments ------------------------------------
3034 !scalars
3035  type(littlegroup_t),intent(inout) :: Ltg(:)
3036 
3037 !Local variables-------------------------------
3038  integer :: ipt
3039 
3040 ! *********************************************************************
3041 
3042  do ipt=1,SIZE(Ltg)
3043    call littlegroup_free_0D(Ltg(ipt))
3044  end do
3045 
3046 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
 Kmesh<kmesh_t>
   %nbz=number of points in the full BZ
   %kbz(3,nbz)=points in the full Brillouin Zon
 Cryst<crystal_t>= Info on symmetries and unit cell
   %symrec(3,3,nsym)=symmetry operations in reciprocal space
   %nsym=number of symmetry operations in the space group
   %timrev=if 2 time-reversal can be used; 1 otherwise
   %gmet(3,3)=reciprocal space metric (bohr**-2).
 npwvec=number of G vectors
 gvec(3,npwvec) coordinates of G vectors
 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.
 use_umklp=flag to include umklapp G0 vectors in the definition of the little group (0:n0,1:yes)

OUTPUT

 Ltg% <littlegroup_t_datatype>.
  %ibzq(nbz)= 1 if the kpoint belongs to the IBZ defined by ext_pt, 0 otherwise
  %bz2ibz(nbz)= sequential index of the point in the IBZ defined by ext_pt
  %ibz2bz(nibz_Ltg) For each nibz_Ltg the correspondind index in the BZ array
  %igmG0(npwepstimrev,nsym)= index of the uklapp vector G_o in the FFT array
  %flag_umklp(timrev,nsym)= flag for umklapp processes
    1 if operation (IS) requires a G_o to preserve ext_pt, 0 otherwise
  %tab(nbz)=table giving, for each k-point in the BZ (kBZ), the corresponding
   irreducible point (kIBZ) in the irreducible zone defined by the little group of ext_pt,
   i.e kBZ= (IS) kIBZ where I is either the inversion or the identity and S is an
   operation in the little group defined by ext_pt
  %tabo(nbz)=the symmetry operation S in the little group that takes kIBZ to each kBZ
  %tabi(nbz)= defines whether inversion has to be considered in the
   relation kBZ=(IS) kIBZ (1 => only S; -1 => -S)
  %preserve(2,nsym)= 1 if ISq=q, 0 otherwise, the first index is for the identity or the time reversal symmetry,
  %wtksym(2,nsym,nbz)= for each kpoint is equal to 1 if the symmetry operation (with or without time reversal)
   must be considered in the calculation of \chi_o, 0 otherwise

PARENTS

      cchi0q0_intraband,setup_screening,sigma

CHILDREN

SOURCE

2614 subroutine littlegroup_init(ext_pt,Kmesh,Cryst,use_umklp,Ltg,npwe,gvec)
2615 
2616 
2617 !This section has been created automatically by the script Abilint (TD).
2618 !Do not modify the following lines by hand.
2619 #undef ABI_FUNC
2620 #define ABI_FUNC 'littlegroup_init'
2621 !End of the abilint section
2622 
2623  implicit none
2624 
2625 !Arguments ------------------------------------
2626 !scalars
2627  integer,intent(in) :: npwe,use_umklp
2628  type(crystal_t),target,intent(in) :: Cryst
2629  type(kmesh_t),intent(in) :: Kmesh
2630  type(littlegroup_t),intent(inout) :: Ltg
2631 !arrays
2632  !integer,optional,intent(in) :: gvec(3,npwvec)
2633  integer,optional,intent(in) :: gvec(:,:)
2634  real(dp),intent(in) :: ext_pt(3)
2635 
2636 !Local variables-------------------------------
2637 !scalars
2638  integer :: dummy_timrev,enough,idx,ige,igpw,ik,ind,iold,iout,isym,itest,itim
2639  integer :: nbz,nkibzq,nsym,nsym_Ltg,ntest,timrev,ierr,npwvec
2640  real(dp) :: G0len,kin,mG0len,max_kin
2641  logical :: found,found_identity,ltest,use_antiferro
2642  character(len=500) :: msg
2643 !arrays
2644  integer :: g0(3),gg(3),gmG0(3),identity(3,3),nop(Cryst%timrev),nopg0(2)
2645  integer :: symxpt(4,2,Cryst%nsym)
2646  integer,allocatable :: indkpt1(:),symafm_ltg(:),symrec_Ltg(:,:,:)
2647  integer,pointer :: symafm(:),symrec(:,:,:)
2648  real(dp) :: knew(3)
2649  real(dp),allocatable :: ktest(:,:),wtk(:),wtk_folded(:)
2650 
2651 !************************************************************************
2652 
2653  DBG_ENTER("COLL")
2654  !
2655  ! !@littlegroup_t
2656  ! === Initial check ====
2657  ltest=(Cryst%timrev==1.or.Cryst%timrev==2)
2658  ABI_CHECK(ltest,'Wrong value for timrev')
2659  !
2660  ! === Get useful data ===
2661  nsym          =  Cryst%nsym
2662  !timrev       = 1
2663  timrev        =  Cryst%timrev
2664  symrec        => Cryst%symrec
2665  symafm        => Cryst%symafm
2666  use_antiferro =  Cryst%use_antiferro
2667 
2668  nbz =  Kmesh%nbz
2669  !
2670  ! === Destroy structure if it already exists ===
2671  call littlegroup_free(Ltg)
2672  !
2673  ! === Store dimensions and useful info ===
2674  Ltg%nsym_sg  =nsym
2675  Ltg%timrev   =timrev
2676  Ltg%nbz      =nbz
2677  !Ltg%use_umklp=use_umklp ! 1 if umklapp processes are used
2678  Ltg%ext_pt(:)=ext_pt(:)
2679 
2680  ABI_MALLOC(Ltg%G0,(3,timrev,nsym))
2681  ABI_MALLOC(Ltg%ibzq,(nbz))
2682  ABI_MALLOC(Ltg%bz2ibz,(nbz))
2683  ABI_MALLOC(Ltg%preserve,(timrev,nsym))
2684  ABI_MALLOC(Ltg%wtksym,(timrev,nsym,nbz))
2685  ABI_MALLOC(Ltg%tab,(nbz))
2686  ABI_MALLOC(Ltg%tabi,(nbz))
2687  ABI_MALLOC(Ltg%tabo,(nbz))
2688  ABI_MALLOC(Ltg%flag_umklp,(timrev,nsym))
2689  !
2690  ! In the old GW implementation we were removing symmetries related by time-reversal and
2691  ! sometimes it happened that only the inversion was reported in the KSS file (see outkss.F90).
2692  identity(:,:)=RESHAPE((/1,0,0,0,1,0,0,0,1/),(/3,3/)) ; found_identity=.FALSE.
2693  do isym=1,nsym
2694    if (ALL(symrec(:,:,isym)==identity)) then
2695     found_identity=.TRUE. ; EXIT
2696    end if
2697  end do
2698  if (.not.found_identity) then
2699    write(msg,'(5a)')&
2700 &    'Only the inversion was found in the set of symmetries read from the KSS file ',ch10,&
2701 &    'Likely you are using a KSS file generated with an old version of Abinit, ',ch10,&
2702 &    'To run a GW calculation with an old KSS file, use version < 5.5 '
2703    MSG_ERROR(msg)
2704  end if
2705 
2706  ! Find operations in the little group as well as umklapp vectors G0 ===
2707  call littlegroup_q(nsym,ext_pt,symxpt,symrec,symafm,dummy_timrev,prtvol=0)
2708 
2709  Ltg%preserve(:,:)=0; Ltg%g0(:,:,:)=0; Ltg%flag_umklp(:,:)=0; mG0len=zero
2710 
2711  do itim=1,timrev
2712    do isym=1,nsym
2713      if (symafm(isym)==-1) CYCLE
2714 
2715      if (symxpt(4,itim,isym)==1) then  !\pm Sq = q+g0
2716        if (ANY(symxpt(1:3,itim,isym)/=0).and.use_umklp==0) CYCLE ! Exclude non zero G0 vectors
2717        Ltg%preserve(itim,isym)=1
2718        g0(:)=symxpt(1:3,itim,isym); Ltg%g0(:,itim,isym)=g0(:)
2719        if (ANY(Ltg%g0(:,itim,isym)/=0)) Ltg%flag_umklp(itim,isym)=1
2720        ! Max radius to be considered to include all G0s
2721        G0len=normv(g0,Cryst%gmet,'G')
2722        mG0len=MAX(mG0len,G0len)
2723      end if
2724    end do
2725  end do
2726 
2727  nop(:)=0; nopg0(:)=0
2728  do itim=1,timrev
2729    nop  (itim)=SUM(Ltg%preserve  (itim,:))
2730    nopg0(itim)=SUM(Ltg%flag_umklp(itim,:))
2731  end do
2732  nsym_Ltg=SUM(nop(:))
2733 
2734  ! Store little group operations, include time-reversal if present ===
2735  Ltg%nsym_Ltg=nsym_Ltg
2736  ABI_MALLOC(symrec_Ltg,(3,3,Ltg%nsym_Ltg))
2737 
2738  ind=1
2739  do itim=1,timrev
2740    do isym=1,nsym
2741      if (Ltg%preserve(itim,isym)==1) then
2742       if (itim==1) symrec_Ltg(:,:,ind)= symrec(:,:,isym)
2743       if (itim==2) symrec_Ltg(:,:,ind)=-symrec(:,:,isym)
2744       ind=ind+1
2745      end if
2746    end do
2747  end do
2748 
2749  ! Check the closure of the (ferromagnetic) little group ===
2750  ABI_MALLOC(symafm_ltg,(Ltg%nsym_Ltg))
2751  symafm_ltg(:)=1
2752  call chkgrp(Ltg%nsym_Ltg,symafm_ltg,symrec_Ltg,ierr)
2753  ABI_CHECK(ierr==0,"Error in group closure")
2754 
2755  ABI_FREE(symafm_ltg)
2756  !
2757  ! Find the irreducible zone associated to ext_pt
2758  ! Do not use time-reversal since it has been manually introduced previously
2759  ABI_MALLOC(indkpt1,(nbz))
2760  ABI_MALLOC(wtk_folded,(nbz))
2761  ABI_MALLOC(wtk,(nbz))
2762  wtk=one; iout=0; dummy_timrev=0
2763 
2764  call symkpt(0,Cryst%gmet,indkpt1,iout,Kmesh%bz,nbz,nkibzq,Ltg%nsym_Ltg,symrec_Ltg,dummy_timrev,wtk,wtk_folded)
2765 
2766  ABI_FREE(indkpt1)
2767  ABI_FREE(wtk)
2768 
2769  Ltg%nibz_Ltg=nkibzq
2770  !
2771  ! === Set up table in the BZ ===
2772  ! * 0 if the point does not belong to IBZ_xpt, 1 otherwise
2773  ABI_MALLOC(Ltg%ibz2bz,(nkibzq))
2774  Ltg%ibzq(:)=0; Ltg%bz2ibz(:)=0; Ltg%ibz2bz(:)=0
2775 
2776  ind=0; enough=0
2777  do ik=1,nbz
2778    if (wtk_folded(ik)>tol8) then
2779      ind=ind+1
2780      Ltg%ibzq(ik)=1
2781      Ltg%bz2ibz(ik) =ind
2782      Ltg%ibz2bz(ind)=ik
2783    end if
2784  end do
2785  ABI_CHECK(ind==Ltg%nibz_Ltg," BUG ind/=Ltg%nibz_Ltg")
2786  !
2787  ! Reconstruct full BZ starting from IBZ_q.
2788  ! Calculate appropriate weight for each item (point,symmetry operation,time-reversal)
2789  Ltg%tab=0; Ltg%tabo=0; Ltg%tabi=0; Ltg%wtksym(:,:,:)=0
2790 
2791  ! Start with zero no. of k-points found
2792  ntest=0
2793  ABI_MALLOC(ktest,(3,nbz))
2794  ktest=zero
2795 
2796  do ik=1,nbz
2797    if (Ltg%ibzq(ik)/=1) CYCLE
2798    ! * Loop over symmetry operations S and time-reversal.
2799    ! * Use spatial inversion instead of time reversal whenever possible.
2800    do itim=1,timrev
2801      do isym=1,nsym
2802 
2803       ! Form IS k only for (IS) pairs in the (ferromagnetic) little group.
2804       if (symafm(isym)==-1) CYCLE
2805       if (Ltg%preserve(itim,isym)==0) CYCLE
2806       knew(:)=(3-2*itim)*MATMUL(symrec(:,:,isym),Kmesh%bz(:,ik))
2807       !
2808       ! === Check whether it has already been found (to within a RL vector) ===
2809       iold=0
2810       do itest=1,ntest
2811         if (isamek(knew(:),ktest(:,itest),gg)) iold=iold+1
2812       end do
2813 
2814       if (iold==0) then
2815         ! == Found new BZ point ===
2816         ! For this point the operation (isym,itim) must be considered to reconstruct the full BZ
2817         Ltg%wtksym(itim,isym,ik)=1
2818         ntest=ntest+1
2819         ktest(:,ntest)=knew(:)
2820         !
2821         ! === Now find knew in the BZ array ===
2822         found=.FALSE.
2823         do idx=1,nbz
2824           if (isamek(knew(:),Kmesh%bz(:,idx),gg)) then ! They are the same within a RL vector
2825             Ltg%tab (idx)=ik
2826             Ltg%tabo(idx)=isym
2827             Ltg%tabi(idx)=3-2*itim
2828             found=.TRUE.; EXIT
2829           end if
2830         end do
2831         if (.not.found) then
2832           write(msg,'(a,3f12.6,a)')'Not able to find the ',knew(:),' in the array BZ '
2833           MSG_ERROR(msg)
2834         end if
2835       end if
2836 
2837      end do !isym
2838    end do !itim
2839 
2840  end do !nbz
2841 
2842  ABI_FREE(ktest)
2843 
2844  if (ntest/=nbz) then
2845    MSG_BUG(sjoin('ntest - nbz = ',itoa(ntest-nbz)))
2846  end if
2847 
2848  if (SUM(Ltg%wtksym)/=nbz) then
2849    MSG_BUG(sjoin('sum(Ltg%wtksym)-nbz = ', itoa(SUM(Ltg%wtksym)-nbz)))
2850  end if
2851 
2852  Ltg%max_kin_gmG0=zero
2853 
2854  if (npwe>0.and.PRESENT(gvec)) then
2855    npwvec = SIZE(gvec,DIM=2)
2856    ! This correspond to the case in which we need to know the index of G-Go in the gvec array
2857    ! where G is one of the npwe vectors. This is required in screening but not in sigma.
2858    ! The drawback is that the effective G sphere used to calculate the oscillators must be smaller
2859    ! that gvec if we want to avoid possible aliasing effects. Lifting this constraint would require
2860    ! a lot of boring coding. (no need to do this if ext_pt=zero, but oh well)
2861    ABI_MALLOC(Ltg%igmG0,(npwe,timrev,nsym))
2862    Ltg%igmG0(:,:,:)=0
2863    max_kin=zero
2864 
2865    ! Loop over symmetry operations S and time-reversal
2866    do itim=1,timrev
2867      do isym=1,nsym
2868        ! Form IS k only for (IS) pairs in the little group
2869        if (symafm(isym)==-1) CYCLE
2870        if (Ltg%preserve(itim,isym)/=0) then
2871          g0(:)=Ltg%g0(:,itim,isym)
2872          do ige=1,npwe
2873            gmG0(:)=gvec(:,ige)-g0(:)
2874            kin=half*normv(gmG0,Cryst%gmet,'G')**2
2875            max_kin=MAX(max_kin,kin)
2876 
2877            found=.FALSE.
2878            do igpw=1,npwvec
2879              if (ALL(gvec(:,igpw)-gmG0(:)==0)) then
2880                Ltg%igmG0(ige,itim,isym)=igpw
2881                found=.TRUE.; EXIT
2882              end if
2883            end do
2884            if (.not.found) then
2885              write(msg,'(5a,f8.3,2a,3i5)')&
2886 &              'Not able to found G-G0 in the largest G-spere ',ch10,&
2887 &              'Decrease the size of epsilon or, if possible, increase ecutwfn (>ecuteps) ',ch10,&
2888 &              'Minimum required cutoff energy for G-G0 sphere= ',kin,ch10,&
2889 &              'G0 = ',g0(:)
2890              MSG_ERROR(msg)
2891            end if
2892          end do
2893        end if
2894      end do
2895    end do
2896    Ltg%max_kin_gmG0=max_kin
2897  end if
2898  ABI_FREE(symrec_Ltg)
2899 
2900 #ifdef DEBUG_MODE
2901  do ik=1,nbz
2902    if (ABS(SUM(Ltg%wtksym(1,:,ik)+Ltg%wtksym(2,:,ik))-wtk_folded(ik))>tol6) then
2903      write(std_out,*)' sum(Ltg%wtksym,ik)-wtk_folded(ik) = ',sum(Ltg%wtksym(1,:,ik)+Ltg%wtksym(2,:,ik))-wtk_folded(ik)
2904      write(std_out,*)Ltg%wtksym(1,:,ik),Ltg%wtksym(2,:,ik),wtk_folded(ik)
2905      write(std_out,*)ik,Kmesh%bz(:,ik)
2906      MSG_BUG("Wrong weight")
2907    end if
2908  end do
2909  do ik=1,nbz
2910    knew = Ltg%tabi(ik) * MATMUL(symrec(:,:,Ltg%tabo(ik)),Kmesh%bz(:,Ltg%tab(ik)))
2911    if (.not.isamek(knew,Kmesh%bz(:,ik),gg)) then
2912      write(std_out,*)knew,Kmesh%bz(:,ik)
2913      write(std_out,*)Ltg%tabo(ik),Ltg%tabi(ik),Ltg%tab(ik)
2914      MSG_BUG("Wrong tables")
2915    end if
2916  end do
2917 #endif
2918 
2919  ABI_FREE(wtk_folded)
2920 
2921  DBG_EXIT("COLL")
2922 
2923 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

PARENTS

      calc_sigc_me,calc_sigx_me,cchi0,cchi0q0,cchi0q0_intraband,cohsex_me

CHILDREN

SOURCE

3074 subroutine littlegroup_print(Ltg,unit,prtvol,mode_paral)
3075 
3076 
3077 !This section has been created automatically by the script Abilint (TD).
3078 !Do not modify the following lines by hand.
3079 #undef ABI_FUNC
3080 #define ABI_FUNC 'littlegroup_print'
3081 !End of the abilint section
3082 
3083  implicit none
3084 
3085 !Arguments ------------------------------------
3086 !scalars
3087  integer,optional,intent(in) :: prtvol,unit
3088  character(len=4),optional,intent(in) :: mode_paral
3089  type(littlegroup_t),intent(in) :: Ltg
3090 
3091 !Local variables-------------------------------
3092 !scalars
3093  integer :: itim,my_unt,my_prtvol
3094  character(len=4) :: my_mode
3095  character(len=500) :: msg
3096 !arrays
3097  integer :: nop(Ltg%timrev),nopg0(Ltg%timrev)
3098 
3099 ! *********************************************************************
3100 
3101  my_unt   =std_out; if (PRESENT(unit      )) my_unt   =unit
3102  my_prtvol=0      ; if (PRESENT(prtvol    )) my_prtvol=prtvol
3103  my_mode  ='COLL' ; if (PRESENT(mode_paral)) my_mode  =mode_paral
3104 
3105  write(msg,'(4a,3es16.8,2a,i5,a,i5,2a,i3,a,i3)')ch10,&
3106 &  ' ==== Little Group Info ==== ',ch10,&
3107 &  '  External point ',Ltg%ext_pt,ch10,&
3108 &  '  Number of points in the IBZ defined by little group  ',Ltg%nibz_Ltg,'/',Ltg%nbz,ch10,&
3109 &  '  Number of operations in the little group : ',Ltg%nsym_Ltg,'/',Ltg%nsym_sg
3110  call wrtout(my_unt,msg,my_mode)
3111 
3112  nop=0 ; nopg0=0
3113  do itim=1,Ltg%timrev
3114    nop  (itim)=SUM(Ltg%preserve  (itim,:))
3115    nopg0(itim)=SUM(Ltg%flag_umklp(itim,:))
3116  end do
3117 
3118  do itim=1,Ltg%timrev
3119    if (itim==1) then
3120      write(msg,'(a,2(a,i2,a))')ch10,&
3121 &      '  No time-reversal symmetry with zero umklapp: ',nop(1)-nopg0(1),ch10,&
3122 &      '  No time-reversal symmetry with non-zero umklapp: ',nopg0(1),ch10
3123      call wrtout(my_unt,msg,my_mode)
3124    else if (itim==2) then
3125      write(msg,'(a,2(a,i2,a))')ch10,&
3126 &      '  time-reversal symmetry with zero umklapp: ',nop(2)-nopg0(2),ch10,&
3127 &      '  time-reversal symmetry with non-zero umklapp: ',nopg0(2),ch10
3128      call wrtout(my_unt,msg,my_mode)
3129    end if
3130  end do
3131 
3132 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

286  type,public :: littlegroup_t
287 
288   integer :: npw             ! No. of planewaves used to describe the wavefuntion, used to dimension igmG0
289   integer :: nsym_sg         ! No. of operations in the space group (*NOT* the little group)
290   integer :: nsym_ltg        ! No. of symmetry operations in the little group (time-reversal is included, if can be used)
291   integer :: timrev          ! 2 if time-reversal is considered, 1 otherwise
292   integer :: nbz             ! No. of kpoints in the full BZ
293   integer :: nibz_ltg        ! No. of points in the irreducible wedge defined by the little group
294   !integer :: use_umklp      ! 1 if umklapp processes are included
295 
296   real(dp) :: max_kin_gmG0
297   ! Max kinetic energy of G-G0 in case of umklapp.
298 
299   integer,allocatable :: G0(:,:,:)
300   ! g0(3,timrev,nsym_sg)
301   ! Reduced coordinates of the umklapp G0 vector.
302 
303   integer,allocatable :: ibzq(:)
304   ! ibzq(nbz)
305   ! 1 if the point belongs to the IBZ_q defined by ext_pt, 0 otherwise.
306 
307   integer,allocatable :: bz2ibz(:)
308   ! bz2ibz(nbz)
309   ! Index of the point in the irreducible wedge defined by the little group, 0 otherwise.
310 
311   integer,allocatable :: ibz2bz(:)
312   ! ibz2bz(nibz_ltg)
313   ! The correspondind index in the BZ array
314 
315   integer,allocatable :: igmG0(:,:,:)
316   ! iumklp(npw,timrev,nsym_sg)
317   ! Index of G-G0 in the FFT array for each operations IS (I=\pm 1).
318 
319   integer,allocatable :: flag_umklp(:,:)
320   ! flag_umklp(timrev,nsym_sg)
321   ! 1 if the operation IS requires a non null G0 vector to preserve q, 0 otherwise.
322 
323   integer,allocatable :: preserve(:,:)
324   ! preserve(timrev,nsym_sg)
325   ! preserve(1,S) is 1 if the operation S in rec space preserves the external q-point i.e Sq=q+G0
326   ! preserve(2,S) is 1 if -Sq=q+G0. G0 is a reciprocal lattice vector also called "umklapp vector".
327 
328   integer,allocatable :: tab(:)
329   ! tab(nbz)
330   ! For each point in BZ, the index of the irreducible point (kIBZ_q) in the irreducible
331   ! wedge defined by the little group of q. kBZ= (IS) kIBZ where I is the inversion or the identity.
332 
333   integer,allocatable :: tabo(:)
334   ! tabo(nbz)
335   ! The index of the operation S in the little group that rotates kIBZ_q into \pm kBZ.
336 
337   integer,allocatable :: tabi(:)
338   ! tabi(nbz)
339   ! for each k-point in the BZ defines whether inversion has to be
340   ! considered in the relation kBZ= IS kIBZ_q (1 => only S; -1 => -S).
341 
342   integer,allocatable :: wtksym(:,:,:)
343   ! wtksym(timrev,nsym_sg,kbz)
344   ! 1 if IS belongs to the little group, 0 otherwise !(should invert firt dimensions.
345 
346   real(dp) :: ext_pt(3)
347   ! The external point defining the little group.
348 
349  end type littlegroup_t
350 
351  public :: littlegroup_init
352  public :: littlegroup_free
353  public :: littlegroup_print

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.

PARENTS

      m_shirley,setup_bse,setup_bse_interp

CHILDREN

SOURCE

1460 subroutine make_mesh(Kmesh,Cryst,kptopt,kptrlatt,nshiftk,shiftk,&
1461 &  vacuum,break_symmetry)  ! Optional
1462 
1463 
1464 !This section has been created automatically by the script Abilint (TD).
1465 !Do not modify the following lines by hand.
1466 #undef ABI_FUNC
1467 #define ABI_FUNC 'make_mesh'
1468 !End of the abilint section
1469 
1470  implicit none
1471 
1472 !Arguments -------------------------------
1473 !scalars
1474  integer,intent(in) :: nshiftk,kptopt
1475  logical,optional,intent(in) :: break_symmetry
1476  type(kmesh_t),intent(inout) :: Kmesh
1477  type(crystal_t),intent(in) :: Cryst
1478 !arrays
1479  integer,intent(inout) :: kptrlatt(3,3)
1480  integer,optional,intent(in) :: vacuum(3)
1481  real(dp),intent(in) :: shiftk(3,nshiftk)
1482 
1483 !Local variables -------------------------
1484 !scalars
1485  integer,parameter :: chksymbreak0=0
1486  integer :: iscf,nkbz,nkibz,nkpt_computed,my_nshiftk
1487  real(dp) :: kptrlen
1488  logical :: my_break_symmetry
1489 !arrays
1490  integer :: my_vacuum(3)
1491  real(dp),allocatable :: kibz(:,:),wtk(:),my_shiftk(:,:),ref_kbz(:,:)
1492 
1493 ! *************************************************************************
1494 
1495  DBG_ENTER("COLL")
1496 
1497  !@kmesh_t
1498  !if (ALL(kptopt/=(/1,2,3,4/))) then
1499  if (ALL(kptopt/=(/1,3/))) then
1500    MSG_WARNING(sjoin(" Not allowed value for kptopt: ", itoa(kptopt)))
1501  end if
1502  !
1503  ! ======================================================================
1504  ! ==== First call to getkgrid to obtain nkibz as well as the BZ set ====
1505  ! ======================================================================
1506  iscf=7  ! use for the Weights in NSCF calculation. check it more carefully.
1507  nkibz=0 ! Compute number of k-points in the BZ and IBZ
1508 
1509  my_vacuum = (/0,0,0/); if (PRESENT(vacuum)) my_vacuum=vacuum
1510 
1511  my_nshiftk = nshiftk
1512  ABI_CHECK(my_nshiftk>0.and.my_nshiftk<=210,"Wrong nshiftk")
1513  ABI_MALLOC(my_shiftk,(3,210))
1514  my_shiftk=zero; my_shiftk(:,1:nshiftk) = shiftk(:,:)
1515 
1516  !write(std_out,*)" In make_mesh"
1517  !write(std_out,*)" kptopt   ",kptopt," kptrlatt ",kptrlatt
1518  !write(std_out,*)" nshiftk  ",nshiftk," shiftk   ",shiftk
1519 
1520  ABI_MALLOC(kibz,(3,nkibz))
1521  ABI_MALLOC(wtk,(nkibz))
1522 
1523  call getkgrid(chksymbreak0,0,iscf,kibz,kptopt,kptrlatt,kptrlen,Cryst%nsym,0,nkibz,my_nshiftk,&
1524 &  Cryst%nsym,Cryst%rprimd,my_shiftk,Cryst%symafm,Cryst%symrel,my_vacuum,wtk,fullbz=ref_kbz)
1525 
1526  nkbz = SIZE(ref_kbz,DIM=2)
1527 
1528  ABI_FREE(kibz)
1529  ABI_FREE(wtk)
1530 
1531  !write(std_out,*)" after getkgrid1: nkbz = ",nkbz," nkibz=",nkibz
1532  !write(std_out,*)" ref_kbz = ",ref_kbz
1533 
1534  ! ==============================================================
1535  ! ==== Recall getkgrid to get kibz(3,nkibz) and wtk(nkibz) =====
1536  ! ==============================================================
1537 
1538  ABI_MALLOC(kibz,(3,nkibz))
1539  ABI_MALLOC(wtk,(nkibz))
1540 
1541  call getkgrid(chksymbreak0,0,iscf,kibz,kptopt,kptrlatt,kptrlen,Cryst%nsym,nkibz,nkpt_computed,my_nshiftk,&
1542 &  Cryst%nsym,Cryst%rprimd,my_shiftk,Cryst%symafm,Cryst%symrel,my_vacuum,wtk)
1543 
1544  ! Store quantities that cannot be easily (and safely) calculated if we only know the IBZ.
1545  Kmesh%nshift   = my_nshiftk
1546  Kmesh%kptrlatt = kptrlatt
1547 
1548  ! Call the main creation method to get the tables tabo, tabi, tabp, umklp...
1549  ! kmesh_init will reconstruct the BZ from kibz but pruning the k-points not in ref_bz
1550  ! TODO: solve problem with timrev
1551  my_break_symmetry=.FALSE.; if (PRESENT(break_symmetry)) my_break_symmetry=break_symmetry
1552  call kmesh_init(Kmesh,Cryst,nkibz,kibz,kptopt,ref_bz=ref_kbz,break_symmetry=my_break_symmetry)
1553 
1554  ABI_FREE(ref_kbz)
1555  ABI_FREE(kibz)
1556  ABI_FREE(wtk)
1557 
1558  ABI_MALLOC(Kmesh%shift,(3,my_nshiftk))
1559  Kmesh%shift=my_shiftk(:,1:my_nshiftk)
1560  ! Init Kmesh is breaking nshiftk
1561  Kmesh%nshift=my_nshiftk
1562 
1563  ABI_FREE(my_shiftk)
1564 
1565  DBG_EXIT("COLL")
1566 
1567 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.

PARENTS

      m_bz_mesh,m_nesting,m_phonons,mkph_linwid

CHILDREN

SOURCE

2023 subroutine make_path(nbounds,bounds,met,space,ndivsm,ndivs,npts,path,unit)
2024 
2025 
2026 !This section has been created automatically by the script Abilint (TD).
2027 !Do not modify the following lines by hand.
2028 #undef ABI_FUNC
2029 #define ABI_FUNC 'make_path'
2030 !End of the abilint section
2031 
2032  implicit none
2033 
2034 !Arguments ------------------------------------
2035 !scalars
2036  integer,intent(in) :: nbounds,ndivsm
2037  integer,optional,intent(in) :: unit
2038  integer,intent(out) :: npts
2039  character(len=1),intent(in) :: space
2040 !arrays
2041  integer,intent(out) :: ndivs(nbounds-1)
2042  real(dp),intent(in) :: bounds(3,nbounds),met(3,3)
2043  real(dp),allocatable,intent(out) :: path(:,:)
2044 
2045 !Local variables-------------------------------
2046 !scalars
2047  integer,parameter :: prtvol=0
2048  integer :: idx,ii,jp,ount
2049  real(dp) :: nfact
2050  character(len=500) :: msg
2051 !arrays
2052  real(dp) :: diff(3),lng(nbounds-1)
2053 
2054 ! *************************************************************************
2055 
2056  ABI_CHECK(ndivsm > 0, sjoin('ndivsm', itoa(ndivsm)))
2057 
2058  ount = dev_null; if (present(unit)) ount = unit
2059 
2060  do ii=1,nbounds-1
2061    diff(:)=bounds(:,ii+1)-bounds(:,ii)
2062    lng(ii) = normv(diff,met,space)
2063  end do
2064 
2065  ! Avoid division by zero if any k(:,i+1)=k(:,i).
2066  nfact=MINVAL(lng)
2067  if (ABS(nfact)<tol6) then
2068    write(msg,'(3a)')&
2069 &    'Found two equivalent consecutive points in the path ',ch10,&
2070 &    'This is not allowed, modify the path in your input file'
2071    MSG_ERROR(msg)
2072  end if
2073 
2074  nfact=nfact/ndivsm
2075  ndivs(:)=NINT(lng(:)/nfact)
2076  npts=SUM(ndivs)+1 !1 for the first point
2077 
2078  write(msg,'(2a,i0,2a)')ch10,&
2079 &  ' Total number of points in the path: ',npts,ch10,&
2080 &  ' Number of divisions for each segment of the normalized path: '
2081  call wrtout(ount,msg,'COLL')
2082 
2083  do ii=1,nbounds-1
2084    write(msg,'(2(3f8.5,a),i0,a)')bounds(:,ii),' ==> ',bounds(:,ii+1),' ( ndivs : ',ndivs(ii),' )'
2085    call wrtout(ount,msg,'COLL')
2086  end do
2087  call wrtout(ount,ch10,'COLL')
2088 
2089  ! Allocate and construct the path.
2090  ABI_MALLOC(path,(3,npts))
2091 
2092  if (prtvol > 0) call wrtout(ount,' Normalized Path: ','COLL')
2093  idx=0
2094  do ii=1,nbounds-1
2095    do jp=1,ndivs(ii)
2096      idx=idx+1
2097      path(:,idx)=bounds(:,ii)+(jp-1)*(bounds(:,ii+1)-bounds(:,ii))/ndivs(ii)
2098      if (prtvol > 0) then
2099        write(msg,'(i4,4x,3(f8.5,1x))')idx,path(:,idx)
2100        call wrtout(ount,msg,'COLL')
2101      end if
2102    end do
2103  end do
2104  path(:,npts)=bounds(:,nbounds)
2105 
2106  if (prtvol > 0) then
2107    write(msg,'(i0,4x,3(f8.5,1x))')npts,path(:,npts)
2108    call wrtout(ount,msg,'COLL')
2109  end if
2110 
2111 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

PARENTS

      m_bz_mesh

CHILDREN

SOURCE

776 subroutine setup_k_rotation(nsym,timrev,symrec,nbz,kbz,gmet,krottb,krottbm1)
777 
778 
779 !This section has been created automatically by the script Abilint (TD).
780 !Do not modify the following lines by hand.
781 #undef ABI_FUNC
782 #define ABI_FUNC 'setup_k_rotation'
783 !End of the abilint section
784 
785  implicit none
786 
787 !Arguments ------------------------------------
788 !scalars
789  integer,intent(in) :: nbz,nsym,timrev
790 !arrays
791  integer,intent(in) :: symrec(3,3,nsym)
792  integer,intent(out) :: krottb(nbz,timrev,nsym),krottbm1(nbz,timrev,nsym)
793  real(dp),intent(in) :: kbz(3,nbz),gmet(3,3)
794 
795 !Local variables ------------------------------
796 !scalars
797  integer :: ik,ikp,isym,itim,nsh,ik_st !,sh_start
798  real(dp),parameter :: KTOL=tol6
799  real(dp) :: norm_old,norm !,norm_rot !norm_base
800  logical :: found,isok
801  character(len=500) :: msg
802 !arrays
803  integer :: g0(3)
804  integer :: iperm(nbz),shlim(nbz+1)
805  real(dp) :: shift(3),kbase(3),krot(3),knorm(nbz),kwrap(3),shlen(nbz+1)
806 
807 !************************************************************************
808 
809  DBG_ENTER("COLL")
810 
811  ! Sort the k-points according to their norm to speed up the search below.
812  do ik=1,nbz
813    call wrap2_pmhalf(kbz(1,ik),kwrap(1),shift(1))
814    call wrap2_pmhalf(kbz(2,ik),kwrap(2),shift(2))
815    call wrap2_pmhalf(kbz(3,ik),kwrap(3),shift(3))
816    knorm(ik) = normv(kwrap,gmet,"G")
817    iperm(ik)= ik
818  end do
819 
820  call sort_dp(nbz,knorm,iperm,KTOL)
821  !
822  ! The index of the initial (sorted) k-point in each shell
823  nsh=1; norm_old=knorm(1)
824 
825  shlim(1)=1; shlen(1) = norm_old
826  do ik_st=2,nbz
827    norm = knorm(ik_st)
828    if (ABS(norm-norm_old) > KTOL) then
829      norm_old   = norm
830      nsh        = nsh+1
831      shlim(nsh) = ik_st
832      shlen(nsh) = norm
833    end if
834  end do
835  shlim(nsh+1)=nbz+1
836  shlen(nsh+1)=HUGE(one)
837  !
838  ! === Set up k-rotation tables ===
839  ! * Use spatial inversion instead of time reversal whenever possible.
840  !call wrtout(std_out," Begin sorting ","COLL")
841 
842  isok=.TRUE.
843  do ik=1,nbz
844    kbase(:)=kbz(:,ik)
845    !
846    do itim=1,timrev
847      do isym=1,nsym
848        krot(:)=(3-2*itim)*MATMUL(symrec(:,:,isym),kbase)
849 
850        found=.FALSE.
851 #if 1
852       ! Old code
853        do ikp=1,nbz
854          if (isamek(krot,kbz(:,ikp),g0)) then
855            found=.TRUE.
856            krottb  (ik ,itim,isym)=ikp
857            krottbm1(ikp,itim,isym)=ik
858            EXIT
859          end if
860        end do
861 #else
862        !
863        ! Locate the shell index with bisection.
864        call wrap2_pmhalf(krot(1),kwrap(1),shift(1))
865        call wrap2_pmhalf(krot(2),kwrap(2),shift(2))
866        call wrap2_pmhalf(krot(3),kwrap(3),shift(3))
867        norm_rot = normv(kwrap,gmet,"G")
868        sh_start = bisect(shlen(1:nsh+1),norm_rot)
869 
870        do ik_st=shlim(sh_start),nbz
871          ikp = iperm(ik_st)
872          if (isamek(krot,kbz(:,ikp),g0)) then
873            found=.TRUE.
874            krottb  (ik ,itim,isym)=ikp
875            krottbm1(ikp,itim,isym)=ik
876            !write(std_out,*)ik_st,shlim(sh_start),nbz
877            EXIT
878          end if
879        end do
880 #endif
881        if (.not.found) then
882          isok=.FALSE.
883          !write(std_out,*)" norm_base,norm_rot ",norm_base,norm_rot
884          !write(std_out,*)normv(kbase,gmet,"G"),normv(krot,gmet,"G")
885          write(msg,'(2(a,i4),2x,2(3f12.6,2a),i3,a,i2)')&
886 &          'Initial k-point ',ik,'/',nbz,kbase(:),ch10,&
887 &          'Rotated k-point (not found) ',krot(:),ch10,&
888 &          'Through symmetry operation ',isym,' and itim ',itim
889          MSG_ERROR(msg)
890        end if
891 
892      end do
893    end do
894  end do
895 
896  if (.not.isok) then
897    MSG_ERROR('k-mesh not closed')
898  end if
899 
900  DBG_EXIT("COLL")
901 
902 end subroutine setup_k_rotation