TABLE OF CONTENTS


ABINIT/m_kpts [ Modules ]

[ Top ] [ Modules ]

NAME

  m_kpts

FUNCTION

COPYRIGHT

 Copyright (C) 2008-2018 ABINIT group (XG, MG, MJV, DRH, DCA, JCC, MM)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_kpts
27 
28  use defs_basis
29  use m_errors
30  use m_abicore
31  use m_crystal
32  use m_sort
33  use m_kptrank
34  use m_xmpi
35 
36  use m_time,           only : timab
37  use m_symtk,          only : mati3inv, mati3det, matr3inv, smallprim
38  use m_fstrings,       only : sjoin, itoa, ltoa
39  use m_numeric_tools,  only : wrap2_pmhalf
40  use m_geometry,       only : metric
41  use m_tetrahedron,    only : t_tetrahedron, init_tetra, destroy_tetra
42  use m_symkpt,     only : symkpt
43 
44  implicit none
45 
46  private
47 
48  public :: kpts_timrev_from_kptopt   ! Returns the value of timrev from kptopt
49  public :: kpts_ibz_from_kptrlatt    ! Determines the IBZ, the weights and the BZ from kptrlatt
50  public :: tetra_from_kptrlatt       ! Create an instance of `t_tetrahedron` from kptrlatt and shiftk
51  public :: symkchk                   ! Checks that the set of k points has the full space group symmetry,
52                                      ! modulo time reversal if appropriate.
53  public :: listkk                    ! Find correspondence between two set of k-points.
54  public :: getkgrid                  ! Compute the grid of k points in the irreducible Brillouin zone.
55  !FIXME: Deprecated
56  public :: get_full_kgrid            ! Create full grid of kpoints and find equivalent irred ones.
57                                      ! Duplicates work in getkgrid, but need all outputs of kpt_fullbz, and indkpt
58 
59  ! FIXME: Deprecated. Used in invars2 but call is commented
60  !public :: get_kpt_fullbz            ! Create full grid of kpoints from kptrlatt and shiftk
61 
62  public :: smpbz                      ! Generate a set of special k (or q) points which samples in a homogeneous way the BZ
63  public :: testkgrid                  ! Test different grids of k points.
64 
65  ! FIXME: deprecated
66  public :: mknormpath

m_kpts/get_full_kgrid [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 get_full_kgrid

FUNCTION

 Create full grid of kpoints and find equivalent
 irred ones. Duplicates work in getkgrid, but need all outputs of kpt_fullbz, and indkpt

INPUTS

  kpt(3,nkpt)=irreducible kpoints
  kptrlatt(3,3)=lattice vectors for full kpoint grid
  nkpt=number of irreducible kpoints
  nkpt_fullbz=number of kpoints in full brillouin zone
  nshiftk=number of kpoint grid shifts
  nsym=number of symmetries
  shiftk(3,nshiftk)=kpoint shifts
  symrel(3,3,nsym)=symmetry matrices in real space

OUTPUT

  indkpt(nkpt_fullbz)=non-symmetrized indices of the k-points (see symkpt.f)
  kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone

NOTES

  MG: The present inplementation always assumes kptopt==1 !!!!

 TODO: This routine should be removed

PARENTS

      m_phonons

CHILDREN

      destroy_kptrank,get_kpt_fullbz,get_rank_1kpt,mati3inv,mkkptrank

SOURCE

1435 subroutine get_full_kgrid(indkpt,kpt,kpt_fullbz,kptrlatt,nkpt,&
1436 & nkpt_fullbz,nshiftk,nsym,shiftk,symrel)
1437 
1438 
1439 !This section has been created automatically by the script Abilint (TD).
1440 !Do not modify the following lines by hand.
1441 #undef ABI_FUNC
1442 #define ABI_FUNC 'get_full_kgrid'
1443 !End of the abilint section
1444 
1445  implicit none
1446 
1447 !Arguments ------------------------------------
1448 !scalars
1449  integer,intent(in) :: nkpt,nkpt_fullbz,nshiftk,nsym
1450 !arrays
1451  integer,intent(in) :: kptrlatt(3,3),symrel(3,3,nsym)
1452  integer,intent(out) :: indkpt(nkpt_fullbz)
1453  real(dp),intent(in) :: kpt(3,nkpt),shiftk(3,nshiftk)
1454  real(dp),intent(out) :: kpt_fullbz(3,nkpt_fullbz)
1455 
1456 !Local variables-------------------------------
1457 !scalars
1458  integer :: ikpt,isym,itim,timrev
1459  integer :: symrankkpt
1460  character(len=500) :: message
1461  type(kptrank_type) :: kptrank_t
1462 
1463 !arrays
1464  integer :: inv_symrel(3,3,nsym)
1465  real(dp) :: k2(3)
1466 
1467 ! *********************************************************************
1468 
1469 !Invert symrels => gives symrels for kpoints
1470 
1471  do isym=1,nsym
1472    call mati3inv (symrel(:,:,isym),inv_symrel(:,:,isym))
1473  end do
1474 
1475  call get_kpt_fullbz(kpt_fullbz,kptrlatt,nkpt_fullbz,nshiftk,shiftk)
1476 
1477 !make full k-point rank arrays
1478  call mkkptrank (kpt,nkpt,kptrank_t)
1479 
1480 !
1481 !find equivalence to irred kpoints in kpt
1482 !
1483  indkpt(:) = 0
1484  timrev=1 ! includes the time inversion symmetry
1485  do ikpt=1,nkpt_fullbz
1486    do isym=1,nsym
1487      do itim=1,(1-2*timrev),-2
1488 
1489        k2(:) = itim*(inv_symrel(:,1,isym)*kpt_fullbz(1,ikpt) + &
1490 &       inv_symrel(:,2,isym)*kpt_fullbz(2,ikpt) + &
1491 &       inv_symrel(:,3,isym)*kpt_fullbz(3,ikpt))
1492 
1493        call get_rank_1kpt (k2,symrankkpt,kptrank_t)
1494        if (kptrank_t%invrank(symrankkpt) /= -1) indkpt(ikpt) = kptrank_t%invrank(symrankkpt)
1495 
1496      end do ! loop time reversal symmetry
1497    end do !  loop sym ops
1498 
1499    if (indkpt(ikpt) == 0) then
1500      write (message,'(a,i0)')' indkpt(ikpt) is still 0: no irred kpoint is equiv to ikpt ',ikpt
1501      MSG_BUG(message)
1502    end if
1503  end do !  loop full kpts
1504 
1505  call destroy_kptrank (kptrank_t)
1506 
1507 end subroutine get_full_kgrid

m_kpts/get_kpt_fullbz [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 get_kpt_fullbz

FUNCTION

 Create full grid of kpoints from kptrlatt and shiftk

INPUTS

  kptrlatt(3,3)=lattice vectors for full kpoint grid
  nkpt_fullbz=number of kpoints in full brillouin zone
  nshiftk=number of kpoint grid shifts
  shiftk(3,nshiftk)=kpoint shifts

OUTPUT

  kpt_fullbz(3,nkpt_fullbz)=kpoints in full brillouin zone

PARENTS

      get_full_kgrid,invars2

CHILDREN

      mati3det,matr3inv,wrap2_pmhalf

SOURCE

1534 subroutine get_kpt_fullbz(kpt_fullbz,kptrlatt,nkpt_fullbz,nshiftk,shiftk)
1535 
1536 
1537 !This section has been created automatically by the script Abilint (TD).
1538 !Do not modify the following lines by hand.
1539 #undef ABI_FUNC
1540 #define ABI_FUNC 'get_kpt_fullbz'
1541 !End of the abilint section
1542 
1543  implicit none
1544 
1545 !Arguments ------------------------------------
1546 !scalars
1547  integer,intent(in) :: nkpt_fullbz,nshiftk
1548 !arrays
1549  integer,intent(in) :: kptrlatt(3,3)
1550  real(dp),intent(in) :: shiftk(3,nshiftk)
1551  real(dp),intent(out) :: kpt_fullbz(3,nkpt_fullbz)
1552 
1553 !Local variables-------------------------------
1554 !scalars
1555  integer, parameter :: max_number_of_prime=47,mshiftk=210
1556  integer :: det,ii,ikshft,iprim,jj,kk,nn
1557  character(len=500) :: message
1558 
1559 !arrays
1560  integer :: boundmax(3),boundmin(3),common_factor(3)
1561  integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,&
1562 &  29,31,37,41,43, 47,53,59,61,67,&
1563 &  71,73,79,83,89, 97,101,103,107,109,&
1564 &  113,127,131,137,139, 149,151,157,163,167,&
1565 &  173,179,181,191,193, 197,199/)
1566  real(dp) :: k1(3),k2(3),klatt(3,3),rlatt(3,3),shift(3),test_rlatt(3,3)
1567 
1568 ! *********************************************************************
1569 
1570 !Identify first factors that can be used to rescale the three kptrlatt vectors
1571 !Only test a large set of prime factors, though ...
1572  do jj=1,3
1573    common_factor(jj)=1
1574    rlatt(:,jj)=kptrlatt(:,jj)
1575    do iprim=1,max_number_of_prime
1576      test_rlatt(:,jj)=rlatt(:,jj)/dble(prime_factor(iprim))
1577 !    If one of the components is lower than 1 in absolute value, then it is not worth to continue the search.
1578      if(minval(abs(abs(test_rlatt(:,jj))-half))<half-tol8)exit
1579      do
1580        if(sum(abs(test_rlatt(:,jj)-nint(test_rlatt(:,jj)) ))<tol8)then
1581          common_factor(jj)=prime_factor(iprim)*common_factor(jj)
1582          rlatt(:,jj)=rlatt(:,jj)/dble(prime_factor(iprim))
1583          test_rlatt(:,jj)=test_rlatt(:,jj)/dble(prime_factor(iprim))
1584        else
1585          exit
1586        end if
1587      end do
1588    end do
1589  end do
1590  call mati3det(kptrlatt,det)
1591  det=det/(common_factor(1)*common_factor(2)*common_factor(3))
1592 
1593  rlatt(:,:)=kptrlatt(:,:)
1594  call matr3inv(rlatt,klatt)
1595 !Now, klatt contains the three primitive vectors of the k lattice,
1596 !in reduced coordinates. One builds all k vectors that
1597 !are contained in the first Brillouin zone, with coordinates
1598 !in the interval [0,1[ . First generate boundaries of a big box.
1599 !In order to generate all possible vectors in the reciprocal space,
1600 !one must consider all multiples of the primitive ones, until a vector with only integers is found.
1601 !The maximum bound is the scale of the corresponding kptrlatt vector, times the determinant of kptrlatt. Also consider negative vectors.
1602 !On this basis, compute the bounds.
1603  do jj=1,3
1604 !  To accomodate the shifts, boundmin starts from -1
1605 !  Well, this is not a complete solution ...
1606    boundmin(jj)=-1-common_factor(jj)*abs(det)
1607    boundmax(jj)=common_factor(jj)*abs(det)
1608  end do
1609 
1610  nn=1
1611  do kk=boundmin(3),boundmax(3)
1612    do jj=boundmin(2),boundmax(2)
1613      do ii=boundmin(1),boundmax(1)
1614        do ikshft=1,nshiftk
1615 
1616 !        Coordinates of the trial k point with respect to the k primitive lattice
1617          k1(1)=ii+shiftk(1,ikshft)
1618          k1(2)=jj+shiftk(2,ikshft)
1619          k1(3)=kk+shiftk(3,ikshft)
1620 
1621 !        Reduced coordinates of the trial k point
1622          k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3)
1623 
1624 !        Eliminate the point if outside [0,1[
1625          if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle
1626          if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle
1627          if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle
1628 
1629 !        Wrap the trial values in the interval ]-1/2,1/2] .
1630          call wrap2_pmhalf(k2(1),k1(1),shift(1))
1631          call wrap2_pmhalf(k2(2),k1(2),shift(2))
1632          call wrap2_pmhalf(k2(3),k1(3),shift(3))
1633          if(nn > nkpt_fullbz) then
1634            write (message,'(a,i0)')' nkpt_fullbz mis-estimated, exceed nn=',nn
1635            MSG_BUG(message)
1636          end if
1637          kpt_fullbz(:,nn)=k1(:)
1638          nn=nn+1
1639        end do
1640      end do
1641    end do
1642  end do
1643  nn = nn-1
1644 
1645  if (nn /= nkpt_fullbz) then
1646    write (message,'(2(a,i0),a,a)')' nkpt_fullbz= ',nkpt_fullbz,' underestimated  nn=',nn,&
1647 &   ch10, "Perhaps your k grid or shifts do not correspond to the symmetry?"
1648    MSG_BUG(message)
1649  end if
1650 
1651 end subroutine get_kpt_fullbz

m_kpts/getkgrid [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 getkgrid

FUNCTION

 Compute the grid of k points in the irreducible Brillouin zone.
 Note that nkpt (and nkpthf) can be computed by calling this routine with nkpt=0, provided that kptopt/=0.
 If downsampling is present, also compute a downsampled k grid.

INPUTS

 chksymbreak= if 1, will check whether the k point grid is symmetric (for kptopt=1,2 and 4), and stop if not.
 iout=unit number for echoed output . 0 if no output is wished.
 iscf= ( <= 0 =>non-SCF), >0 => SCF)  MG: FIXME I don't understand why we have to pass the value iscf.
 kptopt=option for the generation of k points
   (defines whether spatial symmetries and/or time-reversal can be used)
 msym=default maximal number of symmetries
 nsym=number of symmetries
 rprimd(3,3)=dimensional real space primitive translations (bohr)
 symafm(nsym)=(anti)ferromagnetic part of symmetry operations
 symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
 vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
 [downsampling(3) = input variable that governs the downsampling]

OUTPUT

 kptrlen=length of the smallest real space supercell vector associated with the lattice of k points.
 nkpt_computed=number of k-points in the IBZ computed in the present routine
 If nkpt/=0  the following are also output :
   kpt(3,nkpt)=reduced coordinates of k points.
   wtk(nkpt)=weight assigned to each k point.
 [fullbz(3,nkpt_fullbz)]=k-points generated in the full Brillouin zone.
   In output: allocated array with the list of k-points in the BZ.
 [kpthf(3,nkpthf)]=k-points generated in the full Brillouin zone, possibly downsampled (for Fock).

NOTES

  msym not needed since nsym is the last index.

SIDE EFFECTS

 Input/Output
 nkpt=number of k points (might be zero, see output description)
 kptrlatt(3,3)=k-point lattice specification
 nshiftk=actual number of k-point shifts in shiftk
 shiftk(3,210)=shift vectors for k point generation
 [nkpthf = number of k points in the full BZ, for the Fock operator]

PARENTS

      ep_setupqpt,getshell,inkpts,inqpt,m_ab7_kpoints,m_bz_mesh,m_kpts
      nonlinear,testkgrid,thmeig

CHILDREN

      mati3inv,matr3inv,metric,smallprim,smpbz,symkpt

SOURCE

 968 subroutine getkgrid(chksymbreak,iout,iscf,kpt,kptopt,kptrlatt,kptrlen,&
 969 & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk,&
 970 & fullbz,nkpthf,kpthf,downsampling) ! optional
 971 
 972 
 973 !This section has been created automatically by the script Abilint (TD).
 974 !Do not modify the following lines by hand.
 975 #undef ABI_FUNC
 976 #define ABI_FUNC 'getkgrid'
 977 !End of the abilint section
 978 
 979  implicit none
 980 
 981 !Arguments ------------------------------------
 982 !scalars
 983  integer,intent(in) :: chksymbreak,iout,iscf,kptopt,msym,nkpt,nsym
 984  integer,intent(inout),optional :: nkpthf
 985  integer,intent(inout) :: nshiftk
 986  integer,intent(inout) :: nkpt_computed !vz_i
 987  real(dp),intent(out) :: kptrlen
 988 !arrays
 989  integer,intent(in) :: symafm(msym),symrel(3,3,msym),vacuum(3)
 990  integer,optional,intent(in) :: downsampling(3)
 991  integer,intent(inout) :: kptrlatt(3,3)
 992  real(dp),intent(in) :: rprimd(3,3)
 993  real(dp),intent(inout) :: shiftk(3,210)
 994  real(dp),intent(inout) :: kpt(3,nkpt) !vz_i
 995  real(dp),intent(inout) :: wtk(nkpt)
 996  real(dp),optional,allocatable,intent(out) :: fullbz(:,:)
 997  real(dp),optional,intent(out) :: kpthf(:,:)
 998 
 999 !Local variables-------------------------------
1000 !scalars
1001  integer, parameter :: max_number_of_prime=47,mshiftk=210
1002  integer :: brav,decreased,found,ii,ikpt,iprime,ishiftk,isym,jshiftk,kshiftk,mkpt,mult
1003  integer :: nkpthf_computed,nkpt_fullbz,nkptlatt,nshiftk2,nsym_used,option
1004  integer :: test_prime,timrev
1005  real(dp) :: length2,ucvol,ucvol_super
1006  character(len=500) :: message
1007 !arrays
1008  integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,&
1009 &  29,31,37,41,43, 47,53,59,61,67,&
1010 &  71,73,79,83,89, 97,101,103,107,109,&
1011 &  113,127,131,137,139, 149,151,157,163,167,&
1012 &  173,179,181,191,193, 197,199/)
1013  integer :: kptrlatt2(3,3)
1014  integer,allocatable :: belong_chain(:),generator(:),indkpt(:),number_in_chain(:)
1015  integer,allocatable :: repetition_factor(:),symrec(:,:,:)
1016 ! real(dp) :: cart(3,3)
1017  real(dp) :: dijk(3),delta_dmult(3),dmult(3),fact_vacuum(3),gmet(3,3)
1018  real(dp) :: gmet_super(3,3),gprimd(3,3),gprimd_super(3,3),klatt2(3,3)
1019  real(dp) :: klatt3(3,3),kptrlattr(3,3),ktransf(3,3),ktransf_invt(3,3)
1020  real(dp) :: metmin(3,3),minim(3,3),rmet(3,3),rmet_super(3,3),rprimd_super(3,3)
1021  real(dp),allocatable :: deltak(:,:),kpt_fullbz(:,:),shiftk2(:,:),shiftk3(:,:),spkpt(:,:),wtk_folded(:),wtk_fullbz(:)
1022 
1023 ! *************************************************************************
1024 
1025  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
1026 
1027  if (kptopt==1.or.kptopt==4) then
1028    ! Cannot use antiferromagnetic symmetry operations to decrease the number of k points
1029    nsym_used=0
1030    do isym=1,nsym
1031      if(symafm(isym)==1)nsym_used=nsym_used+1
1032    end do
1033    ABI_ALLOCATE(symrec,(3,3,nsym_used))
1034    nsym_used=0
1035    do isym=1,nsym ! Get the symmetry matrices in terms of reciprocal basis
1036      if(symafm(isym)==1)then
1037        nsym_used=nsym_used+1
1038        call mati3inv(symrel(:,:,isym),symrec(:,:,nsym_used))
1039      end if
1040    end do
1041  else if (kptopt==2) then
1042    !Use only the time-reversal
1043    nsym_used=1
1044    ABI_ALLOCATE(symrec,(3,3,1))
1045    symrec(1:3,1:3,1)=0
1046    do ii=1,3
1047      symrec(ii,ii,1)=1
1048    end do
1049  end if
1050 
1051  kptrlatt2(:,:)=kptrlatt(:,:)
1052  nshiftk2=nshiftk
1053  ABI_ALLOCATE(shiftk2,(3,mshiftk))
1054  ABI_ALLOCATE(shiftk3,(3,mshiftk))
1055  shiftk2(:,:)=shiftk(:,:)
1056 
1057 !Find a primitive k point lattice, if possible, by decreasing the number of shifts.
1058  if(nshiftk2/=1)then
1059 
1060    do ! Loop to be repeated if there has been a successful reduction of nshiftk2
1061 
1062      ABI_ALLOCATE(deltak,(3,nshiftk2))
1063      ABI_ALLOCATE(repetition_factor,(nshiftk2))
1064      ABI_ALLOCATE(generator,(nshiftk2))
1065      ABI_ALLOCATE(belong_chain,(nshiftk2))
1066      ABI_ALLOCATE(number_in_chain,(nshiftk2))
1067 
1068      decreased=0
1069      deltak(1,1:nshiftk2)=shiftk2(1,1:nshiftk2)-shiftk2(1,1)
1070      deltak(2,1:nshiftk2)=shiftk2(2,1:nshiftk2)-shiftk2(2,1)
1071      deltak(3,1:nshiftk2)=shiftk2(3,1:nshiftk2)-shiftk2(3,1)
1072      deltak(:,:)=deltak(:,:)-floor(deltak(:,:)+tol8)
1073 
1074 !    Identify for each shift, the smallest repetition prime factor that yields a reciprocal lattice vector.
1075      repetition_factor(:)=0
1076      repetition_factor(1)=1
1077      do ishiftk=2,nshiftk2
1078        do iprime=1,max_number_of_prime
1079          test_prime=prime_factor(iprime)
1080          dmult(:)=test_prime*deltak(:,ishiftk)
1081          if(sum(abs( dmult(:)-nint(dmult(:)) ))<tol8)then
1082            repetition_factor(ishiftk)=test_prime
1083            exit
1084          end if
1085        end do
1086      end do
1087 
1088 !    Initialize the selection of tentative generators
1089      generator(:)=1
1090      do ishiftk=1,nshiftk2
1091        if(repetition_factor(ishiftk)==0 .or. repetition_factor(ishiftk)==1)generator(ishiftk)=0
1092      end do
1093 
1094 !    Try different shifts as generators, by order of increasing repetition factor,
1095 !    provided they are equal or bigger than 2
1096      do iprime=1,max_number_of_prime
1097        do ishiftk=2,nshiftk2
1098          ! Note that ishiftk=1 is never a generator. It is the reference starting point.
1099          if(generator(ishiftk)==1 .and. repetition_factor(ishiftk)==prime_factor(iprime))then
1100 !          Test the generator : is it indeed closed ?
1101            if(prime_factor(iprime)/=2)then
1102              do mult=2,prime_factor(iprime)-1
1103                dmult(:)=mult*deltak(:,ishiftk)
1104                found=0
1105                do jshiftk=1,nshiftk2
1106                  delta_dmult(:)=deltak(:,jshiftk)-dmult(:)
1107                  if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then
1108                    found=1
1109                    exit
1110                  end if
1111                end do
1112                if(found==0)exit
1113              end do
1114              if(found==0)generator(ishiftk)=0
1115            end if
1116            if(generator(ishiftk)==0)cycle
1117          else
1118            cycle
1119          end if
1120 !        Now, test whether all k points can be found in all possible chains
1121          belong_chain(:)=0
1122          do jshiftk=1,nshiftk2
1123 !          Initialize a chain starting from a k point not yet in a chain
1124            if(belong_chain(jshiftk)==0)then
1125              number_in_chain(:)=0   ! Not a member of the chain (yet)
1126              number_in_chain(jshiftk)=1   ! The first point in chain
1127              do mult=1,prime_factor(iprime)-1
1128                dmult(:)=mult*deltak(:,ishiftk)
1129                found=0
1130                do kshiftk=jshiftk+1,nshiftk2
1131                  delta_dmult(:)=deltak(:,kshiftk)-deltak(:,jshiftk)-dmult(:)
1132                  if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then
1133                    found=1
1134                    number_in_chain(kshiftk)=mult+1
1135                    exit
1136                  end if
1137                end do
1138                if(found==0)then
1139                  generator(ishiftk)=0
1140                  exit
1141                end if
1142              end do
1143              if(generator(ishiftk)==1)then
1144 !              Store the chain
1145                do kshiftk=1,nshiftk2
1146                  if(number_in_chain(kshiftk)/=0)belong_chain(kshiftk)=number_in_chain(kshiftk)
1147                end do
1148              else
1149                exit
1150              end if
1151            end if
1152          end do
1153 
1154          if(generator(ishiftk)==0)cycle
1155 
1156 !        For the generator based on ishiftk, all the k points have been found to belong to one chain.
1157 !        All the initializing k points in the different chains have belong_chain(:)=1 .
1158 !        They must be kept, and the others thrown away.
1159          ktransf(:,:)=0.0_dp
1160          ktransf(1,1)=1.0_dp
1161          ktransf(2,2)=1.0_dp
1162          ktransf(3,3)=1.0_dp
1163 !        Replace one of the unit vectors by the shift vector deltak(:,ishiftk).
1164 !        However, must pay attention not to make linear combinations.
1165 !        Also, choose positive sign for first-non-zero value.
1166          if(abs(deltak(1,ishiftk)-nint(deltak(1,ishiftk)))>tol8)then
1167            if(deltak(1,ishiftk)>0)ktransf(:,1)= deltak(:,ishiftk)
1168            if(deltak(1,ishiftk)<0)ktransf(:,1)=-deltak(:,ishiftk)
1169          else if(abs(deltak(2,ishiftk)-nint(deltak(2,ishiftk)))>tol8)then
1170            if(deltak(2,ishiftk)>0)ktransf(:,2)= deltak(:,ishiftk)
1171            if(deltak(2,ishiftk)<0)ktransf(:,2)=-deltak(:,ishiftk)
1172          else if(abs(deltak(3,ishiftk)-nint(deltak(3,ishiftk)))>tol8)then
1173            if(deltak(3,ishiftk)>0)ktransf(:,3)= deltak(:,ishiftk)
1174            if(deltak(3,ishiftk)<0)ktransf(:,3)=-deltak(:,ishiftk)
1175          end if
1176 !        Copy the integers to real(dp)
1177          kptrlattr(:,:)=kptrlatt2(:,:)
1178 !        Go to reciprocal space
1179          call matr3inv(kptrlattr,klatt2)
1180 !        Make the transformation
1181          do ii=1,3
1182            klatt3(:,ii)=ktransf(1,ii)*klatt2(:,1)+ktransf(2,ii)*klatt2(:,2)+ktransf(3,ii)*klatt2(:,3)
1183          end do
1184 !        Back to real space
1185          call matr3inv(klatt3,kptrlattr)
1186 !        real(dp) to integer
1187          kptrlatt2(:,:)=nint(kptrlattr(:,:))
1188 !        Prepare the transformation of the shifts
1189          call matr3inv(ktransf,ktransf_invt)
1190          decreased=1
1191          kshiftk=0
1192          do jshiftk=1,nshiftk2
1193            if(belong_chain(jshiftk)==1)then
1194              kshiftk=kshiftk+1
1195 !            Place the shift with index jshiftk in place of the one in kshiftk,
1196 !            also transform the shift from the old to the new coordinate system
1197              shiftk3(:,kshiftk)=ktransf_invt(1,:)*shiftk2(1,jshiftk)+&
1198 &             ktransf_invt(2,:)*shiftk2(2,jshiftk)+&
1199 &             ktransf_invt(3,:)*shiftk2(3,jshiftk)
1200            end if
1201          end do
1202          nshiftk2=nshiftk2/prime_factor(iprime)
1203          shiftk2(:,1:nshiftk2)=shiftk3(:,1:nshiftk2)-floor(shiftk3(:,1:nshiftk2)+tol8)
1204          if(kshiftk/=nshiftk2)then
1205            MSG_BUG('The search for a primitive k point lattice contains a bug.')
1206          end if
1207 
1208 !        If this trial shift was successful, must exit the loop on trial ishiftk,
1209 !        and reinitialize the global loop
1210          if(decreased==1)exit
1211        end do ! ishiftk
1212        if(decreased==1)exit
1213      end do ! iprime
1214 
1215      ABI_DEALLOCATE(belong_chain)
1216      ABI_DEALLOCATE(deltak)
1217      ABI_DEALLOCATE(number_in_chain)
1218      ABI_DEALLOCATE(repetition_factor)
1219      ABI_DEALLOCATE(generator)
1220 
1221      if(decreased==0 .or. nshiftk2==1)exit
1222 
1223    end do ! Infinite loop
1224 
1225  end if !  End nshiftk being 1 or larger
1226 
1227 !Impose shiftk coordinates to be in [0,1[
1228  do ishiftk=1,nshiftk2
1229    do ii=1,3
1230      if(shiftk2(ii,ishiftk)>one-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)-1.0_dp
1231      if(shiftk2(ii,ishiftk)<-tol8)    shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)+1.0_dp
1232    end do
1233  end do
1234 
1235 !Compute the number of k points in the G-space unit cell
1236  nkptlatt=kptrlatt2(1,1)*kptrlatt2(2,2)*kptrlatt2(3,3) &
1237 & +kptrlatt2(1,2)*kptrlatt2(2,3)*kptrlatt2(3,1) &
1238 & +kptrlatt2(1,3)*kptrlatt2(2,1)*kptrlatt2(3,2) &
1239 & -kptrlatt2(1,2)*kptrlatt2(2,1)*kptrlatt2(3,3) &
1240 & -kptrlatt2(1,3)*kptrlatt2(2,2)*kptrlatt2(3,1) &
1241 & -kptrlatt2(1,1)*kptrlatt2(2,3)*kptrlatt2(3,2)
1242 
1243 !Check whether the number of k points is positive,
1244 !otherwise, change the handedness of kptrlatt2
1245  if(nkptlatt<=0)then
1246 !  write(std_out,*)' getkgrid : nkptlatt is negative !'
1247    kptrlatt2(:,3)=-kptrlatt2(:,3)
1248    nkptlatt=-nkptlatt
1249    do ishiftk=1,nshiftk2
1250      shiftk2(3,ishiftk)=-shiftk2(3,ishiftk)
1251    end do
1252  end if
1253 
1254 !Determine the smallest supercell R-vector whose contribution
1255 !is not taken correctly into account in the k point integration.
1256 !Increase enormously the size of the cell when vacuum is present.
1257  fact_vacuum(:)=1
1258  if(vacuum(1)==1)fact_vacuum(1)=1000.0_dp
1259  if(vacuum(2)==1)fact_vacuum(2)=1000.0_dp
1260  if(vacuum(3)==1)fact_vacuum(3)=1000.0_dp
1261  do ii=1,3
1262    rprimd_super(:,ii)=fact_vacuum(1)*rprimd(:,1)*kptrlatt2(1,ii)+&
1263 &   fact_vacuum(2)*rprimd(:,2)*kptrlatt2(2,ii)+&
1264 &   fact_vacuum(3)*rprimd(:,3)*kptrlatt2(3,ii)
1265  end do
1266 
1267  call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super)
1268  call smallprim(metmin,minim,rprimd_super)
1269  length2=min(metmin(1,1),metmin(2,2),metmin(3,3))
1270  kptrlen=sqrt(length2)
1271 
1272  !write(message,'(a,es16.6)' )' getkgrid : length of smallest supercell vector (bohr)=',kptrlen
1273  !call wrtout(std_out,message,'COLL')
1274 ! If the number of shifts has been decreased, determine the set of kptrlatt2 vectors
1275 ! with minimal length (without using fact_vacuum)
1276 ! It is worth to determine the minimal set of vectors so that the kptrlatt that is output
1277 ! does not seem screwy, although correct but surprising.
1278  if(nshiftk/=nshiftk2)then
1279    do ii=1,3
1280      rprimd_super(:,ii)=rprimd(:,1)*kptrlatt2(1,ii)+rprimd(:,2)*kptrlatt2(2,ii)+rprimd(:,3)*kptrlatt2(3,ii)
1281    end do
1282    call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super)
1283 !  Shift vectors in cartesian coordinates (reciprocal space)
1284    do ishiftk=1,nshiftk2
1285      shiftk3(:,ishiftk)=gprimd_super(:,1)*shiftk2(1,ishiftk)+&
1286 &     gprimd_super(:,2)*shiftk2(2,ishiftk)+&
1287 &     gprimd_super(:,3)*shiftk2(3,ishiftk)
1288    end do
1289    call smallprim(metmin,minim,rprimd_super)
1290    call metric(gmet_super,gprimd_super,-1,rmet_super,minim,ucvol_super)
1291 !  This is the new kptrlatt2
1292    do ii=1,3
1293      dijk(:)=gprimd(1,:)*minim(1,ii)+&
1294 &     gprimd(2,:)*minim(2,ii)+&
1295 &     gprimd(3,:)*minim(3,ii)
1296      kptrlatt2(:,ii)=nint(dijk(:))
1297    end do
1298 !  Shifts in the new set of kptrlatt vectors
1299    do ishiftk=1,nshiftk2
1300      shiftk2(:,ishiftk)=minim(1,:)*shiftk3(1,ishiftk)+&
1301 &     minim(2,:)*shiftk3(2,ishiftk)+&
1302 &     minim(3,:)*shiftk3(3,ishiftk)
1303    end do
1304  end if
1305 
1306 !brav=1 is able to treat all bravais lattices.
1307  brav=1
1308  mkpt=nkptlatt*nshiftk2
1309 
1310  ABI_ALLOCATE(spkpt,(3,mkpt))
1311  option=0
1312  if(iout/=0)option=1
1313 
1314  if (PRESENT(downsampling))then
1315    call smpbz(brav,iout,kptrlatt2,mkpt,nkpthf_computed,nshiftk2,option,shiftk2,spkpt,downsampling=downsampling)
1316    if (PRESENT(kpthf) .and. nkpthf/=0) then ! Returns list of k-points in the Full BZ, possibly downsampled for Fock
1317      kpthf = spkpt(:,1:nkpthf)
1318    end if
1319    nkpthf=nkpthf_computed
1320 
1321  end if
1322 
1323  call smpbz(brav,iout,kptrlatt2,mkpt,nkpt_fullbz,nshiftk2,option,shiftk2,spkpt)
1324 
1325  if (PRESENT(fullbz)) then ! Returns list of k-points in the Full BZ.
1326    ABI_ALLOCATE(fullbz,(3,nkpt_fullbz))
1327    fullbz = spkpt(:,1:nkpt_fullbz)
1328  end if
1329 
1330  if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then
1331 
1332    ABI_ALLOCATE(indkpt,(nkpt_fullbz))
1333    ABI_ALLOCATE(kpt_fullbz,(3,nkpt_fullbz))
1334    ABI_ALLOCATE(wtk_fullbz,(nkpt_fullbz))
1335    ABI_ALLOCATE(wtk_folded,(nkpt_fullbz))
1336 
1337    kpt_fullbz(:,:)=spkpt(:,1:nkpt_fullbz)
1338    wtk_fullbz(1:nkpt_fullbz)=1.0_dp/dble(nkpt_fullbz)
1339 
1340    timrev=1;if (kptopt==4) timrev=0
1341 
1342    call symkpt(chksymbreak,gmet,indkpt,iout,kpt_fullbz,nkpt_fullbz,&
1343 &   nkpt_computed,nsym_used,symrec,timrev,wtk_fullbz,wtk_folded)
1344 
1345    ABI_DEALLOCATE(symrec)
1346    ABI_DEALLOCATE(wtk_fullbz)
1347 
1348  else if(kptopt==3)then
1349    nkpt_computed=nkpt_fullbz
1350  end if
1351 
1352 !The number of k points has been computed from kptopt, kptrlatt, nshiftk, shiftk,
1353 !and the eventual symmetries, it is presently called nkpt_computed.
1354 
1355 !Check that the argument nkpt is coherent with nkpt_computed, if nkpt/=0.
1356  if(nkpt/=nkpt_computed .and. nkpt/=0)then
1357    write(message, '(a,i6,5a,i6,7a)') &
1358 &   'The argument nkpt=',nkpt,', does not match',ch10,&
1359 &   'the number of k points generated by kptopt, kptrlatt, shiftk,',ch10,&
1360 &   'and the eventual symmetries, that is, nkpt=',nkpt_computed,'.',ch10,&
1361 &   'However, note that it might be due to the user,',ch10,&
1362 &   'if nkpt is explicitely defined in the input file.',ch10,&
1363 &   'In this case, please check your input file.'
1364    MSG_BUG(message)
1365  end if
1366 
1367  if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then
1368 
1369    if(nkpt/=0)then
1370      do ikpt=1,nkpt
1371        kpt(:,ikpt)=kpt_fullbz(:,indkpt(ikpt))
1372        if(iscf>=0 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(ikpt)=wtk_folded(indkpt(ikpt))
1373      end do
1374    end if
1375 
1376    ABI_DEALLOCATE(indkpt)
1377    ABI_DEALLOCATE(kpt_fullbz)
1378    ABI_DEALLOCATE(spkpt)
1379    ABI_DEALLOCATE(wtk_folded)
1380 
1381  else if(kptopt==3)then
1382 
1383    if(nkpt/=0)then
1384      kpt(:,1:nkpt)=spkpt(:,1:nkpt)
1385      if(iscf>1 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(1:nkpt)=1.0_dp/dble(nkpt)
1386    end if
1387    ABI_DEALLOCATE(spkpt)
1388 
1389  end if
1390 
1391  kptrlatt(:,:)=kptrlatt2(:,:)
1392  nshiftk=nshiftk2
1393  shiftk(:,1:nshiftk)=shiftk2(:,1:nshiftk)
1394  ABI_DEALLOCATE(shiftk2)
1395  ABI_DEALLOCATE(shiftk3)
1396 
1397 end subroutine getkgrid

m_kpts/kpts_ibz_from_kptrlatt [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

  kpts_ibz_from_kptrlatt

FUNCTION

  Determines the irreducible wedge, the corresponding weights and the list
  of k-points in the Brillouin Zone starting from kptrlatt and the set shifts.

INPUTS

  cryst<crystal_t> = crystalline structure with info on symmetries and time-reversal.
  kptopt=option for the generation of k points
    (defines whether spatical symmetries and/or time-reversal can be used)
  kptrlatt(3,3)=integer coordinates of the primitive vectors of the
   lattice reciprocal to the k point lattice to be generated here
   If diagonal, the three values are the Monkhorst-Pack usual values, in case of simple cubic.
  nshiftk= number of shift vectors in the repeated cell
  shiftk(3,nshiftk) = vectors that will be used to determine the shifts from (0. 0. 0.).

OUTPUT

  nkibz,nkbz = Number of points in IBZ and BZ, respectively.
  The following arrays are allocated and returned by the routine:
  kibz(3,nkibz) = k-points in the IBZ.
  wtk(nkibz) = weights of the k-points in the IBZ (normalized to one).
  kbz(3,nkbz) = k-points in the BZ.
  [new_kptrlatt] = New value of kptrlatt returned by getkgrid
  [new_shiftk(3,new_nshiftk)] = New set of shifts returned by getkgrid

PARENTS

      m_dvdb,m_ebands,m_gruneisen,m_ifc,m_kpts,m_phgamma,m_phonons,m_sigmaph

CHILDREN

      getkgrid,init_tetra,kpts_ibz_from_kptrlatt,listkk

SOURCE

149 subroutine kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, nkibz, kibz, wtk, nkbz, kbz, &
150   new_kptrlatt, new_shiftk)  ! Optional
151 
152 
153 !This section has been created automatically by the script Abilint (TD).
154 !Do not modify the following lines by hand.
155 #undef ABI_FUNC
156 #define ABI_FUNC 'kpts_ibz_from_kptrlatt'
157 !End of the abilint section
158 
159  implicit none
160 
161 !Arguments ------------------------------------
162 !scalars
163  integer,intent(in) :: nshiftk,kptopt
164  integer,intent(out) :: nkibz,nkbz
165  type(crystal_t),intent(in) :: cryst
166 !arrays
167  integer,intent(in) :: kptrlatt(3,3)
168  integer,optional,intent(out) :: new_kptrlatt(3,3)
169  real(dp),intent(in) :: shiftk(3,nshiftk)
170  real(dp),allocatable,intent(out) :: wtk(:),kibz(:,:),kbz(:,:)
171  real(dp),optional,allocatable,intent(out) :: new_shiftk(:,:)
172 
173 !Local variables-------------------------------
174 !scalars
175  integer,parameter :: iout0=0,chksymbreak0=0,iscf2=2
176  integer :: my_nshiftk,nkpt_computed
177  real(dp) :: kptrlen
178 !arrays
179  integer,parameter :: vacuum0(3)=[0,0,0]
180  integer :: my_kptrlatt(3,3)
181  real(dp) :: my_shiftk(3,210)
182 
183 ! *********************************************************************
184 
185  ! First call to getkgrid to obtain the number of points in the BZ.
186  ABI_MALLOC(kibz, (3,0))
187  ABI_MALLOC(wtk, (0))
188 
189  ! Copy kptrlatt and shifts because getkgrid can change them
190  ! Be careful as getkgrid expects shiftk(3,210).
191  ABI_CHECK(nshiftk > 0 .and. nshiftk <= 210, "nshiftk must be in [1,210]")
192  my_nshiftk = nshiftk; my_shiftk = zero; my_shiftk(:,1:nshiftk) = shiftk
193  my_kptrlatt = kptrlatt
194 
195  call getkgrid(chksymbreak0,iout0,iscf2,kibz,kptopt,my_kptrlatt,kptrlen,&
196    cryst%nsym,0,nkibz,my_nshiftk,cryst%nsym,cryst%rprimd,my_shiftk,cryst%symafm,cryst%symrel,vacuum0,wtk)
197 
198  ABI_FREE(kibz)
199  ABI_FREE(wtk)
200 
201  ! Recall getkgrid to get kibz and wtk.
202  ABI_MALLOC(kibz, (3, nkibz))
203  ABI_MALLOC(wtk, (nkibz))
204 
205  call getkgrid(chksymbreak0,iout0,iscf2,kibz,kptopt,my_kptrlatt,kptrlen,&
206    cryst%nsym,nkibz,nkpt_computed,my_nshiftk,cryst%nsym,cryst%rprimd,my_shiftk,&
207    cryst%symafm,cryst%symrel,vacuum0,wtk,fullbz=kbz)
208 
209  nkbz = size(kbz, dim=2)
210 
211  ! Optionally, return new shifts and new_kptrlatt
212  if (present(new_shiftk)) then
213    ABI_MALLOC(new_shiftk, (3, my_nshiftk))
214    new_shiftk = my_shiftk(:, 1:my_nshiftk)
215  end if
216  if (present(new_kptrlatt)) new_kptrlatt = my_kptrlatt
217 
218  DBG_CHECK(abs(sum(wtk) - one) < tol10, "sum(wtk) != one")
219 
220 end subroutine kpts_ibz_from_kptrlatt

m_kpts/kpts_timrev_from_kptopt [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

  kpts_timrev_from_kptopt

FUNCTION

  Returns the value of timrev from kptopt
  1 if the use of time-reversal is allowed; 0 otherwise

INPUTS

  kptopt=option for the generation of k points
    (defines whether spatical symmetries and/or time-reversal can be used)

PARENTS

CHILDREN

SOURCE

 92 integer pure function kpts_timrev_from_kptopt(kptopt) result(timrev)
 93 
 94 
 95 !This section has been created automatically by the script Abilint (TD).
 96 !Do not modify the following lines by hand.
 97 #undef ABI_FUNC
 98 #define ABI_FUNC 'kpts_timrev_from_kptopt'
 99 !End of the abilint section
100 
101  implicit none
102 
103 !Arguments ------------------------------------
104 !scalars
105  integer,intent(in) :: kptopt
106 
107 ! *********************************************************************
108 
109  timrev = 1; if (any(kptopt == [3, 4])) timrev = 0
110 
111 end function kpts_timrev_from_kptopt

m_kpts/listkk [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 listkk

FUNCTION

 Given a list of nkpt1 initial k points kptns1 and a list of nkpt2
 final k points kptns2, associates each final k pt with a "closest"
 initial k point (or symmetric thereof, also taking possible umklapp)
 as determined by a metric gmet, that commutes with the symmetry operations.
 The algorithm does not scale as nkpt1 times nkpt2, thanks
 to the ordering of the kptns1 and kptns2 vectors according to their
 lengths, and comparison first between vectors of similar lengths.
 Returns indirect indexing list indkk.

INPUTS

  gmet(3,3)=reciprocal space metric (bohr^-2)
  kptns1(3,nkpt1)=list of initial k points (reduced coordinates)
  kptns2(3,nkpt2)=list of final k points
  nkpt1=number of initial k points
  nkpt2=number of final k points
  nsym=number of symmetry elements in space group
  sppoldbl=if 1, no spin-polarisation doubling
           if 2, spin-polarisation doubling using symafm
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  symmat(3,3,nsym)=symmetry operations (symrel or symrec, depending on
                   value of use_symrec
  timrev=1 if the use of time-reversal is allowed; 0 otherwise
  use_symrec: if present and true, symmat assumed to be symrec, otherwise assumed to be symrel (default)

OUTPUT

  dksqmax=maximal value of the norm**2 of the difference between
    a kpt2 vector and the closest k-point found from the kptns1 set, using symmetries.
  indkk(nkpt2*sppoldbl,6)=describe k point number of kpt1 that allows to
    generate wavefunctions closest to given kpt2
    if sppoldbl=2, use symafm to generate spin down wfs from spin up wfs

    indkk(:,1)=k point number of kptns1
    indkk(:,2)=symmetry operation to be applied to kpt1, to give kpt1a
      (if 0, means no symmetry operation, equivalent to identity )
    indkk(:,3:5)=shift in reciprocal space to be given to kpt1a,
      to give kpt1b, that is the closest to kpt2.
    indkk(:,6)=1 if time-reversal was used to generate kpt1a from kpt1, 0 otherwise

NOTES

  The tolerances tol12 and tol8 aims at giving a machine-independent ordering.
  (this trick is used in bonds.f, listkk.f, prtrhomxmn.f and rsiaf9.f)
  The tolerance tol12 is used for each component of the k vectors,
  and for the length of the vectors
  while the tolerance tol8 is used for the comparison of the squared lengths
  of the separate vectors.

PARENTS

      initberry,initorbmag,inwffil,m_dvdb,m_ebands,m_eprenorms,m_exc_diago
      m_fock,m_fstab,m_haydock,m_ifc,m_kpts,m_phgamma,m_sigmaph,mlwfovlp_qp

CHILDREN

      sort_dp,timab

SOURCE

587 subroutine listkk(dksqmax,gmet,indkk,kptns1,kptns2,nkpt1,nkpt2,nsym,&
588 & sppoldbl,symafm,symmat,timrev,use_symrec)
589 
590 
591 !This section has been created automatically by the script Abilint (TD).
592 !Do not modify the following lines by hand.
593 #undef ABI_FUNC
594 #define ABI_FUNC 'listkk'
595 !End of the abilint section
596 
597  implicit none
598 
599 !Arguments ------------------------------------
600 !scalars
601  integer,intent(in) :: nkpt1,nkpt2,nsym,sppoldbl,timrev
602  real(dp),intent(out) :: dksqmax
603  logical,optional,intent(in) :: use_symrec
604 !arrays
605  integer,intent(in) :: symafm(nsym),symmat(3,3,nsym)
606  integer,intent(out) :: indkk(nkpt2*sppoldbl,6)
607  real(dp),intent(in) :: gmet(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2)
608 
609 !Local variables-------------------------------
610 !scalars
611  integer :: l3,ig1,ig2,ig3,ii,ikpg1,ikpt1,ikpt2,ikpt2_done
612  integer :: ilarger,ismaller,itrial
613  integer :: isppol,isym,itimrev,jkpt1,jsym,jtime,limit
614  integer :: nsym_used,timrev_used,usesym
615  real(dp) :: dksq,dksqmn,lk2,llarger,ldiff,lsmaller,ltrial,min_l
616  character(len=500) :: message
617 !arrays
618  integer :: dkint(3),jdkint(3),k1int(3),k2int(3)
619  integer, allocatable :: isort(:)
620  real(dp) :: tsec(2)
621  real(dp) :: dk(3),kpg1(3),kpt1a(3),k1(3),k2(3)
622 !real(dp) :: kasq,ka(3)
623  real(dp),allocatable :: lkpg1(:),lkpg1_sorted(:)
624 
625 ! *************************************************************************
626 
627 !write(std_out,*)' listkk : nkpt1,nkpt2,nsym=',nkpt1,nkpt2,nsym
628  call timab(1021,1,tsec)
629 
630  if(sppoldbl<1 .or. sppoldbl>2)then
631    write(message, '(a,i4,3a)' )&
632 &   'The value of sppoldbl is',sppoldbl,',',ch10,&
633 &   'but it should be either 1 or 2.'
634    MSG_BUG(message)
635  end if
636 
637 !When usesym=0, the old way of converting the wavefunctions (without
638 !using the symmetries), is recovered.
639  usesym=1
640 
641  nsym_used=nsym
642  timrev_used=timrev
643  if(usesym==0)nsym_used=1
644  if(usesym==0)timrev_used=0
645 
646 !Precompute the length of the kpt1 vectors, also taking into account
647 !possible umpklapp vectors
648  limit=1 ; l3 = (2*limit+1)**3
649  ABI_ALLOCATE(lkpg1,(l3*nkpt1))
650  ABI_ALLOCATE(lkpg1_sorted,(l3*nkpt1))
651  ABI_ALLOCATE(isort,(l3*nkpt1))
652 !write(std_out,*)' List of kpt1 vectors '
653 !write(std_out,*)' Length of the kpt1 vectors :'
654 
655  do ikpt1=1,nkpt1
656    k1(:)=kptns1(:,ikpt1)
657 !  write(std_out,*)ikpt1,k1(:)
658    k1int(:)=nint(k1(:)+tol12)
659    k1(:)=k1(:)-k1int(:)
660    do ig1=-limit,limit
661      kpg1(1)=k1(1)+ig1
662      do ig2=-limit,limit
663        kpg1(2)=k1(2)+ig2
664        do ig3=-limit,limit
665          kpg1(3)=k1(3)+ig3
666 
667          ikpg1=ig1+limit+1 + (2*limit+1)*(ig2+limit) + (2*limit+1)**2*(ig3+limit) + l3*(ikpt1-1)
668 !        Compute the norm of the vector (also taking into account possible umklapp)
669          lkpg1(ikpg1)=sqrt(gmet(1,1)*kpg1(1)**2+gmet(2,2)*kpg1(2)**2+&
670 &         gmet(3,3)*kpg1(3)**2+two*(gmet(2,1)*kpg1(2)*kpg1(1)+&
671 &         gmet(3,2)*kpg1(3)*kpg1(2)+gmet(3,1)*kpg1(3)*kpg1(1)))
672          lkpg1_sorted(ikpg1)=lkpg1(ikpg1)
673          isort(ikpg1)=ikpg1
674 !        write(std_out,*)' ikpt1,ig1,ig2,ig3,lkpg1=',ikpt1,ig1,ig2,ig3,lkpg1(ikpg1)
675        end do
676      end do
677    end do
678  end do
679 
680  call sort_dp( l3*nkpt1,lkpg1_sorted,isort,tol12)
681 
682 !DEBUG
683 !write(std_out,*)' listkk : output list of kpt1 for checking purposes '
684 !write(std_out,*)' ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii)) '
685 !do ii=1,l3*nkpt1
686 !ikpt1=(isort(ii)-1)/l3+1
687 !write(std_out,*)ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii))
688 !enddo
689 !stop
690 !ENDDEBUG
691 
692  dksqmax=zero
693  do isppol=1,sppoldbl
694    do ikpt2=1,nkpt2
695 
696      ikpt2_done=0
697 !    Precompute the length of the kpt2 vector, with the Umklapp vector such that it is the closest to the Gamma point
698      k2(:)=kptns2(:,ikpt2)
699      k2int(:)=nint(k2(:)+tol12)
700      k2(:)=k2(:)-k2int(:)
701      lk2=sqrt(gmet(1,1)*k2(1)**2+gmet(2,2)*k2(2)**2+&
702 &     gmet(3,3)*k2(3)**2+two*(gmet(2,1)*k2(2)*k2(1)+&
703 &     gmet(3,2)*k2(3)*k2(2)+gmet(3,1)*k2(3)*k2(1)))
704 
705 !    DEBUG
706 !    write(std_out, '(a,i4,7es16.6)' )' listkk : ikpt2,kptns2(:,ikpt2),k2(:),lk2=',ikpt2,kptns2(:,ikpt2),k2(:),lk2
707 !    if(ikpt2/=17)cycle
708 !    ENDDEBUG
709 
710 !    Find the kpt1 vector whose length is the most similar to the length of lk2
711 !    up to a tolerance. Use a bissection algorithm.
712      ismaller=0       ; lsmaller=zero
713      ilarger=l3*nkpt1+1 ; llarger=huge(one)
714 !    This loop should never reach l3*nkpt1, since this is a bissection algorithm
715      do ii=1,l3*nkpt1
716        if((ilarger-ismaller)<2 .or. (llarger-lsmaller)<2*tol12)exit
717        itrial=(ilarger+ismaller)/2 ; ltrial=lkpg1_sorted(itrial)
718        if((ltrial-lk2)>tol12)then
719          ilarger=itrial ; llarger=ltrial
720        else if((ltrial-lk2)<-tol12)then
721          ismaller=itrial ; lsmaller=ltrial
722        else
723          ismaller=itrial ; lsmaller=ltrial
724          ilarger=itrial ; llarger=ltrial
725        end if
726      end do
727      itrial=ismaller
728      if(abs(llarger-lk2)<abs(lsmaller-lk2)-tol12)itrial=ilarger
729      if(itrial==0)itrial=ilarger
730      ismaller=itrial ; ilarger=itrial
731 !    write(std_out,*)' listkk : starting search at itrial=',itrial
732 
733      dksqmn=huge(one)
734 
735 !    The ii index is dummy. This avoids an infinite loop.
736      do ii=1,l3*nkpt1
737 !      do ikpt1=1,nkpt1
738 
739 !      If the difference in length between the trial vector and the target vector is bigger
740 !      than the already achieved distance, the search is finished ...
741        ldiff=abs(lkpg1_sorted(itrial)-lk2)
742 
743 
744 !      DEBUG
745 !      write(std_out,*)' listkk : ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,dksqmn=',ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,dksqmn
746 !      ENDDEBUG
747        if(ldiff**2>dksqmn+tol8)exit
748 
749 !      If this k-point has already been examined in a previous batch, skip it
750 !      First, compute the minimum of the difference of length of the sets of associated vectors thanks to Umklapp vectors
751 !      with the target vector
752        ikpt1=(isort(itrial)-1)/l3+1
753        min_l=minval(abs(lkpg1((ikpt1-1)*l3+1:(ikpt1-1)*l3+l3)-lk2))
754 !      Then compare with the current ldiff
755 
756 !      DEBUG
757 !      write(std_out,*)' listkk : ikpt1,min_l,ldiff=',ikpt1,min_l,ldiff
758 !      ENDDEBUG
759 
760        if(min_l > ldiff-tol12)then
761 
762 !        Now, will examine the trial vector, and the symmetric ones
763 !MG FIXME:
764 ! Here there's a possible problem with the order of symmetries because
765 ! in symkpt, time-reversal is the innermost loop. This can create inconsistencies in the symmetry tables.
766          do itimrev=0,timrev_used
767            do isym=1,nsym_used
768 
769 !            Select magnetic characteristic of symmetries
770              if(isppol==1 .and. symafm(isym)==-1)cycle
771              if(isppol==2 .and. symafm(isym)==1)cycle
772 
773 !            Compute symmetric point to kpt1
774              if(usesym==1)then
775 !              original code only used transpose(symrel)
776 !              kpt1a(:)=symrel(1,:,isym)*kptns1(1,ikpt1)+&
777 !              &             symrel(2,:,isym)*kptns1(2,ikpt1)+&
778 !              &             symrel(3,:,isym)*kptns1(3,ikpt1)
779                if (present(use_symrec)) then
780                  if (use_symrec) then
781                    kpt1a(:) = MATMUL(symmat(:,:,isym),kptns1(:,ikpt1))
782                  else
783                    kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1))
784                  end if
785                else
786                  kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1))
787                end if
788                kpt1a(:)=(1-2*itimrev)*kpt1a(:)
789              else
790                kpt1a(:)=kptns1(:,ikpt1)
791              end if
792 
793 !            Compute difference with respect to kpt2, modulo a lattice vector
794              dk(:)=kptns2(:,ikpt2)-kpt1a(:)
795              if(usesym==1)then
796 !              The tolerance insure similar behaviour on different platforms
797 !              XG120418 : Actually, *assumes* that the closest point will have reduced
798 !              coordinates differing by less than 1/2 . There might be elongated
799 !              cells where this is not correct ...
800                dkint(:)=nint(dk(:)+tol12)
801                dk(:)=dk(:)-dkint(:)
802              else
803                dkint(:)=0
804              end if
805 
806 !            Compute norm of the difference vector, and update kpt1 if better.
807              dksq=gmet(1,1)*dk(1)**2+gmet(2,2)*dk(2)**2+&
808 &             gmet(3,3)*dk(3)**2+two*(gmet(2,1)*dk(2)*dk(1)+&
809 &             gmet(3,2)*dk(3)*dk(2)+gmet(3,1)*dk(3)*dk(1))
810 
811              if (dksq<dksqmn+tol8) then
812 
813 !              If exactly the right point (without using symmetries neither umklapp vector), will exit the search
814 !              Note that in this condition, each coordinate is tested separately, without squaring. So, it is a much stronger
815 !              condition than dksqmn<tol12
816                if(sum(abs(kptns2(:,ikpt2)-kptns1(:,ikpt1)))<3*tol12)then
817                  ikpt2_done=1
818                end if
819 
820 !              Update in three cases : either if succeeded to have exactly the vector, or the distance is better,
821 !              or the distance is only slightly worsened so select the lowest itimrev, isym or ikpt1, in order to respect previous ordering
822                if(  ikpt2_done==1 .or. &
823 &               dksq+tol12<dksqmn .or. &
824 &               ( abs(dksq-dksqmn)<tol12 .and. &
825 &               ((itimrev<jtime) .or. &
826 &               (itimrev==jtime .and. isym<jsym) .or. &
827 &               (itimrev==jtime .and. isym==jsym .and. ikpt1<jkpt1))))then
828 
829                  dksqmn=dksq
830                  jkpt1=ikpt1
831                  jsym=isym
832                  jtime=itimrev
833                  jdkint(:)=dkint(:)
834 
835 !                DEBUG
836 !                write(std_out,*)' ikpt1,ikpt2=',ikpt1,ikpt2
837 !                write(std_out,*)' timrev_used=',timrev_used
838 !                write(std_out,*)' Succeeded to lower dskmn,ikpt2_done=',dksqmn,ikpt2_done
839 !                write(std_out,*)' ikpt1,isym,dkint(:),itimrev=',ikpt1,isym,dkint(:),itimrev
840 !                ka(:)=kpt1a(:)+dkint(:)
841 !                write(std_out,*)'        k1=',kpt1a(:)
842 !                write(std_out,*)'     dkint=',dkint(:)
843 !                write(std_out,*)' Actual k1=',ka(:)
844 !                write(std_out,*)'        k2=',kptns2(:,ikpt2)
845 !                kasq=gmet(1,1)*ka(1)**2+gmet(2,2)*ka(2)**2+&
846 !                &                  gmet(3,3)*ka(3)**2+two*(gmet(2,1)*ka(2)*ka(1)+&
847 !                &                  gmet(3,2)*ka(3)*ka(2)+gmet(3,1)*ka(3)*ka(1))
848 !                write(std_out,*)' Actual k1sq=',kasq
849 !                ENDDEBUG
850                end if
851 
852              end if
853              if(ikpt2_done==1)exit
854 
855            end do ! isym
856            if(ikpt2_done==1)exit
857 
858          end do ! itimrev
859          if(ikpt2_done==1)exit
860 
861        end if
862 
863 !      Update the interval that has been explored
864        if(itrial<ismaller)ismaller=itrial
865        if(itrial>ilarger)ilarger=itrial
866 
867 !      Select the next index to be tried (preferably the smaller indices, but this is a bit arbitrary).
868 
869 !      DEBUG
870 !      write(std_out,*)' before choosing the next index :'
871 !      write(std_out,*)' ismaller,itrial,ilarger=',ismaller,itrial,ilarger
872 !      write(std_out,*)' lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1)=',lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1)
873 !      ENDDEBUG
874        if(ismaller>1 .and. ilarger<l3*nkpt1)then
875          if(abs(lkpg1_sorted(ismaller-1)-lk2)<abs(lkpg1_sorted(ilarger+1)-lk2)+tol12)then
876            itrial=ismaller-1
877          else
878            itrial=ilarger+1
879          end if
880        end if
881        if(ismaller==1 .and. ilarger<l3*nkpt1)itrial=ilarger+1
882        if(ismaller>1 .and. ilarger==l3*nkpt1)itrial=ismaller-1
883 !      if(ismaller==1 .and. ilarger==l3*nkpt1), we are done with the loop !
884 
885      end do ! ikpt1
886 
887      indkk(ikpt2+(isppol-1)*nkpt2,1)=jkpt1
888      indkk(ikpt2+(isppol-1)*nkpt2,2)=jsym
889      indkk(ikpt2+(isppol-1)*nkpt2,3:5)=jdkint(:)
890      indkk(ikpt2+(isppol-1)*nkpt2,6)=jtime
891      dksqmax=max(dksqmax,dksqmn)
892 
893      if(dksqmn<-tol12)then
894        write(message, '(a,es16.6)' )'  The minimum square of dk has negative norm: dksqmn=',dksqmn
895        MSG_BUG(message)
896      end if
897 
898 !    DEBUG
899 !    write(std_out,'(a,i6,i2,2x,i6,5i3,es24.14)' )' listkk: ikpt2,isppol,indkk(ikpt2+(isppol-1)*nkpt2,:)=',ikpt2,isppol,indkk(ikpt2+(isppol-1)*nkpt2,:),dksqmn
900 !    if(nkpt1==17)stop
901 !    ENDDEBUG
902 
903    end do ! ikpt2
904  end do ! isppol
905 
906  ABI_DEALLOCATE(isort)
907  ABI_DEALLOCATE(lkpg1)
908  ABI_DEALLOCATE(lkpg1_sorted)
909 
910  call timab(1021,2,tsec)
911 
912 end subroutine listkk

m_kpts/mknormpath [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 mknormpath

FUNCTION

 Please do not use this  routine, use make_normpath instead.
 mknormpath should be removed

  This simple routine generates a normalized path that can be used to plot a band
  structures in an easy way. For normalized path we mean a path where the number
  of division on each segment is proportional to the length of the segment itself.
  To generate the above mentioned path, the subroutine must be called twice.
  The first call reports the total number of divisions in the normalized path, dimension
  that is required to correctly allocate the array.
  The second call calculates the reduced coordinates of the circuit.

COPYRIGHT

  Copyright (C) 2007-2018 ABINIT group (MG)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

INPUTS

 nbounds=number of points defining the path
 ndiv_small=number of points to be used to sample the smallest
  segment defined by bounds(:,1:nbounds)
 bounds(3,nbounds)=points defining the path
 gmet(3,3)=metric

OUTPUT

 ndiv(nbounds-1)= number of divisions for each segment
 npt_tot=total number of points sampled along the circuit
 path(3,npt_tot)= normalized path in reciprocal space

TODO

  Do not use this routine, it is obsolete and should be replaced by make_path in m_bz_mesh.

PARENTS

      inkpts

CHILDREN

      wrtout

SOURCE

3003 subroutine mknormpath(nbounds,bounds,gmet,ndiv_small,ndiv,npt_tot,path)
3004 
3005 
3006 !This section has been created automatically by the script Abilint (TD).
3007 !Do not modify the following lines by hand.
3008 #undef ABI_FUNC
3009 #define ABI_FUNC 'mknormpath'
3010 !End of the abilint section
3011 
3012  implicit none
3013 
3014 !Arguments ------------------------------------
3015  !F95 construct, interface required but we can call mknormpath once
3016  !real(dp),pointer :: path(:,:)
3017 !scalars
3018  integer,intent(in) :: nbounds,ndiv_small
3019  integer,intent(inout) :: npt_tot
3020 !arrays
3021  integer,intent(inout) :: ndiv(nbounds-1)
3022  real(dp),intent(in) :: bounds(3,nbounds),gmet(3,3)
3023  real(dp),intent(out),optional :: path(3,npt_tot)
3024 
3025 !Local variables-------------------------------
3026 !scalars
3027  integer :: idx,ii,jp
3028  real(dp) :: fct
3029  character(len=500) :: message
3030 !arrays
3031  real(dp) :: dd(3),lng(nbounds-1)
3032 
3033 ! *************************************************************************
3034 
3035  if (ndiv_small<=0) then
3036    write(message,'(3a,i0)')&
3037 &   'The argument ndiv_small should be a positive number,',ch10,&
3038 &   'however, ndiv_small=',ndiv_small
3039    MSG_ERROR(message)
3040  end if
3041 
3042  do ii=1,nbounds-1
3043    dd(:)=bounds(:,ii+1)-bounds(:,ii)
3044    lng(ii)= sqrt( dd(1)*gmet(1,1)*dd(1)+ &
3045 &   dd(2)*gmet(2,2)*dd(2)+ &
3046 &   dd(3)*gmet(3,3)*dd(3)+ &
3047 &   2.0d0*(dd(1)*gmet(1,2)*dd(2)+ &
3048 &   dd(1)*gmet(1,3)*dd(3)+ &
3049 &   dd(2)*gmet(2,3)*dd(3)) &
3050 &   )
3051  end do
3052  write(std_out,*)lng
3053  fct=minval(lng)
3054 
3055 !Avoid division by zero if k(:,i+1)=k(:,i)
3056  if (abs(fct)<tol6) then
3057    write(message,'(3a)')&
3058 &   'found two consecutive points in the path which are equal',ch10,&
3059 &   'This is not allowed, please modify the path in your input file'
3060    MSG_ERROR(message)
3061  end if
3062 
3063  fct=fct/ndiv_small
3064  ndiv(:)=nint(lng(:)/fct)
3065 !The 1 stand for the first point
3066  npt_tot=sum(ndiv)+1
3067 
3068 !allocate(path(3,npt_tot)
3069  if (.not.present(path)) then
3070    write(message,'(2a,i8)')ch10,&
3071 &   ' mknormpath : total number of points on the path : ',npt_tot
3072    call wrtout(std_out,message,'COLL')
3073    write(message,'(2a)')ch10,' Number of divisions for each segment of the normalized path : '
3074    call wrtout(std_out,message,'COLL')
3075    do ii=1,nbounds-1
3076      write(message,'(2(3f8.5,a),i5,a)')&
3077      bounds(:,ii),' ==> ',bounds(:,ii+1),' ( ndiv : ',ndiv(ii),' )'
3078      call wrtout(std_out,message,'COLL')
3079    end do
3080    write(message,'(a)')ch10
3081    call wrtout(std_out,message,'COLL')
3082  else
3083    write(message,'(2a)')ch10,' Normalized Path : '
3084    call wrtout(std_out,message,'COLL')
3085    idx=1
3086    do ii=1,nbounds-1
3087      do jp=1,ndiv(ii)
3088        path(:,idx)=bounds(:,ii)+(jp-1)*(path(:,ii+1)-path(:,ii))/ndiv(ii)
3089        write(message,'(i4,4x,3(f8.5,1x))')idx,path(:,idx)
3090        call wrtout(std_out,message,'COLL')
3091        idx=idx+1
3092      end do
3093    end do
3094  end if
3095 
3096 end subroutine mknormpath

m_kpts/smpbz [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 smpbz

FUNCTION

 Generate a set of special k (or q) points which samples in a homogeneous way
 the entire Brillouin zone of a simple lattice, face-centered cubic,
 body-centered lattice and hexagonal lattice.
 If kptrlatt is diagonal, the algorithm used here reduces to the usual
 Monkhorst-Pack set of k points.

INPUTS

  brav = 1 or -1 -> simple lattice; 2 -> face-centered cubic;
   3 -> body-centered lattice; 4 -> hexagonal lattice (D6h)
  downsampling(3) [optional, for brav=1 only]
    Three integer numbers, describing the downsampling of the k grid
    If present, in any case, only the first shiftk is taken into account
    The absolute value of one number gives, for the corresponding k-coordinate, the factor of decrease of the sampling
    If zero, only one point is used to sample along this direction
    The sign has also a meaning :
    - if three numbers are negative, perform a face-centered sampling
    - if two numbers are negative, perform a body-centered sampling
    - if one number is negative, perform a face-centered sampling for the two-dimensional lattice of the other directions
    - if one number is zero and at least one number is negative, perform face-centered sampling for the non-zero directions.
  iout = unit number for output
  kptrlatt(3,3)=integer coordinates of the primitive vectors of the
   lattice reciprocal to the k point lattice to be generated here
   If diagonal, the three values are the Monkhorst-Pack usual values, in case of simple cubic.
  mkpt = maximum number of k points
  nshiftk= number of shift vectors in the repeated cell
  option= Flag defining what will be printed of iout: 0 for k points, anything else for q points.
    Also, for q points, if the Gamma point is present, place it first in the list.
  shiftk(3,nshiftk) = vectors that will be used to determine the shifts from (0. 0. 0.).

OUTPUT

  nkpt = number of k points
  spkpt(3,mkpt) = the nkpt first values contain the special k points
   obtained by the Monkhorst & Pack method, in reduced coordinates.
   These vectors have to be multiplied by the reciprocal basis vectors
   gprimd(3,3) (in cartesian coordinates) to obtain the special k points
   set in cartesian coordinates.

NOTES

  also allows for more than one vector in repeated cell.
  this routine should be rewritten, to use the Wigner-Seitz cell,
  and thus unify the different treatments.
  References :
  H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13, 5188 (1976) [[cite:Monkhorst1976]]
  J.D. Pack and H.J. Monkhorst, Phys. Rev. B 16, 1748 (1977) [[cite:Pack1977]]
  A.H. MacDonald, Phys. Rev. B 18, 5897 (1978) [[cite:MacDonald1978]]
  R.A. Evarestov and V.P. Smirnov, Phys. Stat. Sol. (b) 119, 9 (1983) [[cite:Evarestov1983]]

PARENTS

      ep_setupqpt,getkgrid,harmonic_thermo,initberry,initorbmag,m_fstab,m_ifc
      m_tdep_abitypes

CHILDREN

      matr3inv,wrap2_pmhalf,wrtout

SOURCE

1715 subroutine smpbz(brav,iout,kptrlatt,mkpt,nkpt,nshiftk,option,shiftk,spkpt,downsampling)
1716 
1717 
1718 !This section has been created automatically by the script Abilint (TD).
1719 !Do not modify the following lines by hand.
1720 #undef ABI_FUNC
1721 #define ABI_FUNC 'smpbz'
1722 !End of the abilint section
1723 
1724  implicit none
1725 
1726 !Arguments -------------------------------
1727 !scalars
1728  integer,intent(in) :: brav,iout,mkpt,nshiftk,option
1729  integer,intent(out) :: nkpt
1730 !arrays
1731  integer,intent(in) :: kptrlatt(3,3)
1732  integer,optional,intent(in) :: downsampling(3)
1733  real(dp),intent(in) :: shiftk(3,nshiftk)
1734  real(dp),intent(out) :: spkpt(3,mkpt)
1735 
1736 !Local variables -------------------------
1737 !scalars
1738  integer,parameter :: prtvol=0
1739  integer :: dividedown,ii,ikshft,jj,kk,nkpout,nkptlatt,nn,proddown
1740  real(dp) :: shift
1741  character(len=500) :: message
1742 !arrays
1743  integer :: ads(3),boundmax(3),boundmin(3),cds(3),coord(3),ngkpt(3)
1744  integer, allocatable :: found1(:,:),found2(:,:),found3(:,:)
1745  real(dp) :: k1(3),k2(3),kcar(3),klatt(3,3),ktest(3),rlatt(3,3)
1746 
1747 ! *********************************************************************
1748 
1749 !DEBUG
1750 !write(std_out,*)' smpbz : brav,iout,mkpt,nkpt,option=',brav,iout,mkpt,nkpt,option
1751 !write(std_out,*)' smpbz : kptrlatt(:,:)=',kptrlatt(:,:)
1752 !write(std_out,*)' smpbz : nshiftk=',nshiftk
1753 !write(std_out,*)' smpbz : shiftk(:,:)=',shiftk(:,:)
1754 !write(std_out,*)' smpbz : downsampling(:)=',downsampling(:)
1755 !ENDDEBUG
1756 
1757  if(option/=0)then
1758    call wrtout(iout,'       Homogeneous q point set in the B.Z.  ','COLL')
1759  end if
1760 
1761  if(abs(brav)/=1)then
1762 !  Only generate Monkhorst-Pack lattices
1763    if(kptrlatt(1,2)/=0 .or. kptrlatt(2,1)/=0 .or. &
1764 &   kptrlatt(1,3)/=0 .or. kptrlatt(3,1)/=0 .or. &
1765 &   kptrlatt(2,3)/=0 .or. kptrlatt(3,2)/=0     ) then
1766      write(message, '(2a,a,3i4,a,a,3i4,a,a,3i4)' )&
1767 &     'When abs(brav)/=1, kptrlatt must be diagonal, while it is',ch10,&
1768 &     'kptrlatt(:,1)=',kptrlatt(:,1),ch10,&
1769 &     'kptrlatt(:,2)=',kptrlatt(:,2),ch10,&
1770 &     'kptrlatt(:,3)=',kptrlatt(:,3)
1771      MSG_BUG(message)
1772    end if
1773 
1774    ngkpt(1)=kptrlatt(1,1)
1775    ngkpt(2)=kptrlatt(2,2)
1776    ngkpt(3)=kptrlatt(3,3)
1777 !
1778    if( (ngkpt(1)<=0.or.ngkpt(2)<=0.or.ngkpt(3)<=0) .and. (ngkpt(1)/=0.or.ngkpt(2)/=0.or.ngkpt(3)/=0) ) then
1779      write(message, '(5a,i4,a,a,i4,a,a,i4,a,a)' )&
1780 &     'All ngkpt (or ngqpt) must be strictly positive',ch10,&
1781 &     'or all ngk(q)pt must be zero (for Gamma sampling), but :',ch10,&
1782 &     'ngk(q)pt(1) = ',ngkpt(1),ch10,&
1783 &     'ngk(q)pt(2) = ',ngkpt(2),ch10,&
1784 &     'ngk(q)pt(3) = ',ngkpt(3),ch10,&
1785 &     'Action: correct ngkpt or ngqpt in the input file.'
1786      MSG_BUG(message)
1787    end if
1788  end if
1789 
1790 !Just in case the user wants the grid downsampled to the Gamma point, checks that it is present, and possibly exits
1791  if(present(downsampling))then
1792    if(sum(abs(downsampling(:)))==0)then
1793      do ikshft=1,nshiftk
1794        if(sum(abs(shiftk(:,ikshft)))>tol12)cycle
1795        nkpt=1
1796        spkpt(:,1)=zero
1797        return
1798      end do
1799    end if
1800  end if
1801 
1802 !*********************************************************************
1803 
1804  if(abs(brav)==1)then
1805 
1806 !  Compute the number of k points in the G-space unit cell
1807 !  (will be multiplied by nshiftk later).
1808    nkptlatt=kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3) &
1809 &   +kptrlatt(1,2)*kptrlatt(2,3)*kptrlatt(3,1) &
1810 &   +kptrlatt(1,3)*kptrlatt(2,1)*kptrlatt(3,2) &
1811 &   -kptrlatt(1,2)*kptrlatt(2,1)*kptrlatt(3,3) &
1812 &   -kptrlatt(1,3)*kptrlatt(2,2)*kptrlatt(3,1) &
1813 &   -kptrlatt(1,1)*kptrlatt(2,3)*kptrlatt(3,2)
1814 
1815    if(present(downsampling))then
1816      if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then
1817        if(nshiftk>1)then
1818          write(message, '(a,3i4,2a,i4,4a)' )&
1819 &         'Real downsampling is activated, with downsampling(1:3)=',downsampling(1:3),ch10,&
1820 &         'However, nshiftk must be 1 in this case, while the input nshiftk=',nshiftk,ch10,&
1821 &         'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,&
1822 &         'or set nshiftk=1.'
1823          MSG_ERROR(message)
1824        end if
1825        proddown=downsampling(1)*downsampling(2)*downsampling(3)
1826        if(proddown/=0)then
1827          dividedown=abs(proddown)
1828          if(minval(downsampling(:))<0)then   ! If there is at least one negative number
1829            dividedown=dividedown*2
1830            if(proddown>0)dividedown=dividedown*2 ! If there are two negative numbers
1831          end if
1832        end if
1833        if(mod(nkptlatt,dividedown)==0)then
1834          nkptlatt=nkptlatt/dividedown
1835        else
1836          write(message, '(a,3i4,2a,i4,4a)' )&
1837 &         'The requested downsampling, with downsampling(1:3)=',downsampling(1:3),ch10,&
1838 &         'is not compatible with kptrlatt=',ch10,&
1839 &         kptrlatt(:,:),ch10,&
1840 &         'that gives nkptlatt=',nkptlatt,ch10,&
1841 &         'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,&
1842 &         'or modify your k-point grid and/or your downsampling in order for them to be compatible.'
1843          MSG_ERROR(message)
1844        end if
1845      end if
1846    end if
1847 
1848 !  Simple Lattice
1849    if (prtvol > 0) call wrtout(std_out,'       Simple Lattice Grid ','COLL')
1850    if (mkpt<nkptlatt*nshiftk) then
1851      write(message, '(a,a,a,i8,a,a,a,a,a)' )&
1852 &     'The value of mkpt is not large enough. It should be',ch10,&
1853 &     'at least',nkptlatt*nshiftk,',',ch10,&
1854 &     'Action: set mkpt to that value in the main routine,',ch10,&
1855 &     'and recompile the code.'
1856      MSG_BUG(message)
1857    end if
1858 
1859 !  Build primitive vectors of the k lattice
1860    rlatt(:,:)=kptrlatt(:,:)
1861    call matr3inv(rlatt,klatt)
1862 
1863 !DEBUG
1864 !        write(std_out,*)' First primitive vector of the k lattice :',klatt(:,1)
1865 !        write(std_out,*)' Second primitive vector of the k lattice :',klatt(:,2)
1866 !        write(std_out,*)' Third primitive vector of the k lattice :',klatt(:,3)
1867 !ENDDEBUG
1868 
1869 !  Now, klatt contains the three primitive vectors of the k lattice,
1870 !  in reduced coordinates. One builds all k vectors that
1871 !  are contained in the first Brillouin zone, with coordinates
1872 !  in the interval [0,1[ . First generate boundaries of a big box.
1873 
1874    do jj=1,3
1875 
1876 !    Mathematically, one has to find the coordinates of the corners of a
1877 !    rectangular paralleliped with integer coordinates, that multiplies the klatt primitive cell and allows
1878 !    it to incorporate completely the [0,1]^3 box. Then take the minimum and maximum
1879 !    of these coordinates, and round them negatively and positively to the next integer.
1880 !    This can be done easily using kptrlatt, considering each coordinate in turn
1881 !    and boils down to enlarging the boundaries for jj by the value of kptrlatt(:,jj),
1882 !    acting on boundmin or boundmax depending on the sign ot kptrlatt(:,jj).
1883 !    XG171020 The coding before 171020 was correct, despite being very simple.
1884      boundmin(jj)=0 ; boundmax(jj)=0
1885      do ii=1,3
1886        if(kptrlatt(ii,jj)<0)boundmin(jj)=boundmin(jj)+kptrlatt(ii,jj)
1887        if(kptrlatt(ii,jj)>0)boundmax(jj)=boundmax(jj)+kptrlatt(ii,jj)
1888      end do
1889 
1890 !    To accomodate the shifts, boundmin and boundmax don't start from 0, but are enlarged by one
1891 !    positively and/or negatively.
1892 !    XG171020 Coding in v8.6.0 and before was not correct. This one is even simpler actually.
1893      boundmin(jj)=boundmin(jj)-ceiling(maxval(shiftk(jj,:))+tol14)
1894      boundmax(jj)=boundmax(jj)-floor(minval(shiftk(jj,:))-tol14)
1895 
1896    end do
1897 
1898    if(present(downsampling))then
1899      ABI_ALLOCATE(found1,(boundmin(2):boundmax(2),boundmin(3):boundmax(3)))
1900      ABI_ALLOCATE(found2,(boundmin(1):boundmax(1),boundmin(3):boundmax(3)))
1901      ABI_ALLOCATE(found3,(boundmin(1):boundmax(1),boundmin(2):boundmax(2)))
1902      found1=0 ; found2=0 ; found3=0
1903    end if
1904 
1905    nn=1
1906    do kk=boundmin(3),boundmax(3)
1907      coord(3)=kk
1908      do jj=boundmin(2),boundmax(2)
1909        coord(2)=jj
1910        do ii=boundmin(1),boundmax(1)
1911          coord(1)=ii
1912 
1913 !        Here, apply the downsampling : skip some of the trials
1914          if(present(downsampling))then
1915 
1916            if(downsampling(1)==0 .and. found1(coord(2),coord(3))==1)cycle
1917            if(downsampling(2)==0 .and. found2(coord(1),coord(3))==1)cycle
1918            if(downsampling(3)==0 .and. found3(coord(1),coord(2))==1)cycle
1919 
1920            ads(:)=abs(downsampling(:))
1921            if(ads(1)>0 .and. mod(coord(1),ads(1))/=0)cycle
1922            if(ads(2)>0 .and. mod(coord(2),ads(2))/=0)cycle
1923            if(ads(3)>0 .and. mod(coord(3),ads(2))/=0)cycle
1924            cds(:)=coord(:)/ads(:)
1925            if(minval(downsampling(:))<0)then   ! If there is at least one negative number
1926 
1927              if(downsampling(1)*downsampling(2)*downsampling(3)/=0)then  ! If there is no zero number
1928 !              Face-centered case
1929                if(downsampling(1)<0 .and. downsampling(2)<0 .and. downsampling(3)<0)then ! All three are negative
1930                  if(mod(sum(cds(:)),2)/=0)cycle
1931 !              One-face-centered case
1932                else if(downsampling(1)*downsampling(2)*downsampling(3)<0)then  ! Only one is negative
1933                  if(downsampling(1)<0 .and. mod(cds(2)+cds(3),2)/=0)cycle
1934                  if(downsampling(2)<0 .and. mod(cds(1)+cds(3),2)/=0)cycle
1935                  if(downsampling(3)<0 .and. mod(cds(1)+cds(2),2)/=0)cycle
1936 !              Body-centered case ! What is left : two are negative
1937                else
1938                  if(sum(mod(cds(:),2))==1 .or. sum(mod(cds(:),2))==2)cycle ! Either all are zero, or all are one, so skip when sum is 1 or 2.
1939                end if
1940              else
1941                if(downsampling(1)==0 .and. mod(cds(2)+cds(3),2)/=0)cycle
1942                if(downsampling(2)==0 .and. mod(cds(1)+cds(3),2)/=0)cycle
1943                if(downsampling(3)==0 .and. mod(cds(1)+cds(2),2)/=0)cycle
1944              end if
1945            end if
1946          end if
1947 
1948          do ikshft=1,nshiftk
1949 
1950 !          Only the first shiftk is taken into account if downsampling
1951 !          if(.false.)then
1952            if(present(downsampling))then
1953              if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then
1954                if(ikshft>1)cycle
1955              end if
1956            end if
1957 
1958 !          Coordinates of the trial k point with respect to the k primitive lattice
1959            k1(1)=ii+shiftk(1,ikshft)
1960            k1(2)=jj+shiftk(2,ikshft)
1961            k1(3)=kk+shiftk(3,ikshft)
1962 !          Reduced coordinates of the trial k point
1963            k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3)
1964 !          Eliminate the point if outside [0,1[
1965            if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle
1966            if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle
1967            if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle
1968 !          Wrap the trial values in the interval ]-1/2,1/2] .
1969            call wrap2_pmhalf(k2(1),k1(1),shift)
1970            call wrap2_pmhalf(k2(2),k1(2),shift)
1971            call wrap2_pmhalf(k2(3),k1(3),shift)
1972            spkpt(:,nn)=k1(:)
1973            nn=nn+1
1974 
1975            if(present(downsampling))then
1976              found1(coord(2),coord(3))=1
1977              found2(coord(1),coord(3))=1
1978              found3(coord(1),coord(2))=1
1979            end if
1980 
1981          end do
1982        end do
1983      end do
1984    end do
1985    nkpt=nn-1
1986 
1987    if(present(downsampling))then
1988      ABI_DEALLOCATE(found1)
1989      ABI_DEALLOCATE(found2)
1990      ABI_DEALLOCATE(found3)
1991    end if
1992 
1993    if(nkpt/=nkptlatt*nshiftk)then
1994      write(message, '(a,i8,a,a,a,i8,a)' )&
1995 &     'The number of k points ',nkpt,'  is not equal to',ch10,&
1996 &     'nkptlatt*nshiftk which is',nkptlatt*nshiftk,'.'
1997 !DEBUG
1998 ! write(std_out,*)' smpbz : brav,iout,mkpt,nkpt,option=',brav,iout,mkpt,nkpt,option
1999 ! write(std_out,*)' smpbz : kptrlatt(:,:)=',kptrlatt(:,:)
2000 ! write(std_out,*)' smpbz : nshiftk=',nshiftk
2001 ! write(std_out,*)' smpbz : shiftk(:,:)=',shiftk(:,:)
2002 ! write(std_out,*)' smpbz : downsampling(:)=',downsampling(:)
2003 ! write(std_out,*)message
2004 ! stop
2005 !ENDDEBUG
2006      MSG_BUG(message)
2007    end if
2008 
2009  else if(brav==2)then
2010 
2011 !  Face-Centered Lattice
2012    if (prtvol > 0) call wrtout(std_out,'       Face-Centered Lattice Grid ','COLL')
2013    if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2) then
2014      write(message, '(a,a,a,i8,a,a,a,a,a)' )&
2015 &     'The value of mkpt is not large enough. It should be',ch10,&
2016 &     'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,',',ch10,&
2017 &     'Action: set mkpt to that value in the main routine,',ch10,&
2018 &     'and recompile the code.'
2019      MSG_BUG(message)
2020    end if
2021    nn=1
2022    if (ngkpt(1)/=ngkpt(2).or.ngkpt(1)/=ngkpt(3)) then
2023      write(message, '(4a,3(a,i6,a),a)' )&
2024 &     'For face-centered lattices, the numbers ngqpt(1:3)',ch10,&
2025 &     'must be equal, while they are :',ch10,&
2026 &     'ngqpt(1) = ',ngkpt(1),ch10,&
2027 &     'ngqpt(2) = ',ngkpt(2),ch10,&
2028 &     'ngqpt(3) = ',ngkpt(3),ch10,&
2029 &     'Action: modify ngqpt(1:3) in the input file.'
2030      MSG_BUG(message)
2031    end if
2032    if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2)) then
2033      write(message, '(4a,3(a,i6,a),a)' )&
2034 &     'For face-centered lattices, the numbers ngqpt(1:3)*nshiftk',ch10,&
2035 &     'must be even, while they are :',ch10,&
2036 &     'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,&
2037 &     'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,&
2038 &     'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,&
2039 &     'Action: modify ngqpt(1:3)*nshiftk in the input file.'
2040      MSG_ERROR(message)
2041    end if
2042    if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then
2043      spkpt(1,1)=0.0_dp
2044      spkpt(2,1)=0.0_dp
2045      spkpt(3,1)=0.0_dp
2046      nkpt=1
2047    else
2048      do kk=1,ngkpt(3)
2049        do jj=1,ngkpt(2)
2050          do ii=1,ngkpt(1)
2051            do ikshft=1,nshiftk
2052              k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1)
2053              k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2)
2054              k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3)
2055 !            Wrap the trial values in the interval ]-1/2,1/2] .
2056              call wrap2_pmhalf(k1(1),k2(1),shift)
2057              call wrap2_pmhalf(k1(2),k2(2),shift)
2058              call wrap2_pmhalf(k1(3),k2(3),shift)
2059 !            Test whether it is inside the FCC BZ.
2060              ktest(1)=2*k2(1)-1.0d-10
2061              ktest(2)=2*k2(2)-2.0d-10
2062              ktest(3)=2*k2(3)-5.0d-10
2063              if (abs(ktest(1))+abs(ktest(2))+abs(ktest(3))<1.5_dp) then
2064                kcar(1)=ktest(1)+1.0d-10
2065                kcar(2)=ktest(2)+2.0d-10
2066                kcar(3)=ktest(3)+5.0d-10
2067                spkpt(1,nn)=0.5_dp*kcar(2)+0.5_dp*kcar(3)
2068                spkpt(2,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(3)
2069                spkpt(3,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(2)
2070                nn=nn+1
2071              end if
2072            end do
2073          end do
2074        end do
2075      end do
2076      nkpt=nn-1
2077      if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2)then
2078        write(message, '(a,i8,a,a,a,i8,a)' )&
2079 &       'The number of k points ',nkpt,'  is not equal to',ch10,&
2080 &       '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2 which is',&
2081 &       (ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,'.'
2082        MSG_BUG(message)
2083      end if
2084    end if
2085 
2086  else if(brav==3)then
2087 
2088 !  Body-Centered Lattice (not mandatory cubic !)
2089    if (prtvol > 0) call wrtout(std_out,'       Body-Centered Lattice Grid ','COLL')
2090    if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/4) then
2091      write(message, '(a,a,a,i8,a,a,a,a,a)' )&
2092 &     'The value of mkpt is not large enough. It should be',ch10,&
2093 &     'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,',',ch10,&
2094 &     'Action: set mkpt to that value in the main routine,',ch10,&
2095 &     'and recompile the code.'
2096      MSG_BUG(message)
2097    end if
2098    nn=1
2099    if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2) .or.&
2100 &   (ngkpt(2)*nshiftk)/=(((ngkpt(2)*nshiftk)/2)*2) .or.&
2101 &   (ngkpt(3)*nshiftk)/=(((ngkpt(3)*nshiftk)/2)*2) ) then
2102      write(message, '(4a,3(a,i6,a),a)' )&
2103 &     'For body-centered lattices, the numbers ngqpt(1:3)',ch10,&
2104 &     'must be even, while they are :',ch10,&
2105 &     'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,&
2106 &     'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,&
2107 &     'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,&
2108 &     'Action: modify ngqpt(1:3) in the input file.'
2109      MSG_ERROR(message)
2110    end if
2111    if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then
2112      spkpt(1,1)=0.0_dp
2113      spkpt(2,1)=0.0_dp
2114      spkpt(3,1)=0.0_dp
2115      nkpt=1
2116    else
2117      do kk=1,ngkpt(3)
2118        do jj=1,ngkpt(2)
2119          do ii=1,ngkpt(1)
2120            do ikshft=1,nshiftk
2121              k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1)
2122              k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2)
2123              k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3)
2124 !            Wrap the trial values in the interval ]-1/2,1/2] .
2125              call wrap2_pmhalf(k1(1),k2(1),shift)
2126              call wrap2_pmhalf(k1(2),k2(2),shift)
2127              call wrap2_pmhalf(k1(3),k2(3),shift)
2128 !            Test whether it is inside the BCC BZ.
2129              ktest(1)=2*k2(1)-1.0d-10
2130              ktest(2)=2*k2(2)-2.0d-10
2131              ktest(3)=2*k2(3)-5.0d-10
2132              if (abs(ktest(1))+abs(ktest(2))<1._dp) then
2133                if (abs(ktest(1))+abs(ktest(3))<1._dp) then
2134                  if (abs(ktest(2))+abs(ktest(3))<1._dp) then
2135                    kcar(1)=ktest(1)+1.0d-10
2136                    kcar(2)=ktest(2)+2.0d-10
2137                    kcar(3)=ktest(3)+5.0d-10
2138                    spkpt(1,nn)=-0.5*kcar(1)+0.5*kcar(2)+0.5*kcar(3)
2139                    spkpt(2,nn)=0.5*kcar(1)-0.5*kcar(2)+0.5*kcar(3)
2140                    spkpt(3,nn)=0.5*kcar(1)+0.5*kcar(2)-0.5*kcar(3)
2141                    nn=nn+1
2142                  end if
2143                end if
2144              end if
2145            end do
2146          end do
2147        end do
2148      end do
2149      nkpt=nn-1
2150      if(nkpt==0)then
2151        write(message, '(3a)' )&
2152 &       'BCC lattice, input ngqpt=0, so no kpt is generated.',ch10,&
2153 &       'Action: modify ngqpt(1:3) in the input file.'
2154        MSG_ERROR(message)
2155      end if
2156      if(nkpt/=(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4)then
2157        write(message, '(a,i8,a,a,a,i8,a)' )&
2158 &       'The number of k points ',nkpt,'  is not equal to',ch10,&
2159 &       '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4 which is',&
2160 &       (ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,'.'
2161        MSG_BUG(message)
2162      end if
2163    end if
2164 
2165  else if(brav==4)then
2166 
2167 !  Hexagonal Lattice  (D6h)
2168    if (prtvol > 0) call wrtout(std_out,'       Hexagonal Lattice Grid ','COLL')
2169    if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)) then
2170      write(message, '(a,a,a,i8,a,a,a,a,a)' )&
2171 &     'The value of mkpt is not large enough. It should be',ch10,&
2172 &     'at least',ngkpt(1)*ngkpt(2)*ngkpt(3),',',ch10,&
2173 &     'Action: set mkpt to that value in the main routine,',ch10,&
2174 &     'and recompile the code.'
2175      MSG_BUG(message)
2176    end if
2177    nn=1
2178    if (ngkpt(1)/=ngkpt(2)) then
2179      write(message, '(4a,2(a,i6,a),a)' )&
2180 &     'For hexagonal lattices, the numbers ngqpt(1:2)',ch10,&
2181 &     'must be equal, while they are :',ch10,&
2182 &     'ngqpt(1) = ',ngkpt(1),ch10,&
2183 &     'ngqpt(2) = ',ngkpt(2),ch10,&
2184 &     'Action: modify ngqpt(1:3) in the input file.'
2185      MSG_ERROR(message)
2186    end if
2187    if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then
2188      write(message, '(3a)' )&
2189 &     'For hexagonal lattices, ngqpt(1:3)=0 is not permitted',ch10,&
2190 &     'Action: modify ngqpt(1:3) in the input file.'
2191      MSG_ERROR(message)
2192    else
2193      do kk=1,ngkpt(3)
2194        do jj=1,ngkpt(2)
2195          do ii=1,ngkpt(1)
2196            do ikshft=1,nshiftk
2197              k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1)
2198              k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2)
2199              k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3)
2200 !            Wrap the trial values in the interval ]-1/2,1/2] .
2201              call wrap2_pmhalf(k1(1),k2(1),shift)
2202              call wrap2_pmhalf(k1(2),k2(2),shift)
2203              call wrap2_pmhalf(k1(3),k2(3),shift)
2204              spkpt(:,nn)=k2(:)
2205              nn=nn+1
2206            end do
2207          end do
2208        end do
2209      end do
2210      nkpt=nn-1
2211      if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)then
2212        write(message, '(a,i8,a,a,a,i8,a)' )&
2213 &       'The number of k points ',nkpt,'  is not equal to',ch10,&
2214 &       'ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk which is',&
2215 &       ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk,'.'
2216        MSG_BUG(message)
2217      end if
2218    end if
2219 
2220  else
2221 
2222    write(message, '(a,i6,a,a,a)' )&
2223 &   'The calling routine asks brav=',brav,'.',ch10,&
2224 &   'but only brav=1 or -1,2,3 or 4 are allowed.'
2225    MSG_BUG(message)
2226  end if
2227 
2228  if (option/=0) then
2229 !  Put the Gamma point first
2230    if(nkpt>1)then
2231      do ii=1,nkpt
2232        if(sum(abs(spkpt(:,ii)))<tol8)then
2233          spkpt(:,ii)=spkpt(:,1)
2234          spkpt(:,1)=zero
2235          exit
2236        end if
2237      end do
2238    end if
2239 
2240    write(message,'(a,i8)')' Grid q points  : ',nkpt
2241    call wrtout(iout,message,'COLL')
2242    nkpout=nkpt
2243    if(nkpt>80)then
2244      call wrtout(iout,' greater than 80, so only write 20 of them ','COLL')
2245      nkpout=20
2246    end if
2247    do ii=1,nkpout
2248      write(message, '(1x,i2,a2,3es16.8)' )ii,') ',spkpt(1,ii),spkpt(2,ii),spkpt(3,ii)
2249      call wrtout(iout,message,'COLL')
2250    end do
2251  end if
2252 
2253 end subroutine smpbz

m_kpts/symkchk [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 symkchk

FUNCTION

 Checks that the set of k points chosen for a response function
 calculation has the full space group symmetry, modulo time reversal if appropriate.
 Returns ierr/=0 with error message if not satisfied
 Currently used only when strain perturbation is treated. Based on symkpt.

INPUTS

 kptns(3,nkpt)= k vectors in reciprocal space
 nkpt = number of k-points whose weights are wtk
 nsym=number of space group symmetries
 symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space)
 timrev: if 1, the time reversal operation has to be taken into account
 if 0, no time reversal symmetry.

OUTPUT

  msg=Error message if ierr /= 0

TODO

  This version should scale badly with the number of k-points. Replace loops with listkk

PARENTS

      respfn

CHILDREN

      wrtout

SOURCE

402 integer function symkchk(kptns,nkpt,nsym,symrec,timrev,errmsg) result(ierr)
403 
404 
405 !This section has been created automatically by the script Abilint (TD).
406 !Do not modify the following lines by hand.
407 #undef ABI_FUNC
408 #define ABI_FUNC 'symkchk'
409 !End of the abilint section
410 
411  implicit none
412 
413 !Arguments -------------------------------
414 !scalars
415  integer,intent(in) :: nkpt,nsym,timrev
416  character(len=*),intent(out) :: errmsg
417 !arrays
418  integer,intent(in) :: symrec(3,3,nsym)
419  real(dp),intent(in) :: kptns(3,nkpt)
420 
421 !Local variables -------------------------
422 !scalars
423  integer :: identi,ii,ikpt,ikpt2,imatch,isym,jj,tident
424  real(dp) :: difk,reduce
425  character(len=500) :: message
426 !arrays
427  real(dp) :: ksym(3)
428 
429 ! *********************************************************************
430  ierr = 0
431 
432  if(timrev/=1 .and. timrev/=0)then
433    write(errmsg, '(3a,i0,a)' )&
434 &   'timrev should be 0 or 1, while',ch10,&
435 &   'it is equal to ',timrev,'.'
436    ierr = 1; return
437  end if
438 
439  if(nsym/=1)then
440 !  Find the identity symmetry operation
441    do isym=1,nsym
442      tident=1
443      do jj=1,3
444        if(symrec(jj,jj,isym)/=1)tident=0
445        do ii=1,3
446          if( ii/=jj .and.&
447 &         symrec(ii,jj,isym)/=0)tident=0
448        end do
449      end do
450      if(tident==1)then
451        identi=isym
452        call wrtout(std_out,sjoin(' symkchk: found identity with number:', itoa(identi)))
453        exit
454      end if
455    end do
456    if(tident==0)then
457      errmsg = 'Did not found the identity operation.'
458      ierr = 1; return
459    end if
460  end if
461 
462 !Here begins the serious business
463 !The length sorting, etc. of symkpt have been dropped because the
464 !computational cost is estimated to be negligible.
465 
466  if(nsym>1 .or. timrev==1)then
467 
468 !  Outer loop over kpts
469    do ikpt=1,nkpt-1
470 
471 !    Loop on the symmetries
472 !    For each k-point and each symmetry transformation, a matching
473 !    k-pointpt must be found, modulo time reversal if appropriate
474      do isym=1,nsym
475 
476 !      Get the symmetric of the vector
477        do ii=1,3
478          ksym(ii)= kptns(1,ikpt)*symrec(ii,1,isym)&
479 &         +kptns(2,ikpt)*symrec(ii,2,isym)&
480 &         +kptns(3,ikpt)*symrec(ii,3,isym)
481        end do
482 
483 !      Second loop k-points
484        do ikpt2=1,nkpt
485 
486 !        Test for match of symmetric and any vector (including original)
487          imatch=1
488          do ii=1,3
489            difk= ksym(ii)-kptns(ii,ikpt2)
490            reduce=difk-anint(difk)
491            if(abs(reduce)>tol8)imatch=0
492          end do
493          if(imatch==1)exit
494 
495 !        Test for match with time reversal
496          if(timrev==1)then
497            imatch=1
498            do ii=1,3
499              difk= ksym(ii)+kptns(ii,ikpt2)
500              reduce=difk-anint(difk)
501              if(abs(reduce)>tol8)imatch=0
502            end do
503            if(imatch==1)exit
504          end if
505 
506        end do ! End secondary loop over k-points
507        if (imatch/=1) then
508          write(errmsg, '(a,a,a,i4,a,i4,a,a,a,a)' )&
509 &         'k-point set must have full space-group symmetry',ch10,&
510 &         'there is no match for kpt',ikpt,' transformed by symmetry',isym,ch10,&
511 &         'Action: change kptopt to 2 or 3 and/or change or use shiftk',ch10,&
512 &         'shiftk = 0 0 0 is always a safe choice.'
513          ierr = 2; return
514        end if
515 
516      end do ! End loop on isym
517    end do ! End primary loop over k-points
518 
519    write(message,'(a)')' symkchk : k-point set has full space-group symmetry.'
520    call wrtout(std_out,message,'COLL')
521    call wrtout(ab_out,message,'COLL')
522  end if
523 
524 end function symkchk

m_kpts/testkgrid [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 testkgrid

FUNCTION

 Test different grids of k points. The algorithm used is based on the idea of testing different
 one-dimensional sets of possible k point grids. It is not exhaustive (other families could be included),
 but should do a respectable job in all cases. The Monkhorst-Pack set of grids (defined with respect to
 symmetry axes, and not primitive axes) is always tested.

INPUTS

  bravais(11): bravais(1)=iholohedry
               bravais(2)=center
               bravais(3:11)=coordinates of rprim in the axes of the conventional bravais lattice (*2 if center/=0)
  iout=unit number for echoed output
  msym=default maximal number of symmetries
  nsym=number of symetries
  prtkpt=if non-zero, will write the characteristics of k grids, then stop
  rprimd(3,3)=dimensional real space primitive translations (bohr)
  symafm(nsym)=(anti)ferromagnetic part of symmetry operations
  symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
  vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum

OUTPUT

  kptrlatt(3,3)=k-point lattice specification
  nshiftk=number of k-point shifts in shiftk (always 1 from this routine)
  shiftk(3,210)=shift vectors for k point generation

SIDE EFFECTS

  kptrlen=length of the smallest real space supercell vector associated with the lattice of k points.

NOTES

 Note that nkpt can be computed by calling this routine with input value nkpt=0
 Note that kptopt is always =1 in this routine.

PARENTS

      inkpts,m_ab7_kpoints

CHILDREN

      getkgrid,abi_abort,matr3inv,metric,smallprim,wrtout,xmpi_abort

SOURCE

2299 subroutine testkgrid(bravais,iout,kptrlatt,kptrlen,&
2300 & msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum)
2301 
2302 
2303 !This section has been created automatically by the script Abilint (TD).
2304 !Do not modify the following lines by hand.
2305 #undef ABI_FUNC
2306 #define ABI_FUNC 'testkgrid'
2307 !End of the abilint section
2308 
2309  implicit none
2310 
2311 !Arguments ------------------------------------
2312 !scalars
2313  integer,intent(in) :: iout,msym,nsym,prtkpt
2314  integer,intent(out) :: nshiftk
2315  real(dp),intent(inout) :: kptrlen
2316 !arrays
2317  integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3)
2318  integer,intent(out) :: kptrlatt(3,3)
2319  real(dp),intent(in) :: rprimd(3,3)
2320  real(dp),intent(inout) :: shiftk(3,210) !vz_i
2321 
2322 !Local variables-------------------------------
2323 !scalars
2324  integer,parameter :: kptopt=1,mkpt_list=100000
2325  integer :: ang90,center,dirvacuum,equal,igrid,igrid_current,iholohedry,ii,init_mult,iscale,iscf
2326  integer :: iset,mult1,mult2,mult3,ndims,nkpt,nkpt_current,nkpt_trial,nset
2327  real(dp) :: buffer_scale,determinant,fact,factor,kptrlen_current,kptrlen_max,kptrlen_target
2328  real(dp) :: kptrlen_trial,length1,length2,length3,length_axis1,length_axis2
2329  real(dp) :: length_axis3,merit_factor,mult1h,mult2h,mult3h,reduceda,reducedb
2330  real(dp) :: sca,scb,scc,surface,ucvol
2331  character(len=500) :: message
2332 !arrays
2333  integer :: kptrlatt_current(3,3),kptrlatt_trial(3,3)
2334  integer,allocatable :: grid_list(:)
2335  real(dp) :: axes(3,3),gmet(3,3),gprimd(3,3),matrix1(3,3),matrix2(3,3)
2336  real(dp) :: metmin(3,3),minim(3,3),r2d(3,3),rmet(3,3),rsuper(3,3)
2337  real(dp) :: shiftk_current(3,210),shiftk_trial(3,210)
2338  real(dp),allocatable :: kpt(:,:),kptrlen_list(:),wtk(:)
2339 
2340 ! *************************************************************************
2341 
2342  kptrlen_target=kptrlen
2343 
2344 !The vacuum array must be made of 0 or 1
2345  do ii=1,3
2346    if(vacuum(ii)/=0 .and. vacuum(ii)/=1)then
2347      write(message,'(a,a,a,i1,a,i3,a,a)')&
2348 &     'The values of vacuum must be 0 or 1.',ch10,&
2349 &     'However, the input vacuum(',ii,') is',vacuum(ii),ch10,&
2350 &     'Action : correct vacuum in your input file.'
2351      MSG_ERROR(message)
2352    end if
2353  end do
2354 
2355 !Specific preparation for 2-dimensional system
2356  if(sum(vacuum(:))==1)then
2357 
2358 !  Make the non-active vector orthogonal to the active vectors,
2359 !  and take it along the z direction
2360    if(vacuum(1)==1)then
2361      r2d(1,3)=rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3)
2362      r2d(2,3)=rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3)
2363      r2d(3,3)=rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3)
2364      r2d(:,1)=rprimd(:,2)
2365      r2d(:,2)=rprimd(:,3)
2366      dirvacuum=1
2367    else if(vacuum(2)==1)then
2368      r2d(1,3)=rprimd(2,3)*rprimd(3,1)-rprimd(3,3)*rprimd(2,1)
2369      r2d(2,3)=rprimd(3,3)*rprimd(1,1)-rprimd(1,3)*rprimd(3,1)
2370      r2d(3,3)=rprimd(1,3)*rprimd(2,1)-rprimd(2,3)*rprimd(1,1)
2371      r2d(:,1)=rprimd(:,3)
2372      r2d(:,2)=rprimd(:,1)
2373      dirvacuum=2
2374    else if(vacuum(3)==1)then
2375      r2d(1,3)=rprimd(2,1)*rprimd(3,2)-rprimd(3,1)*rprimd(2,2)
2376      r2d(2,3)=rprimd(3,1)*rprimd(1,2)-rprimd(1,1)*rprimd(3,2)
2377      r2d(3,3)=rprimd(1,1)*rprimd(2,2)-rprimd(2,1)*rprimd(1,2)
2378      r2d(:,1)=rprimd(:,1)
2379      r2d(:,2)=rprimd(:,2)
2380      dirvacuum=3
2381    end if
2382    surface=sqrt(sum(r2d(:,3)**2))
2383 !  Identify the 2-D Bravais lattice
2384 !  DEBUG
2385 !  write(std_out,*)' r2d=',r2d(:,:)
2386 !  ENDDEBUG
2387    call metric(gmet,gprimd,-1,rmet,r2d,ucvol)
2388    call smallprim(metmin,minim,r2d)
2389 !  DEBUG
2390 !  write(std_out,*)' minim=',minim(:,:)
2391 !  ENDDEBUG
2392    ang90=0 ; equal=0 ; center=0
2393    axes(:,:)=minim(:,:)
2394    if(abs(metmin(1,2))<tol8)ang90=1
2395    if(abs(metmin(1,1)-metmin(2,2))<tol8)equal=1
2396    if(ang90==1)then
2397      if(equal==1)iholohedry=4
2398      if(equal==0)iholohedry=2
2399    else if(equal==1)then
2400      reduceda=metmin(1,2)/metmin(1,1)
2401      if(abs(reduceda+0.5_dp)<tol8)then
2402        iholohedry=3
2403      else if(abs(reduceda-0.5_dp)<tol8)then
2404        iholohedry=3
2405 !      Use conventional axes
2406        axes(:,2)=minim(:,2)-minim(:,1)
2407      else
2408        iholohedry=2 ; center=1
2409        axes(:,1)=minim(:,1)+minim(:,2)
2410        axes(:,2)=minim(:,2)-minim(:,1)
2411      end if
2412    else
2413      reduceda=metmin(1,2)/metmin(1,1)
2414      reducedb=metmin(1,2)/metmin(2,2)
2415      if(abs(reduceda+0.5_dp)<tol8)then
2416        iholohedry=2 ; center=1
2417        axes(:,2)=2.0_dp*minim(:,2)+minim(:,1)
2418      else if(abs(reduceda-0.5_dp)<tol8)then
2419        iholohedry=2 ; center=1
2420        axes(:,2)=2.0_dp*minim(:,2)-minim(:,1)
2421      else if(abs(reducedb+0.5_dp)<tol8)then
2422        iholohedry=2 ; center=1
2423        axes(:,1)=2.0_dp*minim(:,1)+minim(:,2)
2424      else if(abs(reducedb-0.5_dp)<tol8)then
2425        iholohedry=2 ; center=1
2426        axes(:,1)=2.0_dp*minim(:,1)-minim(:,2)
2427      else
2428        iholohedry=1
2429      end if
2430    end if
2431 !  Make sure that axes form a right-handed coordinate system
2432    determinant=axes(1,1)*axes(2,2)*axes(3,3) &
2433 &   +axes(1,2)*axes(2,3)*axes(3,1) &
2434 &   +axes(1,3)*axes(3,2)*axes(2,1) &
2435 &   -axes(1,1)*axes(3,2)*axes(2,3) &
2436 &   -axes(1,3)*axes(2,2)*axes(3,1) &
2437 &   -axes(1,2)*axes(2,1)*axes(3,3)
2438    if(determinant<zero)then
2439      axes(:,1)=-axes(:,1)
2440    end if
2441 !  Prefer symmetry axes on the same side as the primitive axes
2442    sca=axes(1,1)*r2d(1,1)+axes(2,1)*r2d(2,1)+axes(3,1)*r2d(3,1)
2443    scb=axes(1,2)*r2d(1,2)+axes(2,2)*r2d(2,2)+axes(3,2)*r2d(3,2)
2444    scc=axes(1,3)*rprimd(1,dirvacuum)&
2445 &   +axes(2,3)*rprimd(2,dirvacuum)&
2446 &   +axes(3,3)*rprimd(3,dirvacuum)
2447    if(sca<-tol8 .and. scb<-tol8)then
2448      axes(:,1)=-axes(:,1) ; sca=-sca
2449      axes(:,2)=-axes(:,2) ; scb=-scb
2450    end if
2451 !  Doing this might change the angle between vectors, so that
2452 !  the cell is not conventional anymore
2453 !  if(sca<-tol8 .and. scc<-tol8)then
2454 !  axes(:,1)=-axes(:,1) ; sca=-sca
2455 !  axes(:,3)=-axes(:,3) ; scc=-scc
2456 !  end if
2457 !  if(scb<-tol8 .and. scc<-tol8)then
2458 !  axes(:,2)=-axes(:,2) ; scb=-scb
2459 !  axes(:,3)=-axes(:,3) ; scc=-scc
2460 !  end if
2461    length_axis1=sqrt(axes(1,1)**2+axes(2,1)**2+axes(3,1)**2)
2462    length_axis2=sqrt(axes(1,2)**2+axes(2,2)**2+axes(3,2)**2)
2463 
2464 !  DEBUG
2465 !  write(std_out,*)' testkgrid : iholohedry, center =',iholohedry,center
2466 !  write(std_out,*)' testkgrid : axis 1=',axes(:,1)
2467 !  write(std_out,*)' testkgrid : axis 2=',axes(:,2)
2468 !  write(std_out,*)' testkgrid : axis 3=',axes(:,3)
2469 !  write(std_out,*)' testkgrid : length_axis=',length_axis1,length_axis2
2470 !  ENDDEBUG
2471 
2472 !  End special treatment of 2-D case
2473  end if
2474 
2475 !3-dimensional system
2476  if(sum(vacuum(:))==0)then
2477    iholohedry=bravais(1)
2478    center=bravais(2)
2479    fact=1.0_dp
2480    if(center/=0)fact=0.5_dp
2481    matrix1(:,1)=bravais(3:5)*fact
2482    matrix1(:,2)=bravais(6:8)*fact
2483    matrix1(:,3)=bravais(9:11)*fact
2484    call matr3inv(matrix1,matrix2)
2485    do ii=1,3
2486      axes(:,ii)=rprimd(:,1)*matrix2(ii,1)+rprimd(:,2)*matrix2(ii,2)+rprimd(:,3)*matrix2(ii,3)
2487    end do
2488    length_axis1=sqrt(axes(1,1)**2+axes(2,1)**2+axes(3,1)**2)
2489    length_axis2=sqrt(axes(1,2)**2+axes(2,2)**2+axes(3,2)**2)
2490    length_axis3=sqrt(axes(1,3)**2+axes(2,3)**2+axes(3,3)**2)
2491 !  DEBUG
2492 !  write(std_out,*)' testkgrid : axes=',axes(:,:)
2493 !  write(std_out,*)' length_axis=',length_axis1,length_axis2,length_axis3
2494 !  ENDDEBUG
2495  end if
2496 
2497 !This routine examine only primitive k lattices.
2498  nshiftk=1
2499 
2500 !If prtkpt/=0, will examine more grids than strictly needed
2501  buffer_scale=one
2502  if(prtkpt/=0)buffer_scale=two
2503 
2504  if(prtkpt/=0)then
2505    write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
2506 &   ' testkgrid : will perform the analysis of a series of k-grids.',ch10,&
2507 &   '  Note that kptopt=1 in this analysis, irrespective of its input value.',ch10,ch10,&
2508 &   ' Grid#    kptrlatt         shiftk         kptrlen       nkpt  iset',ch10
2509    call wrtout(std_out,message,'COLL')
2510    call wrtout(iout,message,'COLL')
2511    ABI_ALLOCATE(grid_list,(mkpt_list))
2512    ABI_ALLOCATE(kptrlen_list,(mkpt_list))
2513    grid_list(:)=0
2514    kptrlen_list(:)=0.0_dp
2515  end if
2516 
2517  if(sum(vacuum(:))==3)then
2518 
2519    kptrlatt(:,:)=0
2520    kptrlatt(1,1)=1
2521    kptrlatt(2,2)=1
2522    kptrlatt(3,3)=1
2523    shiftk(:,1)=0.0_dp
2524    kptrlen=1000.0_dp
2525    nkpt_current=1
2526    igrid_current=1
2527 
2528    if(prtkpt/=0)then
2529      write(message,&
2530 &     '(a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )&
2531 &     '    1  ',kptrlatt(:,1),'  ',shiftk(1,1),'  ',kptrlen,1,1,ch10,&
2532 &     '       ',kptrlatt(:,2),'  ',shiftk(2,1),ch10,&
2533 &     '       ',kptrlatt(:,3),'  ',shiftk(3,1),ch10
2534      call wrtout(std_out,message,'COLL')
2535      call wrtout(iout,message,'COLL')
2536 !    The unit cell volume is fake
2537      ucvol=kptrlen**3
2538    end if
2539 
2540  else
2541 
2542    nkpt=0 ; nkpt_current=0 ; iscf=1 ; iset=1
2543    kptrlen_current=0.0_dp
2544    mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1
2545    ABI_ALLOCATE(kpt,(3,nkpt))
2546    ABI_ALLOCATE(wtk,(nkpt))
2547    call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2548 
2549 !  Loop on different grids, the upper limit is only to avoid an infinite loop
2550    do igrid=1,1000
2551 
2552      kptrlatt_trial(:,:)=0
2553      kptrlatt_trial(1,1)=1
2554      kptrlatt_trial(2,2)=1
2555      kptrlatt_trial(3,3)=1
2556      shiftk_trial(:,1)=0.0_dp
2557 
2558 !    1-dimensional system
2559      if(sum(vacuum(:))==2)then
2560        if(vacuum(1)==0)then
2561          kptrlatt_trial(1,1)=2*igrid ; shiftk_trial(1,1)=0.5_dp
2562        else if(vacuum(2)==0)then
2563          kptrlatt_trial(2,2)=2*igrid ; shiftk_trial(2,1)=0.5_dp
2564        else if(vacuum(3)==0)then
2565          kptrlatt_trial(3,3)=2*igrid ; shiftk_trial(3,1)=0.5_dp
2566        end if
2567      end if
2568 
2569 !    2-dimensional system
2570      if(sum(vacuum(:))==1)then
2571 
2572 !      Treat hexagonal holohedries separately
2573        if(iholohedry==3)then
2574 
2575 !        write(std_out,*)' testkgrid : 2D, hexagonal'
2576 
2577          mult1=mult1+1
2578          nset=4
2579          if(iset==1)then
2580            rsuper(:,1)=axes(:,1)*mult1
2581            rsuper(:,2)=axes(:,2)*mult1
2582            shiftk_trial(:,1)=0.0_dp
2583          else if(iset==2)then
2584            rsuper(:,1)=(axes(:,1)-axes(:,2))  *mult1
2585            rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
2586            shiftk_trial(1,1)=1.0_dp/3.0_dp
2587            shiftk_trial(2,1)=1.0_dp/3.0_dp
2588          else if(iset==3)then
2589            rsuper(:,1)=(axes(:,1)-axes(:,2))  *mult1
2590            rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
2591            shiftk_trial(:,1)=0.0_dp
2592          else if(iset==4)then
2593            rsuper(:,1)=axes(:,1)*mult1
2594            rsuper(:,2)=axes(:,2)*mult1
2595            shiftk_trial(1,1)=0.5_dp
2596            shiftk_trial(2,1)=0.5_dp
2597          end if
2598 
2599        else
2600 !        Now treat all other holohedries
2601          length1=length_axis1*mult1
2602          length2=length_axis2*mult2
2603 !        DEBUG
2604 !        write(std_out,*)' testkgrid : (2d) length=',length1,length2
2605 !        ENDDEBUG
2606          if(abs(length1-length2)<tol8)then
2607            mult1=mult1+1
2608            mult2=mult2+1
2609          else if(length1>length2)then
2610            mult2=mult2+1
2611          else if(length2>length1)then
2612            mult1=mult1+1
2613          end if
2614          nset=4
2615 !        iset==5 and 6 are allowed only for centered lattice
2616          if(center==1)nset=6
2617          if(iset==1 .or. iset==2)then
2618            rsuper(:,1)=axes(:,1)*mult1
2619            rsuper(:,2)=axes(:,2)*mult2
2620          else if(iset==3 .or. iset==4)then
2621            rsuper(:,1)=axes(:,1)*mult1-axes(:,2)*mult2
2622            rsuper(:,2)=axes(:,1)*mult1+axes(:,2)*mult2
2623          else if(iset==5 .or. iset==6)then
2624            rsuper(:,1)=axes(:,1)*(mult1-0.5_dp)-axes(:,2)*(mult2-0.5_dp)
2625            rsuper(:,2)=axes(:,1)*(mult1-0.5_dp)+axes(:,2)*(mult2-0.5_dp)
2626          end if
2627 !        This was the easiest way to code all even mult1 and mult2 pairs :
2628 !        make separate series for this possibility.
2629          if(iset==2 .or. iset==4 .or. iset==6)then
2630            rsuper(:,1)=2.0_dp*rsuper(:,1)
2631            rsuper(:,2)=2.0_dp*rsuper(:,2)
2632          end if
2633          shiftk_trial(1,1)=0.5_dp
2634          shiftk_trial(2,1)=0.5_dp
2635 
2636        end if
2637 
2638 !      Put back the inactive direction
2639        if(dirvacuum==1)then
2640          rsuper(:,3)=rsuper(:,1)
2641          shiftk_trial(3,1)=shiftk_trial(1,1)
2642          rsuper(:,1)=rprimd(:,1)
2643          shiftk_trial(1,1)=0.0_dp
2644        else if(dirvacuum==2)then
2645          rsuper(:,3)=rsuper(:,1)
2646          shiftk_trial(3,1)=shiftk_trial(1,1)
2647          rsuper(:,1)=rsuper(:,2)
2648          shiftk_trial(1,1)=shiftk_trial(2,1)
2649          rsuper(:,2)=rprimd(:,2)
2650          shiftk_trial(2,1)=0.0_dp
2651        else if(dirvacuum==3)then
2652          rsuper(:,3)=rprimd(:,3)
2653          shiftk_trial(3,1)=0.0_dp
2654        end if
2655 
2656 !      The supercell and the corresponding shift have been generated !
2657 !      Convert cartesian coordinates into kptrlatt_trial
2658        do ii=1,3
2659          kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+&
2660 &         gprimd(2,:)*rsuper(2,ii)+&
2661 &         gprimd(3,:)*rsuper(3,ii)  )
2662        end do
2663 
2664 !      End of 2-dimensional system
2665      end if
2666 
2667 !    3-dimensional system
2668      if(sum(vacuum(:))==0)then
2669 !      Treat hexagonal holohedries separately
2670        if(iholohedry==6)then
2671          length1=length_axis1*mult1
2672          length3=length_axis3*mult3
2673 !        DEBUG
2674 !        write(std_out,*)' testkgrid : (hex) lengths=',length1,length2
2675 !        ENDDEBUG
2676          if(abs(length1-length3)<tol8)then
2677            mult1=mult1+1
2678            mult3=mult3+1
2679          else if(length1>length3)then
2680            mult3=mult3+1
2681          else if(length3>length1)then
2682            mult1=mult1+1
2683          end if
2684          nset=4
2685          if(iset==1)then
2686            rsuper(:,1)=axes(:,1)*mult1
2687            rsuper(:,2)=axes(:,2)*mult1
2688            rsuper(:,3)=axes(:,3)*mult3
2689            shiftk_trial(:,1)=0.0_dp
2690            shiftk_trial(3,1)=0.5_dp
2691          else if(iset==2)then
2692            rsuper(:,1)=(axes(:,1)-axes(:,2))  *mult1
2693            rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
2694            rsuper(:,3)=axes(:,3)*mult3
2695            shiftk_trial(1,1)=1.0_dp/3.0_dp
2696            shiftk_trial(2,1)=1.0_dp/3.0_dp
2697            shiftk_trial(3,1)=0.5_dp
2698          else if(iset==3)then
2699            rsuper(:,1)=(axes(:,1)-axes(:,2))  *mult1
2700            rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1
2701            rsuper(:,3)=axes(:,3)*mult3
2702            shiftk_trial(:,1)=0.0_dp
2703            shiftk_trial(3,1)=0.5_dp
2704          else if(iset==4)then
2705            rsuper(:,1)=axes(:,1)*mult1
2706            rsuper(:,2)=axes(:,2)*mult1
2707            rsuper(:,3)=axes(:,3)*mult3
2708            shiftk_trial(:,1)=0.5_dp
2709          end if
2710 
2711        else
2712 !        Now treat all other holohedries
2713          length1=length_axis1*mult1
2714          length2=length_axis2*mult2
2715          length3=length_axis3*mult3
2716 !        DEBUG
2717 !        write(std_out,*)' testkgrid : length=',length1,length2,length3
2718 !        ENDDEBUG
2719          if(length2>length1+tol8 .and. length3>length1+tol8)then
2720            mult1=mult1+1
2721          else if(length1>length2+tol8 .and. length3>length2+tol8)then
2722            mult2=mult2+1
2723          else if(length1>length3+tol8 .and. length2>length3+tol8)then
2724            mult3=mult3+1
2725          else if(abs(length2-length3)<tol8 .and. &
2726 &           abs(length1-length3)<tol8 .and. &
2727 &           abs(length1-length2)<tol8        )then
2728            mult1=mult1+1 ; mult2=mult2+1 ; mult3=mult3+1
2729          else if(abs(length1-length2)<tol8)then
2730            mult1=mult1+1 ; mult2=mult2+1
2731          else if(abs(length1-length3)<tol8)then
2732            mult1=mult1+1 ; mult3=mult3+1
2733          else if(abs(length2-length3)<tol8)then
2734            mult2=mult2+1 ; mult3=mult3+1
2735          end if
2736          nset=6
2737          if(center==-1 .or. center==-3)nset=8
2738          if(iset==1 .or. iset==2)then
2739 !          Simple lattice of k points
2740            rsuper(:,1)=axes(:,1)*mult1
2741            rsuper(:,2)=axes(:,2)*mult2
2742            rsuper(:,3)=axes(:,3)*mult3
2743            shiftk_trial(:,1)=0.5_dp
2744          else if(iset==3 .or. iset==4)then
2745 !          FCC lattice of k points = BCC lattice in real space
2746            rsuper(:,1)=-axes(:,1)*mult1+axes(:,2)*mult2+axes(:,3)*mult3
2747            rsuper(:,2)= axes(:,1)*mult1-axes(:,2)*mult2+axes(:,3)*mult3
2748            rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2-axes(:,3)*mult3
2749            shiftk_trial(:,1)=0.5_dp
2750          else if(iset==5 .or. iset==6)then
2751 !          BCC lattice of k points = FCC lattice in real space
2752            rsuper(:,1)=                 axes(:,2)*mult2+axes(:,3)*mult3
2753            rsuper(:,2)= axes(:,1)*mult1                +axes(:,3)*mult3
2754            rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2
2755 !          The BCC lattice has no empty site with full symmetry
2756            shiftk_trial(:,1)=0.0_dp
2757          else if(iset==7 .or. iset==8)then
2758 !          iset==7 and 8 are allowed only for centered lattice
2759            mult1h=mult1-0.5_dp
2760            mult2h=mult2-0.5_dp
2761            mult3h=mult3-0.5_dp
2762            if(center==-1)then
2763 !            FCC lattice of k points = BCC lattice in real space
2764              rsuper(:,1)=-axes(:,1)*mult1h+axes(:,2)*mult2h+axes(:,3)*mult3h
2765              rsuper(:,2)= axes(:,1)*mult1h-axes(:,2)*mult2h+axes(:,3)*mult3h
2766              rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h-axes(:,3)*mult3h
2767              shiftk_trial(:,1)=0.5_dp
2768            else if(center==-3)then
2769 !            BCC lattice of k points = FCC lattice in real space
2770              rsuper(:,1)=                  axes(:,2)*mult2h+axes(:,3)*mult3h
2771              rsuper(:,2)= axes(:,1)*mult1h                 +axes(:,3)*mult3h
2772              rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h
2773 !            The BCC lattice has no empty site with full symmetry
2774              shiftk_trial(:,1)=0.0_dp
2775            end if
2776          end if
2777 !        This was the easiest way to code all even mult1, mult2, mult3 triplets :
2778 !        make separate series for this possibility.
2779          if(2*(iset/2)==iset)then
2780            rsuper(:,1)=2.0_dp*rsuper(:,1)
2781            rsuper(:,2)=2.0_dp*rsuper(:,2)
2782            rsuper(:,3)=2.0_dp*rsuper(:,3)
2783          end if
2784        end if
2785 
2786 !      DEBUG
2787 !      write(std_out,*)' testkgrid : gprimd=',gprimd(:,:)
2788 !      write(std_out,*)' testkgrid : rsuper=',rsuper(:,:)
2789 !      write(std_out,*)' testkgrid : iset  =',iset
2790 !      ENDDEBUG
2791 
2792 
2793 !      The supercell and the corresponding shift have been generated !
2794 !      Convert cartesian coordinates into kptrlatt_trial
2795        do ii=1,3
2796          kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+&
2797 &         gprimd(2,:)*rsuper(2,ii)+&
2798 &         gprimd(3,:)*rsuper(3,ii)  )
2799        end do
2800 
2801 !      End of 3-dimensional system
2802      end if
2803 
2804 !    DEBUG
2805 !    write(std_out,*)' testkgrid : before getkgrid'
2806 !    write(std_out,*)' testkgrid : rprimd=',rprimd(:,:)
2807 !    write(std_out,*)' testkgrid : kptrlatt_trial=',kptrlatt_trial(:,:)
2808 !    ENDDEBUG
2809 
2810      call getkgrid(0,0,iscf,kpt,&
2811 &     kptopt,kptrlatt_trial,kptrlen_trial,&
2812 &     msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,&
2813 &     shiftk_trial,symafm,symrel,vacuum,wtk)
2814 
2815 !    DEBUG
2816 !    write(std_out,*)' testkgrid : after getkgrid'
2817 !    ENDDEBUG
2818 
2819 !    In case one does not need the full list of grids, will take a shortcut, and go to one of the last grids of the series,
2820 !    that generates a kptrlen_trial that is just below kptrlen.
2821      if(prtkpt==0 .and. init_mult==1 .and. kptrlen_trial<(half-tol8)*kptrlen )then
2822        iscale=int((one-tol8)*kptrlen/kptrlen_trial)
2823        mult1=mult1*iscale
2824        mult2=mult2*iscale
2825        mult3=mult3*iscale
2826        init_mult=0
2827 !       DEBUG
2828 !       write(std_out,*)' testkgrid : iscale=',iscale
2829 !       ENDDEBUG
2830        kptrlatt_trial(:,:)=kptrlatt_trial(:,:)*iscale
2831        call getkgrid(0,0,iscf,kpt,&
2832 &       kptopt,kptrlatt_trial,kptrlen_trial,&
2833 &       msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,&
2834 &       shiftk_trial,symafm,symrel,vacuum,wtk)
2835      end if
2836 
2837      if( (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_current==0) .or. &
2838 &     (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_trial<nkpt_current) .or. &
2839 &     (nkpt_trial==nkpt_current  .and. kptrlen_trial>kptrlen_current*(1.0_dp+tol8)))then
2840 
2841        kptrlatt_current(:,:)=kptrlatt_trial(:,:)
2842        nkpt_current=nkpt_trial
2843        shiftk_current(:,:)=shiftk_trial(:,:)
2844        kptrlen_current=kptrlen_trial
2845        igrid_current=igrid
2846      end if
2847 
2848      if(prtkpt/=0)then
2849        write(message,'(i5,a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )&
2850 &       igrid,'  ',kptrlatt_trial(:,1),'  ',shiftk_trial(1,1),&
2851 &       '  ',kptrlen_trial,nkpt_trial,iset,ch10,&
2852 &       '       ',kptrlatt_trial(:,2),'  ',shiftk_trial(2,1),ch10,&
2853 &       '       ',kptrlatt_trial(:,3),'  ',shiftk_trial(3,1),ch10
2854        call wrtout(std_out,message,'COLL')
2855        call wrtout(iout,message,'COLL')
2856 
2857 !      Keep track of this grid, if it is worth
2858        if(kptrlen_trial > kptrlen_list(nkpt_trial)*(1.0_dp+tol8))then
2859          grid_list(nkpt_trial)=igrid
2860          kptrlen_list(nkpt_trial)=kptrlen_trial
2861        end if
2862      end if
2863 
2864 !    Treat 1-D case
2865      if( sum(vacuum(:))==2 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )exit
2866 
2867 !    Treat 2-D case or 3-D case
2868      if( sum(vacuum(:))<=1 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )then
2869 !      The present set of sets of k points is finished :
2870 !      either it was the last, or one has to go to the next one
2871        if(iset==nset)exit
2872        iset=iset+1
2873        mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1
2874      end if
2875 
2876    end do ! igrid=1,1000
2877 
2878    ABI_DEALLOCATE(kpt)
2879    ABI_DEALLOCATE(wtk)
2880 
2881    kptrlatt(:,:)=kptrlatt_current(:,:)
2882    shiftk(:,:)=shiftk_current(:,:)
2883    kptrlen=kptrlen_current
2884 
2885  end if ! test on the number of dimensions
2886 
2887  if(prtkpt/=0)then
2888 
2889 !  sqrt(1/2) comes from the FCC packing, the best one
2890    factor=sqrt(0.5_dp)/ucvol/dble(nsym)
2891    ndims=3
2892    if(sum(vacuum(:))/=0)then
2893      if(sum(vacuum(:))==1)then
2894 !      sqrt(3/4) comes from the hex packing, the best one
2895 !      one multiplies by 2 because nsym is likely twice the number
2896 !      of symmetries that can be effectively used in 2D
2897        ndims=2 ; factor=sqrt(0.75_dp)/surface/dble(nsym)*2
2898        write(message,'(2a)' )ch10,' Note that the system is bi-dimensional.'
2899      else if(sum(vacuum(:))==2)then
2900        ndims=1 ; factor=1/ucvol
2901        write(message,'(2a)' )ch10,' Note that the system is uni-dimensional.'
2902      else if(sum(vacuum(:))==3)then
2903        ndims=0
2904        write(message,'(2a)' )ch10,' Note that the system is zero-dimensional.'
2905      end if
2906      call wrtout(std_out,message,'COLL')
2907      call wrtout(iout,message,'COLL')
2908    end if
2909 
2910 !  The asymptotic value of the merit factor is determined
2911 !  by the set of symmetries : in 3D, if it includes the
2912 !  inversion symmetry, the limit will be 1, if not, it
2913 !  will be two. In 2D, if it includes the inversion symmetry
2914 !  and an operation that maps z on -z, it will tend to one,
2915 !  while if only one of these operations is present,
2916 !  it will tend to two, and if none is present, it will tend to four.
2917    write(message,'(11a)' )ch10,&
2918 &   ' List of best grids, ordered by nkpt.',ch10,&
2919 &   '  (stop at a value of kptrlen 20% larger than the target value).',ch10,&
2920 &   '  (the merit factor will tend to one or two in 3 dimensions)',ch10,&
2921 &   '  (and to one, two or four in 2 dimensions)',ch10,ch10,&
2922 &   '    nkpt   kptrlen    grid#  merit_factor'
2923    call wrtout(std_out,message,'COLL')
2924    call wrtout(iout,message,'COLL')
2925 
2926    kptrlen_max=0.0_dp
2927    do ii=1,mkpt_list
2928      if(kptrlen_list(ii)>kptrlen_max*(1.0_dp+tol8))then
2929        kptrlen_max=kptrlen_list(ii)
2930        merit_factor=kptrlen_max**ndims/dble(ii)*factor
2931        write(message, '(i6,es14.4,i6,f12.4)' )ii,kptrlen_max,grid_list(ii),merit_factor
2932        call wrtout(std_out,message,'COLL')
2933        call wrtout(iout,message,'COLL')
2934      end if
2935      if(kptrlen_max>1.2_dp*(1.0_dp-tol8)*kptrlen_target)exit
2936    end do
2937 
2938    write(message,'(a,a,es14.4,a,a,i6,a,a,a,es14.4,a,i6)' )ch10,&
2939 &   ' For target kptrlen=',kptrlen_target,',',&
2940 &   ' the selected grid is number',igrid_current,',',ch10,&
2941 &   '     giving kptrlen=',kptrlen_current,' with nkpt=',nkpt_current
2942    call wrtout(std_out,message,'COLL')
2943    call wrtout(iout,message,'COLL')
2944 
2945    write(message,'(a,a,a,a)' )ch10,&
2946 &   ' testkgrid : stop after analysis of a series of k-grids.',ch10,&
2947 &   '  For usual production runs, set prtkpt back to 0 (the default).'
2948    call wrtout(std_out,message,'COLL',do_flush=.True.)
2949    call wrtout(iout,message,'COLL',do_flush=.True.)
2950 
2951    call abi_abort('PERS',exit_status=0,print_config=.false.)
2952  end if
2953 
2954 end subroutine testkgrid

m_kpts/tetra_from_kptrlatt [ Functions ]

[ Top ] [ m_kpts ] [ Functions ]

NAME

 tetra_from_kptrlatt

FUNCTION

  Create an instance of `t_tetrahedron` from kptrlatt and shiftk

INPUTS

  cryst<cryst_t>=Crystalline structure.
  kptopt=Option for the k-point generation.
  kptrlatt(3,3)=k-point lattice specification
  nshiftk= number of shift vectors.
  shiftk(3,nshiftk)=shift vectors for k point generation
  nkibz=Number of points in the IBZ
  kibz(3,nkibz)=Reduced coordinates of the k-points in the IBZ.

OUTPUT

  tetra<t_tetrahedron>=Tetrahedron object, fully initialized if ierr == 0.
  msg=Error message if ierr /= 0
  ierr=Exit status

PARENTS

      gstate,wfk_analyze

CHILDREN

      init_tetra,listkk,smpbz

SOURCE

254 type(t_tetrahedron) function tetra_from_kptrlatt( &
255 &  cryst, kptopt, kptrlatt, nshiftk, shiftk, nkibz, kibz, msg, ierr) result (tetra)
256 
257 
258 !This section has been created automatically by the script Abilint (TD).
259 !Do not modify the following lines by hand.
260 #undef ABI_FUNC
261 #define ABI_FUNC 'tetra_from_kptrlatt'
262 !End of the abilint section
263 
264  implicit none
265 
266 !Arguments ------------------------------------
267 !scalars
268  integer,intent(in) :: kptopt,nshiftk,nkibz
269  integer,intent(out) :: ierr
270  character(len=*),intent(out) :: msg
271  type(crystal_t),intent(in) :: cryst
272 !arrays
273  integer,intent(in) :: kptrlatt(3,3)
274  real(dp),intent(in) :: shiftk(3,nshiftk),kibz(3,nkibz)
275 
276 !Local variables-------------------------------
277 !scalars
278  integer :: nkfull,timrev,sppoldbl,my_nkibz,new_nshiftk
279  real(dp) :: dksqmax
280  character(len=80) :: errorstring
281 !arrays
282  integer :: new_kptrlatt(3,3)
283  integer,allocatable :: indkk(:,:)
284  real(dp) :: rlatt(3,3),klatt(3,3)
285  real(dp),allocatable :: kfull(:,:),my_kibz(:,:),my_wtk(:),new_shiftk(:,:)
286 
287 ! *************************************************************************
288 
289  ierr = 0
290 
291  ! Refuse only 1 kpoint: the algorithms are no longer valid. DOH!
292  if (nkibz == 1) then
293    msg = 'You need at least 2 kpoints to use the tetrahedron method.'
294    ierr = 1; goto 10
295  end if
296  if (all(kptrlatt == 0)) then
297    msg = 'Cannot generate tetrahedron because input kptrlatt == 0.'
298    ierr = 1; goto 10
299  end if
300  if (kptopt <= 0) then
301    msg = sjoin("Cannot generate tetrahedron because input kptopt:", itoa(kptopt))
302    ierr = 1; goto 10
303  end if
304 
305  call kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, &
306    my_nkibz, my_kibz, my_wtk, nkfull, kfull, new_kptrlatt=new_kptrlatt, new_shiftk=new_shiftk)
307 
308  ABI_FREE(my_kibz)
309  ABI_FREE(my_wtk)
310  new_nshiftk = size(new_shiftk, dim=2)
311 
312  if (my_nkibz /= nkibz) then
313    msg = sjoin("Input nkibz:", itoa(nkibz), "does not agree with computed value:", itoa(my_nkibz))
314    ierr = 1; goto 10
315  end if
316 
317  ! Do not support new_nshiftk > 1: lattice must be decomposed into boxes
318  ! and this is not always possible (I think) with bizzare shiftks
319  ! normally at this point we have incorporated everything into
320  ! new_kptrlatt, and only 1 shift is needed (in particular for MP grids).
321  if (new_nshiftk > 1) then
322    write(msg, "(9a)") &
323      'Cannot create tetrahedron object...',ch10, &
324      'Only simple lattices are supported. Action: use nshiftk=1.',ch10, &
325      'new_shiftk: ', trim(ltoa(reshape(new_shiftk, [3*new_nshiftk]))),ch10, &
326      'new_kptrlatt: ', trim(ltoa(reshape(new_kptrlatt, [9])))
327    ierr = 2; goto 10
328  end if
329 
330  ! Costruct full BZ and create mapping BZ --> IBZ
331  ! Note:
332  !   - we don't change the value of nsppol hence sppoldbl is set to 1
333  !   - we use symrec (operations in reciprocal space)
334  !
335  sppoldbl = 1; timrev = kpts_timrev_from_kptopt(kptopt)
336  ABI_MALLOC(indkk, (nkfull*sppoldbl,6))
337 
338  ! Compute k points from input file closest to the output file
339  call listkk(dksqmax,cryst%gmet,indkk,kibz,kfull,nkibz,nkfull,cryst%nsym,&
340     sppoldbl,cryst%symafm,cryst%symrec,timrev,use_symrec=.True.)
341 
342  if (dksqmax > tol12) then
343    write(msg, '(3a,es16.6,6a)' )&
344    'At least one of the k points could not be generated from a symmetrical one.',ch10,&
345    'dksqmax=',dksqmax,ch10,&
346    'new_kptrkatt= ',trim(ltoa(reshape(new_kptrlatt, [9]))),ch10,&
347    'new_shiftk= ',trim(ltoa(reshape(new_shiftk, [3*new_nshiftk])))
348    ierr = 2; goto 10
349  end if
350 
351  rlatt = new_kptrlatt; call matr3inv(rlatt, klatt)
352 
353  call init_tetra(indkk(:,1), cryst%gprimd, klatt, kfull, nkfull, tetra, ierr, errorstring)
354  if (ierr /= 0) msg = errorstring
355 
356  10 continue
357  if (allocated(indkk)) then
358    ABI_FREE(indkk)
359  end if
360  if (allocated(kfull)) then
361    ABI_FREE(kfull)
362  end if
363  if (allocated(new_shiftk)) then
364    ABI_FREE(new_shiftk)
365  end if
366 
367 end function tetra_from_kptrlatt