TABLE OF CONTENTS
- ABINIT/m_kpts
- m_kpts/get_full_kgrid
- m_kpts/get_kpt_fullbz
- m_kpts/getkgrid
- m_kpts/getkgrid_low
- m_kpts/kpts_ibz_from_kptrlatt
- m_kpts/kpts_map
- m_kpts/kpts_map_print
- m_kpts/kpts_pack_in_stars
- m_kpts/kpts_sort
- m_kpts/kpts_timrev_from_kptopt
- m_kpts/listkk
- m_kpts/mknormpath
- m_kpts/smpbz
- m_kpts/symkchk
- m_kpts/testkgrid
- m_kpts/tetra_from_kptrlatt
ABINIT/m_kpts [ Modules ]
NAME
m_kpts
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 #include "abi_common.h" 20 21 module m_kpts 22 23 use defs_basis 24 use m_errors 25 use m_abicore 26 use m_crystal 27 use m_sort 28 use m_krank 29 use m_htetra 30 use m_xmpi 31 32 use m_time, only : timab, cwtime, cwtime_report 33 use m_copy, only : alloc_copy 34 use m_symtk, only : mati3inv, mati3det, matr3inv, smallprim 35 use m_fstrings, only : sjoin, itoa, ftoa, ltoa, ktoa 36 use m_numeric_tools, only : wrap2_pmhalf 37 use m_geometry, only : metric 38 use m_symkpt, only : symkpt, symkpt_new 39 40 implicit none 41 42 private 43 44 public :: kpts_timrev_from_kptopt ! Returns the value of timrev from kptopt 45 public :: kpts_ibz_from_kptrlatt ! Determines the IBZ, the weights and the BZ from kptrlatt 46 public :: tetra_from_kptrlatt ! Create an instance of htetra_t from kptrlatt and shiftk 47 public :: symkchk ! Checks that the set of k points has the full space group symmetry, 48 ! modulo time reversal if appropriate. 49 public :: kpts_sort ! Order list of k-points according to the norm. 50 public :: kpts_pack_in_stars ! 51 public :: kpts_map ! Compute symmetry table. 52 public :: kpts_map_print ! Print the symmetry table bz2ibz to a list of units with header. 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 private :: get_kpt_fullbz ! Create full grid of kpoints from kptrlatt and shiftk 60 61 public :: smpbz ! Generate a set of special k (or q) points which samples in a homogeneous way the BZ 62 public :: testkgrid ! Test different grids of k points. 63 64 ! FIXME: deprecated 65 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 implementation always assumes kptopt==1 !!!! TODO: This routine should be removed
SOURCE
1742 subroutine get_full_kgrid(indkpt,kpt,kpt_fullbz,kptrlatt,nkpt,nkpt_fullbz,nshiftk,nsym,shiftk,symrel) 1743 1744 !Arguments ------------------------------------ 1745 !scalars 1746 integer,intent(in) :: nkpt,nkpt_fullbz,nshiftk,nsym 1747 !arrays 1748 integer,intent(in) :: kptrlatt(3,3),symrel(3,3,nsym) 1749 integer,intent(out) :: indkpt(nkpt_fullbz) 1750 real(dp),intent(in) :: kpt(3,nkpt),shiftk(3,nshiftk) 1751 real(dp),intent(out) :: kpt_fullbz(3,nkpt_fullbz) 1752 1753 !Local variables------------------------------- 1754 !scalars 1755 integer :: ikpt,isym,itim,timrev 1756 integer :: symrankkpt 1757 character(len=500) :: msg 1758 type(krank_t) :: krank 1759 !arrays 1760 integer :: inv_symrel(3,3,nsym) 1761 real(dp) :: k2(3) 1762 1763 ! ********************************************************************* 1764 1765 !Invert symrels => gives symrels for kpoints 1766 1767 do isym=1,nsym 1768 call mati3inv (symrel(:,:,isym),inv_symrel(:,:,isym)) 1769 end do 1770 1771 call get_kpt_fullbz(kpt_fullbz,kptrlatt,nkpt_fullbz,nshiftk,shiftk) 1772 1773 ! make full k-point rank arrays 1774 krank = krank_new(nkpt, kpt) 1775 1776 !find equivalence to irred kpoints in kpt 1777 indkpt(:) = 0 1778 timrev=1 ! includes the time inversion symmetry 1779 do ikpt=1,nkpt_fullbz 1780 do isym=1,nsym 1781 do itim=1,(1-2*timrev),-2 1782 1783 k2(:) = itim*(inv_symrel(:,1,isym)*kpt_fullbz(1,ikpt) + & 1784 inv_symrel(:,2,isym)*kpt_fullbz(2,ikpt) + & 1785 inv_symrel(:,3,isym)*kpt_fullbz(3,ikpt)) 1786 1787 symrankkpt = krank%get_rank(k2) 1788 if (krank%invrank(symrankkpt) /= -1) indkpt(ikpt) = krank%invrank(symrankkpt) 1789 1790 end do ! loop time reversal symmetry 1791 end do ! loop sym ops 1792 1793 if (indkpt(ikpt) == 0) then 1794 write(msg,'(a,i0)')' indkpt(ikpt) is still 0: no irred kpoint is equiv to ikpt ',ikpt 1795 ABI_BUG(msg) 1796 end if 1797 end do ! loop full kpts 1798 1799 call krank%free() 1800 1801 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
SOURCE
1822 subroutine get_kpt_fullbz(kpt_fullbz,kptrlatt,nkpt_fullbz,nshiftk,shiftk) 1823 1824 !Arguments ------------------------------------ 1825 !scalars 1826 integer,intent(in) :: nkpt_fullbz,nshiftk 1827 !arrays 1828 integer,intent(in) :: kptrlatt(3,3) 1829 real(dp),intent(in) :: shiftk(3,nshiftk) 1830 real(dp),intent(out) :: kpt_fullbz(3,nkpt_fullbz) 1831 1832 !Local variables------------------------------- 1833 !scalars 1834 integer, parameter :: max_number_of_prime=47 1835 integer :: det,ii,ikshft,iprim,jj,kk,nn 1836 character(len=500) :: msg 1837 !arrays 1838 integer :: boundmax(3),boundmin(3),common_factor(3) 1839 integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,& 1840 & 29,31,37,41,43, 47,53,59,61,67,& 1841 & 71,73,79,83,89, 97,101,103,107,109,& 1842 & 113,127,131,137,139, 149,151,157,163,167,& 1843 & 173,179,181,191,193, 197,199/) 1844 real(dp) :: k1(3),k2(3),klatt(3,3),rlatt(3,3),shift(3),test_rlatt(3,3) 1845 1846 ! ********************************************************************* 1847 1848 !Identify first factors that can be used to rescale the three kptrlatt vectors 1849 !Only test a large set of prime factors, though ... 1850 do jj=1,3 1851 common_factor(jj)=1 1852 rlatt(:,jj)=kptrlatt(:,jj) 1853 do iprim=1,max_number_of_prime 1854 test_rlatt(:,jj)=rlatt(:,jj)/dble(prime_factor(iprim)) 1855 ! If one of the components is lower than 1 in absolute value, then it is not worth to continue the search. 1856 if(minval(abs(abs(test_rlatt(:,jj))-half))<half-tol8)exit 1857 do 1858 if(sum(abs(test_rlatt(:,jj)-nint(test_rlatt(:,jj)) ))<tol8)then 1859 common_factor(jj)=prime_factor(iprim)*common_factor(jj) 1860 rlatt(:,jj)=rlatt(:,jj)/dble(prime_factor(iprim)) 1861 test_rlatt(:,jj)=test_rlatt(:,jj)/dble(prime_factor(iprim)) 1862 else 1863 exit 1864 end if 1865 end do 1866 end do 1867 end do 1868 call mati3det(kptrlatt,det) 1869 det=det/(common_factor(1)*common_factor(2)*common_factor(3)) 1870 1871 rlatt(:,:)=kptrlatt(:,:) 1872 call matr3inv(rlatt,klatt) 1873 !Now, klatt contains the three primitive vectors of the k lattice, 1874 !in reduced coordinates. One builds all k vectors that 1875 !are contained in the first Brillouin zone, with coordinates 1876 !in the interval [0,1[ . First generate boundaries of a big box. 1877 !In order to generate all possible vectors in the reciprocal space, 1878 !one must consider all multiples of the primitive ones, until a vector with only integers is found. 1879 !The maximum bound is the scale of the corresponding kptrlatt vector, times the determinant of kptrlatt. Also consider negative vectors. 1880 !On this basis, compute the bounds. 1881 do jj=1,3 1882 ! To accomodate the shifts, boundmin starts from -1 1883 ! Well, this is not a complete solution ... 1884 boundmin(jj)=-1-common_factor(jj)*abs(det) 1885 boundmax(jj)=common_factor(jj)*abs(det) 1886 end do 1887 1888 nn=1 1889 do kk=boundmin(3),boundmax(3) 1890 do jj=boundmin(2),boundmax(2) 1891 do ii=boundmin(1),boundmax(1) 1892 do ikshft=1,nshiftk 1893 1894 ! Coordinates of the trial k point with respect to the k primitive lattice 1895 k1(1)=ii+shiftk(1,ikshft) 1896 k1(2)=jj+shiftk(2,ikshft) 1897 k1(3)=kk+shiftk(3,ikshft) 1898 1899 ! Reduced coordinates of the trial k point 1900 k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3) 1901 1902 ! Eliminate the point if outside [0,1[ 1903 if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle 1904 if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle 1905 if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle 1906 1907 ! Wrap the trial values in the interval ]-1/2,1/2] . 1908 call wrap2_pmhalf(k2(1),k1(1),shift(1)) 1909 call wrap2_pmhalf(k2(2),k1(2),shift(2)) 1910 call wrap2_pmhalf(k2(3),k1(3),shift(3)) 1911 if(nn > nkpt_fullbz) then 1912 write (msg,'(a,i0)')' nkpt_fullbz mis-estimated, exceed nn=',nn 1913 ABI_BUG(msg) 1914 end if 1915 kpt_fullbz(:,nn)=k1(:) 1916 nn=nn+1 1917 end do 1918 end do 1919 end do 1920 end do 1921 nn = nn-1 1922 1923 if (nn /= nkpt_fullbz) then 1924 write (msg,'(2(a,i0),a,a)')' nkpt_fullbz= ',nkpt_fullbz,' underestimated nn=',nn,& 1925 & ch10, "Perhaps your k grid or shifts do not correspond to the symmetry?" 1926 ABI_BUG(msg) 1927 end if 1928 1929 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,MAX_NSHIFTK)=shift vectors for k point generation [nkpthf] = number of k points in the full BZ, for the Fock operator.
SOURCE
1156 subroutine getkgrid(chksymbreak,iout,iscf,kpt,kptopt,kptrlatt,kptrlen,& 1157 & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk,& 1158 & fullbz,nkpthf,kpthf,downsampling) ! optional 1159 1160 !Arguments ------------------------------------ 1161 !scalars 1162 integer,intent(in) :: chksymbreak,iout,iscf,kptopt,msym,nkpt,nsym 1163 integer,intent(inout),optional :: nkpthf 1164 integer,intent(inout) :: nshiftk 1165 integer,intent(inout) :: nkpt_computed !vz_i 1166 real(dp),intent(out) :: kptrlen 1167 !arrays 1168 integer,intent(in) :: symafm(msym),symrel(3,3,msym),vacuum(3) 1169 integer,optional,intent(in) :: downsampling(3) 1170 integer,intent(inout) :: kptrlatt(3,3) 1171 integer,allocatable :: indkpt(:) 1172 integer,allocatable :: bz2ibz_smap(:,:) 1173 real(dp),intent(in) :: rprimd(3,3) 1174 real(dp),intent(inout) :: shiftk(3,MAX_NSHIFTK) 1175 real(dp),intent(inout) :: kpt(3,nkpt) !vz_i 1176 real(dp),intent(inout) :: wtk(nkpt) 1177 real(dp),optional,allocatable,intent(out) :: fullbz(:,:) 1178 real(dp),optional,intent(out) :: kpthf(:,:) 1179 1180 !Local variables------------------------------- 1181 real(dp),allocatable :: kpt_tmp(:,:), wtk_tmp(:) 1182 1183 call getkgrid_low(chksymbreak,iout,iscf,kpt_tmp,kptopt,kptrlatt,kptrlen,& 1184 msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk_tmp,indkpt,bz2ibz_smap,& 1185 fullbz,nkpthf,kpthf,downsampling) 1186 1187 if (nkpt > 0) then 1188 kpt(:,1:nkpt) = kpt_tmp(:,1:nkpt) 1189 wtk(1:nkpt) = wtk_tmp(1:nkpt) 1190 end if 1191 1192 ABI_SFREE(kpt_tmp) 1193 ABI_SFREE(wtk_tmp) 1194 ABI_SFREE(indkpt) 1195 ABI_SFREE(bz2ibz_smap) 1196 1197 end subroutine getkgrid
m_kpts/getkgrid_low [ Functions ]
[ Top ] [ m_kpts ] [ Functions ]
NAME
getkgrid_low
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. bz2ibz_smap(nkbz, 6)= Mapping BZ --> IBZ. [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,MAX_NSHIFTK)=shift vectors for k point generation [nkpthf] = number of k points in the full BZ, for the Fock operator.
SOURCE
1246 subroutine getkgrid_low(chksymbreak,iout,iscf,kpt,kptopt,kptrlatt,kptrlen,& 1247 & msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,shiftk,symafm,symrel,vacuum,wtk,indkpt,bz2ibz_smap,& 1248 & fullbz,nkpthf,kpthf,downsampling) ! optional 1249 1250 !Arguments ------------------------------------ 1251 !scalars 1252 integer,intent(in) :: chksymbreak,iout,iscf,kptopt,msym,nkpt,nsym 1253 integer,intent(inout),optional :: nkpthf 1254 integer,intent(inout) :: nshiftk 1255 integer,intent(inout) :: nkpt_computed !vz_i 1256 real(dp),intent(out) :: kptrlen 1257 !arrays 1258 integer,intent(in) :: symafm(msym),symrel(3,3,msym),vacuum(3) 1259 integer,optional,intent(in) :: downsampling(3) 1260 integer,intent(inout) :: kptrlatt(3,3) 1261 real(dp),intent(in) :: rprimd(3,3) 1262 real(dp),intent(inout) :: shiftk(3,MAX_NSHIFTK) 1263 integer,allocatable,intent(out) :: indkpt(:) 1264 integer,allocatable,intent(out) :: bz2ibz_smap(:,:) 1265 real(dp),allocatable,intent(out) :: kpt(:,:) !vz_i 1266 real(dp),allocatable,intent(out) :: wtk(:) 1267 real(dp),optional,allocatable,intent(out) :: fullbz(:,:) 1268 real(dp),optional,intent(out) :: kpthf(:,:) 1269 1270 !Local variables------------------------------- 1271 !scalars 1272 integer, parameter :: max_number_of_prime=47 1273 integer :: brav,decreased,found,ii,ikpt,iprime,ishiftk,isym,jshiftk,kshiftk,mkpt,mult 1274 integer :: nkpthf_computed,nkpt_fullbz,nkptlatt,nshiftk2,nsym_used,option 1275 integer :: test_prime,timrev 1276 integer :: nkpt_use 1277 real(dp) :: length2,ucvol,ucvol_super 1278 character(len=500) :: msg 1279 !arrays 1280 integer, parameter :: prime_factor(max_number_of_prime)=(/2,3,5,7,9, 11,13,17,19,23,& 1281 & 29,31,37,41,43, 47,53,59,61,67,& 1282 & 71,73,79,83,89, 97,101,103,107,109,& 1283 & 113,127,131,137,139, 149,151,157,163,167,& 1284 & 173,179,181,191,193, 197,199/) 1285 integer :: kptrlatt2(3,3) 1286 integer,allocatable :: belong_chain(:),generator(:),number_in_chain(:) 1287 integer,allocatable :: repetition_factor(:),symrec(:,:,:) 1288 ! real(dp) :: cart(3,3) 1289 real(dp) :: dijk(3),delta_dmult(3),dmult(3),fact_vacuum(3),gmet(3,3) 1290 real(dp) :: gmet_super(3,3),gprimd(3,3),gprimd_super(3,3),klatt2(3,3) 1291 real(dp) :: klatt3(3,3),kptrlattr(3,3),ktransf(3,3),ktransf_invt(3,3) 1292 real(dp) :: metmin(3,3),minim(3,3),rmet(3,3),rmet_super(3,3),rprimd_super(3,3) 1293 real(dp),allocatable :: deltak(:,:),kpt_fullbz(:,:),shiftk2(:,:),shiftk3(:,:),spkpt(:,:),wtk_folded(:),wtk_fullbz(:) 1294 1295 ! ************************************************************************* 1296 1297 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 1298 1299 !call cwtime(cpu, wall, gflops, "start") 1300 if (kptopt==1.or.kptopt==4) then 1301 ! Cannot use antiferromagnetic symmetry operations to decrease the number of k points 1302 !XG20191123: now, antiferromagnetic symmetry operations can be used to decrease the number of k points for kptopt==4 1303 nsym_used=0 1304 do isym=1,nsym 1305 if(symafm(isym)==1 .or. kptopt==4)nsym_used=nsym_used+1 1306 end do 1307 ABI_MALLOC(symrec,(3,3,nsym_used)) 1308 nsym_used=0 1309 do isym=1,nsym ! Get the symmetry matrices in terms of reciprocal basis 1310 if(symafm(isym)==1 .or. kptopt==4)then 1311 nsym_used=nsym_used+1 1312 call mati3inv(symrel(:,:,isym),symrec(:,:,nsym_used)) 1313 end if 1314 end do 1315 else if (kptopt==2) then 1316 !Use only the time-reversal 1317 nsym_used=1 1318 ABI_MALLOC(symrec,(3,3,1)) 1319 symrec(1:3,1:3,1)=0 1320 do ii=1,3 1321 symrec(ii,ii,1)=1 1322 end do 1323 end if 1324 1325 kptrlatt2(:,:)=kptrlatt(:,:) 1326 nshiftk2=nshiftk 1327 ABI_MALLOC(shiftk2,(3,MAX_NSHIFTK)) 1328 ABI_MALLOC(shiftk3,(3,MAX_NSHIFTK)) 1329 shiftk2(:,:)=shiftk(:,:) 1330 1331 !Find a primitive k point lattice, if possible, by decreasing the number of shifts. 1332 if(nshiftk2/=1)then 1333 1334 do 1335 ! Loop to be repeated if there has been a successful reduction of nshiftk2 1336 ABI_MALLOC(deltak,(3,nshiftk2)) 1337 ABI_MALLOC(repetition_factor,(nshiftk2)) 1338 ABI_MALLOC(generator,(nshiftk2)) 1339 ABI_MALLOC(belong_chain,(nshiftk2)) 1340 ABI_MALLOC(number_in_chain,(nshiftk2)) 1341 1342 decreased=0 1343 deltak(1,1:nshiftk2)=shiftk2(1,1:nshiftk2)-shiftk2(1,1) 1344 deltak(2,1:nshiftk2)=shiftk2(2,1:nshiftk2)-shiftk2(2,1) 1345 deltak(3,1:nshiftk2)=shiftk2(3,1:nshiftk2)-shiftk2(3,1) 1346 deltak(:,:)=deltak(:,:)-floor(deltak(:,:)+tol8) 1347 1348 ! Identify for each shift, the smallest repetition prime factor that yields a reciprocal lattice vector. 1349 repetition_factor(:)=0 1350 repetition_factor(1)=1 1351 do ishiftk=2,nshiftk2 1352 do iprime=1,max_number_of_prime 1353 test_prime=prime_factor(iprime) 1354 dmult(:)=test_prime*deltak(:,ishiftk) 1355 if(sum(abs( dmult(:)-nint(dmult(:)) ))<tol8)then 1356 repetition_factor(ishiftk)=test_prime 1357 exit 1358 end if 1359 end do 1360 end do 1361 1362 ! Initialize the selection of tentative generators 1363 generator(:)=1 1364 do ishiftk=1,nshiftk2 1365 if(repetition_factor(ishiftk)==0 .or. repetition_factor(ishiftk)==1)generator(ishiftk)=0 1366 end do 1367 1368 ! Try different shifts as generators, by order of increasing repetition factor, 1369 ! provided they are equal or bigger than 2 1370 do iprime=1,max_number_of_prime 1371 do ishiftk=2,nshiftk2 1372 ! Note that ishiftk=1 is never a generator. It is the reference starting point. 1373 if(generator(ishiftk)==1 .and. repetition_factor(ishiftk)==prime_factor(iprime))then 1374 ! Test the generator : is it indeed closed ? 1375 if(prime_factor(iprime)/=2)then 1376 do mult=2,prime_factor(iprime)-1 1377 dmult(:)=mult*deltak(:,ishiftk) 1378 found=0 1379 do jshiftk=1,nshiftk2 1380 delta_dmult(:)=deltak(:,jshiftk)-dmult(:) 1381 if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then 1382 found=1 1383 exit 1384 end if 1385 end do 1386 if(found==0)exit 1387 end do 1388 if(found==0)generator(ishiftk)=0 1389 end if 1390 if(generator(ishiftk)==0)cycle 1391 else 1392 cycle 1393 end if 1394 ! Now, test whether all k points can be found in all possible chains 1395 belong_chain(:)=0 1396 do jshiftk=1,nshiftk2 1397 ! Initialize a chain starting from a k point not yet in a chain 1398 if(belong_chain(jshiftk)==0)then 1399 number_in_chain(:)=0 ! Not a member of the chain (yet) 1400 number_in_chain(jshiftk)=1 ! The first point in chain 1401 do mult=1,prime_factor(iprime)-1 1402 dmult(:)=mult*deltak(:,ishiftk) 1403 found=0 1404 do kshiftk=jshiftk+1,nshiftk2 1405 delta_dmult(:)=deltak(:,kshiftk)-deltak(:,jshiftk)-dmult(:) 1406 if(sum(abs(delta_dmult(:)-nint(delta_dmult(:)) ))<tol8)then 1407 found=1 1408 number_in_chain(kshiftk)=mult+1 1409 exit 1410 end if 1411 end do 1412 if(found==0)then 1413 generator(ishiftk)=0 1414 exit 1415 end if 1416 end do 1417 if(generator(ishiftk)==1)then 1418 ! Store the chain 1419 do kshiftk=1,nshiftk2 1420 if(number_in_chain(kshiftk)/=0)belong_chain(kshiftk)=number_in_chain(kshiftk) 1421 end do 1422 else 1423 exit 1424 end if 1425 end if 1426 end do 1427 1428 if(generator(ishiftk)==0)cycle 1429 1430 ! For the generator based on ishiftk, all the k points have been found to belong to one chain. 1431 ! All the initializing k points in the different chains have belong_chain(:)=1 . 1432 ! They must be kept, and the others thrown away. 1433 ktransf(:,:)=0.0_dp 1434 ktransf(1,1)=1.0_dp 1435 ktransf(2,2)=1.0_dp 1436 ktransf(3,3)=1.0_dp 1437 ! Replace one of the unit vectors by the shift vector deltak(:,ishiftk). 1438 ! However, must pay attention not to make linear combinations. 1439 ! Also, choose positive sign for first-non-zero value. 1440 if(abs(deltak(1,ishiftk)-nint(deltak(1,ishiftk)))>tol8)then 1441 if(deltak(1,ishiftk)>0)ktransf(:,1)= deltak(:,ishiftk) 1442 if(deltak(1,ishiftk)<0)ktransf(:,1)=-deltak(:,ishiftk) 1443 else if(abs(deltak(2,ishiftk)-nint(deltak(2,ishiftk)))>tol8)then 1444 if(deltak(2,ishiftk)>0)ktransf(:,2)= deltak(:,ishiftk) 1445 if(deltak(2,ishiftk)<0)ktransf(:,2)=-deltak(:,ishiftk) 1446 else if(abs(deltak(3,ishiftk)-nint(deltak(3,ishiftk)))>tol8)then 1447 if(deltak(3,ishiftk)>0)ktransf(:,3)= deltak(:,ishiftk) 1448 if(deltak(3,ishiftk)<0)ktransf(:,3)=-deltak(:,ishiftk) 1449 end if 1450 ! Copy the integers to real(dp) 1451 kptrlattr(:,:)=kptrlatt2(:,:) 1452 ! Go to reciprocal space 1453 call matr3inv(kptrlattr,klatt2) 1454 ! Make the transformation 1455 do ii=1,3 1456 klatt3(:,ii)=ktransf(1,ii)*klatt2(:,1)+ktransf(2,ii)*klatt2(:,2)+ktransf(3,ii)*klatt2(:,3) 1457 end do 1458 ! Back to real space 1459 call matr3inv(klatt3,kptrlattr) 1460 ! real(dp) to integer 1461 kptrlatt2(:,:)=nint(kptrlattr(:,:)) 1462 ! Prepare the transformation of the shifts 1463 call matr3inv(ktransf,ktransf_invt) 1464 decreased=1 1465 kshiftk=0 1466 do jshiftk=1,nshiftk2 1467 if(belong_chain(jshiftk)==1)then 1468 kshiftk=kshiftk+1 1469 ! Place the shift with index jshiftk in place of the one in kshiftk, 1470 ! also transform the shift from the old to the new coordinate system 1471 shiftk3(:,kshiftk)=ktransf_invt(1,:)*shiftk2(1,jshiftk)+& 1472 & ktransf_invt(2,:)*shiftk2(2,jshiftk)+& 1473 & ktransf_invt(3,:)*shiftk2(3,jshiftk) 1474 end if 1475 end do 1476 nshiftk2=nshiftk2/prime_factor(iprime) 1477 shiftk2(:,1:nshiftk2)=shiftk3(:,1:nshiftk2)-floor(shiftk3(:,1:nshiftk2)+tol8) 1478 if(kshiftk/=nshiftk2)then 1479 ABI_BUG('The search for a primitive k point lattice contains a bug.') 1480 end if 1481 1482 ! If this trial shift was successful, must exit the loop on trial ishiftk, 1483 ! and reinitialize the global loop 1484 if(decreased==1)exit 1485 end do ! ishiftk 1486 if(decreased==1)exit 1487 end do ! iprime 1488 1489 ABI_FREE(belong_chain) 1490 ABI_FREE(deltak) 1491 ABI_FREE(number_in_chain) 1492 ABI_FREE(repetition_factor) 1493 ABI_FREE(generator) 1494 1495 if(decreased==0 .or. nshiftk2==1)exit 1496 1497 end do ! Infinite loop 1498 1499 end if ! End nshiftk being 1 or larger 1500 1501 !Impose shiftk coordinates to be in [0,1[ 1502 do ishiftk=1,nshiftk2 1503 do ii=1,3 1504 if(shiftk2(ii,ishiftk)>one-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)-1.0_dp 1505 if(shiftk2(ii,ishiftk)<-tol8) shiftk2(ii,ishiftk)=shiftk2(ii,ishiftk)+1.0_dp 1506 end do 1507 end do 1508 1509 !Compute the number of k points in the G-space unit cell 1510 nkptlatt=kptrlatt2(1,1)*kptrlatt2(2,2)*kptrlatt2(3,3) & 1511 & +kptrlatt2(1,2)*kptrlatt2(2,3)*kptrlatt2(3,1) & 1512 & +kptrlatt2(1,3)*kptrlatt2(2,1)*kptrlatt2(3,2) & 1513 & -kptrlatt2(1,2)*kptrlatt2(2,1)*kptrlatt2(3,3) & 1514 & -kptrlatt2(1,3)*kptrlatt2(2,2)*kptrlatt2(3,1) & 1515 & -kptrlatt2(1,1)*kptrlatt2(2,3)*kptrlatt2(3,2) 1516 1517 !Check whether the number of k points is positive, otherwise, change the handedness of kptrlatt2 1518 if(nkptlatt<=0)then 1519 ! write(std_out,*)' getkgrid : nkptlatt is negative !' 1520 kptrlatt2(:,3)=-kptrlatt2(:,3) 1521 nkptlatt=-nkptlatt 1522 do ishiftk=1,nshiftk2 1523 shiftk2(3,ishiftk)=-shiftk2(3,ishiftk) 1524 end do 1525 end if 1526 1527 !Determine the smallest supercell R-vector whose contribution 1528 !is not taken correctly into account in the k point integration. 1529 !Increase enormously the size of the cell when vacuum is present. 1530 fact_vacuum(:)=1 1531 if(vacuum(1)==1)fact_vacuum(1)=1000.0_dp 1532 if(vacuum(2)==1)fact_vacuum(2)=1000.0_dp 1533 if(vacuum(3)==1)fact_vacuum(3)=1000.0_dp 1534 do ii=1,3 1535 rprimd_super(:,ii)=fact_vacuum(1)*rprimd(:,1)*kptrlatt2(1,ii)+& 1536 & fact_vacuum(2)*rprimd(:,2)*kptrlatt2(2,ii)+& 1537 & fact_vacuum(3)*rprimd(:,3)*kptrlatt2(3,ii) 1538 end do 1539 1540 call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super) 1541 call smallprim(metmin,minim,rprimd_super) 1542 length2=min(metmin(1,1),metmin(2,2),metmin(3,3)) 1543 kptrlen=sqrt(length2) 1544 1545 !write(msg,'(a,es16.6)' )' getkgrid : length of smallest supercell vector (bohr)=',kptrlen 1546 !call wrtout(std_out,msg,'COLL') 1547 ! If the number of shifts has been decreased, determine the set of kptrlatt2 vectors 1548 ! with minimal length (without using fact_vacuum) 1549 ! It is worth to determine the minimal set of vectors so that the kptrlatt that is output 1550 ! does not seem screwy, although correct but surprising. 1551 if(nshiftk/=nshiftk2)then 1552 do ii=1,3 1553 rprimd_super(:,ii)=rprimd(:,1)*kptrlatt2(1,ii)+rprimd(:,2)*kptrlatt2(2,ii)+rprimd(:,3)*kptrlatt2(3,ii) 1554 end do 1555 call metric(gmet_super,gprimd_super,-1,rmet_super,rprimd_super,ucvol_super) 1556 ! Shift vectors in cartesian coordinates (reciprocal space) 1557 do ishiftk=1,nshiftk2 1558 shiftk3(:,ishiftk)=gprimd_super(:,1)*shiftk2(1,ishiftk)+& 1559 & gprimd_super(:,2)*shiftk2(2,ishiftk)+& 1560 & gprimd_super(:,3)*shiftk2(3,ishiftk) 1561 end do 1562 call smallprim(metmin,minim,rprimd_super) 1563 call metric(gmet_super,gprimd_super,-1,rmet_super,minim,ucvol_super) 1564 ! This is the new kptrlatt2 1565 do ii=1,3 1566 dijk(:)=gprimd(1,:)*minim(1,ii)+& 1567 & gprimd(2,:)*minim(2,ii)+& 1568 & gprimd(3,:)*minim(3,ii) 1569 kptrlatt2(:,ii)=nint(dijk(:)) 1570 end do 1571 ! Shifts in the new set of kptrlatt vectors 1572 do ishiftk=1,nshiftk2 1573 shiftk2(:,ishiftk)=minim(1,:)*shiftk3(1,ishiftk)+& 1574 & minim(2,:)*shiftk3(2,ishiftk)+& 1575 & minim(3,:)*shiftk3(3,ishiftk) 1576 end do 1577 end if 1578 1579 !brav=1 is able to treat all bravais lattices. 1580 brav=1 1581 mkpt=nkptlatt*nshiftk2 1582 1583 ABI_MALLOC(spkpt,(3,mkpt)) 1584 option=0 1585 if(iout/=0)option=1 1586 1587 !call cwtime_report(' shifts', cpu, wall, gflops) 1588 1589 if (present(downsampling))then 1590 call smpbz(brav,iout,kptrlatt2,mkpt,nkpthf_computed,nshiftk2,option,shiftk2,spkpt,downsampling=downsampling) 1591 if (present(kpthf) .and. nkpthf/=0) then 1592 ! Returns list of k-points in the Full BZ, possibly downsampled for Fock 1593 kpthf = spkpt(:,1:nkpthf) 1594 end if 1595 nkpthf=nkpthf_computed 1596 end if 1597 1598 call smpbz(brav,iout,kptrlatt2,mkpt,nkpt_fullbz,nshiftk2,option,shiftk2,spkpt) 1599 !call cwtime_report(' smpbz', cpu, wall, gflops) 1600 1601 if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then 1602 1603 ABI_MALLOC(indkpt,(nkpt_fullbz)) 1604 ABI_MALLOC(kpt_fullbz,(3,nkpt_fullbz)) 1605 ABI_MALLOC(bz2ibz_smap, (6, nkpt_fullbz)) 1606 #if 1 1607 ABI_MALLOC(wtk_fullbz,(nkpt_fullbz)) 1608 ABI_MALLOC(wtk_folded,(nkpt_fullbz)) 1609 1610 kpt_fullbz(:,:)=spkpt(:,1:nkpt_fullbz) 1611 wtk_fullbz(1:nkpt_fullbz)=1.0_dp/dble(nkpt_fullbz) 1612 1613 timrev=1;if (kptopt==4) timrev=0 1614 1615 indkpt = 0 1616 call symkpt(chksymbreak,gmet,indkpt,iout,kpt_fullbz,nkpt_fullbz,& 1617 & nkpt_computed,nsym_used,symrec,timrev,wtk_fullbz,wtk_folded,bz2ibz_smap,xmpi_comm_self) 1618 1619 ABI_FREE(symrec) 1620 ABI_FREE(wtk_fullbz) 1621 1622 !do ikpt=1,nkpt_fullbz 1623 ! write(*,*) ikpt, indkpt(ikpt), bz2ibz_smap(1,ikpt), indkpt(bz2ibz_smap(1,ikpt)) 1624 !end do 1625 #else 1626 kpt_fullbz(:,:)=spkpt(:,1:nkpt_fullbz) 1627 1628 timrev=1;if (kptopt==4) timrev=0 1629 1630 call symkpt_new(chksymbreak,gmet,indkpt,iout,kpt_fullbz,nkpt_fullbz,& 1631 & nkpt_computed,nsym_used,symrec,timrev,bz2ibz_smap,xmpi_comm_self) 1632 1633 ABI_FREE(symrec) 1634 ABI_CALLOC(wtk_folded,(nkpt_fullbz)) 1635 do ii=1,nkpt_fullbz 1636 ikpt = indkpt(bz2ibz_smap(1,ii)) 1637 wtk_folded(ikpt) = wtk_folded(ikpt) + one 1638 end do 1639 wtk_folded = wtk_folded / nkpt_fullbz 1640 #endif 1641 1642 else if(kptopt==3)then 1643 ABI_ICALLOC(bz2ibz_smap, (6, nkpt_fullbz)) 1644 bz2ibz_smap(1,:) = [(ii,ii=1,nkpt_fullbz)] 1645 bz2ibz_smap(2,:) = 1 !isym 1646 nkpt_computed=nkpt_fullbz 1647 end if 1648 !call cwtime_report(' symkpt', cpu, wall, gflops) 1649 1650 !The number of k points has been computed from kptopt, kptrlatt, nshiftk, shiftk, 1651 !and the eventual symmetries, it is presently called nkpt_computed. 1652 nkpt_use = nkpt 1653 if (nkpt<0) nkpt_use = nkpt_computed 1654 1655 !Check that the argument nkpt is coherent with nkpt_computed, if nkpt/=0. 1656 if(nkpt_use/=nkpt_computed .and. nkpt/=0)then 1657 write(msg, '(a,i0,5a,i0,7a)') & 1658 'The argument nkpt = ',nkpt_use,', does not match',ch10,& 1659 'the number of k points generated by kptopt, kptrlatt, shiftk,',ch10,& 1660 'and the eventual symmetries, that is, nkpt= ',nkpt_computed,'.',ch10,& 1661 'However, note that it might be due to the user,',ch10,& 1662 'if nkpt is explicitely defined in the input file.',ch10,& 1663 'In this case, please check your input file.' 1664 ABI_BUG(msg) 1665 end if 1666 1667 ABI_MALLOC(kpt,(3,nkpt_use)) 1668 ABI_MALLOC(wtk,(nkpt_use)) 1669 1670 if(kptopt==1 .or. kptopt==2 .or. kptopt==4)then 1671 1672 if(nkpt_use/=0)then 1673 do ikpt=1,nkpt_use 1674 kpt(:,ikpt)=kpt_fullbz(:,indkpt(ikpt)) 1675 if(iscf>=0 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(ikpt)=wtk_folded(indkpt(ikpt)) 1676 end do 1677 end if 1678 1679 if (present(fullbz)) then 1680 ! Returns list of k-points in the Full BZ. 1681 ABI_MOVE_ALLOC(kpt_fullbz,fullbz) 1682 else 1683 ABI_FREE(kpt_fullbz) 1684 end if 1685 1686 ABI_FREE(wtk_folded) 1687 1688 else if(kptopt==3)then 1689 1690 if(nkpt_use/=0)then 1691 kpt(:,1:nkpt_use)=spkpt(:,1:nkpt_use) 1692 if(iscf>1 .or. iscf==-3 .or. iscf==-1.or.iscf==-2)wtk(1:nkpt_use)=1.0_dp/dble(nkpt_use) 1693 end if 1694 1695 if (present(fullbz)) then 1696 ! Returns list of k-points in the Full BZ. 1697 ABI_MALLOC(fullbz,(3,nkpt_fullbz)) 1698 fullbz = spkpt(:,1:nkpt_fullbz) 1699 end if 1700 1701 end if 1702 1703 ABI_FREE(spkpt) 1704 kptrlatt(:,:)=kptrlatt2(:,:) 1705 nshiftk=nshiftk2 1706 shiftk(:,1:nshiftk)=shiftk2(:,1:nshiftk) 1707 ABI_FREE(shiftk2) 1708 ABI_FREE(shiftk3) 1709 1710 end subroutine getkgrid_low
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 spatial 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 [bz2ibz(6,nkbz)]=Mapping BZ --> IBZ
SOURCE
130 subroutine kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, nkibz, kibz, wtk, nkbz, kbz, & 131 new_kptrlatt, new_shiftk, bz2ibz) ! Optional 132 133 !Arguments ------------------------------------ 134 !scalars 135 integer,intent(in) :: nshiftk,kptopt 136 integer,intent(out) :: nkibz,nkbz 137 type(crystal_t),intent(in) :: cryst 138 !arrays 139 integer,intent(in) :: kptrlatt(3,3) 140 integer,optional,allocatable,intent(out) :: bz2ibz(:,:) 141 integer,optional,intent(out) :: new_kptrlatt(3,3) 142 real(dp),intent(in) :: shiftk(3,nshiftk) 143 real(dp),allocatable,intent(out) :: wtk(:),kibz(:,:),kbz(:,:) 144 real(dp),optional,allocatable,intent(out) :: new_shiftk(:,:) 145 146 !Local variables------------------------------- 147 !scalars 148 integer,parameter :: iout0 = 0, chksymbreak0 = 0, iscf2 = 2 149 integer :: my_nshiftk 150 real(dp) :: kptrlen 151 !arrays 152 integer,parameter :: vacuum0(3) = [0, 0, 0] 153 integer :: my_kptrlatt(3,3) 154 integer,allocatable :: indkpt(:),bz2ibz_smap(:,:) 155 real(dp) :: my_shiftk(3,MAX_NSHIFTK) 156 157 ! ********************************************************************* 158 159 ! Copy kptrlatt and shifts because getkgrid can change them 160 ! Be careful as getkgrid expects shiftk(3,MAX_NSHIFTK). 161 ABI_CHECK_IRANGE(nshiftk, 1, MAX_NSHIFTK, "Invalid value of nshiftk") 162 my_nshiftk = nshiftk; my_shiftk = zero; my_shiftk(:,1:nshiftk) = shiftk 163 my_kptrlatt = kptrlatt 164 165 call getkgrid_low(chksymbreak0,iout0,iscf2,kibz,kptopt,my_kptrlatt,kptrlen,& 166 cryst%nsym,-1,nkibz,my_nshiftk,cryst%nsym,cryst%rprimd,my_shiftk,cryst%symafm,& 167 cryst%symrel,vacuum0,wtk,indkpt,bz2ibz_smap,fullbz=kbz) 168 169 if (present(bz2ibz)) then 170 ABI_MOVE_ALLOC(bz2ibz_smap, bz2ibz) 171 else 172 ABI_SFREE(bz2ibz_smap) 173 endif 174 ABI_SFREE(indkpt) 175 176 nkbz = size(kbz, dim=2) 177 178 ! Optionally, return new shifts and new_kptrlatt 179 if (present(new_shiftk)) then 180 ABI_MALLOC(new_shiftk, (3, my_nshiftk)) 181 new_shiftk = my_shiftk(:, 1:my_nshiftk) 182 end if 183 if (present(new_kptrlatt)) new_kptrlatt = my_kptrlatt 184 185 DBG_CHECK(abs(sum(wtk) - one) < tol10, "sum(wtk) != one") 186 187 end subroutine kpts_ibz_from_kptrlatt
m_kpts/kpts_map [ Functions ]
[ Top ] [ m_kpts ] [ Functions ]
NAME
kpts_map
FUNCTION
Compute symmetry table.
INPUTS
OUTPUT
SOURCE
601 integer function kpts_map(mode, timrev, cryst, krank, nkpt2, kpt2, map, qpt, dksqmax_tol) result(ierr) 602 603 !Arguments ------------------------------------ 604 !scalars 605 character(len=*),intent(in) :: mode 606 integer,intent(in) :: timrev, nkpt2 607 class(crystal_t),intent(in) :: cryst 608 class(krank_t),intent(inout) :: krank 609 real(dp),optional,intent(in) :: dksqmax_tol 610 !arrays 611 real(dp),intent(in) :: kpt2(3, nkpt2) 612 real(dp),optional,intent(in) :: qpt(3) 613 integer,intent(out) :: map(6, nkpt2) 614 615 !Local variables------------------------------- 616 !scalars 617 real(dp) :: dksqmax, my_tol 618 !arrays 619 real(dp) :: my_qpt(3) 620 621 ! ************************************************************************* 622 623 my_qpt = zero; if (present(qpt)) my_qpt = qpt 624 625 select case (mode) 626 case("symrel") 627 ! Note symrel and use_symrec = .False. 628 ! These are the conventions for the symmetrization of the wavefunctions used in cgtk_rotate. 629 call krank%get_mapping(nkpt2, kpt2, dksqmax, cryst%gmet, map, & 630 cryst%nsym, cryst%symafm, cryst%symrel, timrev, use_symrec=.False., qpt=my_qpt) 631 632 case("symrec") 633 ! Note symrec and use_symrec = .True. 634 ! These are the conventions for the symmetrization of the DVDB as well as the conventions 635 ! used in several BZ routines. We should always use this convention. 636 637 call krank%get_mapping(nkpt2, kpt2, dksqmax, cryst%gmet, map, & 638 cryst%nsym, cryst%symafm, cryst%symrec, timrev, use_symrec=.True., qpt=my_qpt) 639 640 case default 641 ABI_ERROR(sjoin("Invalid mode:", mode)) 642 end select 643 644 my_tol = tol12; if (present(dksqmax_tol)) my_tol = dksqmax_tol 645 646 ierr = merge(1, 0, dksqmax > my_tol) 647 if (ierr /= 0) call wrtout(std_out, sjoin(" CRITICAL WARNING: dksqmax ", ftoa(dksqmax), " > ", ftoa(my_tol))) 648 649 end function kpts_map
m_kpts/kpts_map_print [ Functions ]
[ Top ] [ m_kpts ] [ Functions ]
NAME
kpts_map_print
FUNCTION
Print the symmetry table bz2ibz associated to the BZ bz and the IBZ ibz to a list of units with header. Mode corresponds to the value passed to kpts_map If prtvol is 0, max 20 entries are printed. Use prtvol > 0 to print all k-points.
SOURCE
663 subroutine kpts_map_print(units, header, mode, bz, ibz, bz2ibz, prtvol) 664 665 !Arguments ------------------------------------ 666 !scalars 667 character(len=*),intent(in) :: header, mode 668 integer,intent(in) :: prtvol, units(:), bz2ibz(:,:) 669 real(dp),intent(in) :: bz(:,:), ibz(:,:) 670 671 !Local variables------------------------------- 672 !scalars 673 integer :: ik_ibz, ik_bz, isym_k, trev_k, g0_k(3) 674 logical :: isirr_k 675 character(len=5000) :: msg 676 677 ! ************************************************************************* 678 679 call wrtout(units, " "//trim(header)) 680 681 select case (mode) 682 case ("symrec") 683 call wrtout(units, & 684 " Legend: bz = TS(ibz) + g0 where isym is the index of the symrec operation S and itim is 1 if TR is used.") 685 case ("symrel") 686 call wrtout(units, & 687 " Legend: bz = TS^t(ibz) + g0 where isym is the index of the symrel operation S and itim is 1 if TR is used.") 688 case default 689 ABI_ERROR(sjoin("Invalid mode:", mode)) 690 end select 691 692 ! yes, I'm a barbarian but Fortran string formatting is a pain. 693 msg = " BZ IBZ ibz isym itim g0" 694 call wrtout(units, msg) 695 696 do ik_bz=1,size(bz2ibz, dim=2) 697 if (prtvol == 0 .and. ik_bz > 20) then 698 call wrtout(units, "prtvol = 0, max 20 points are written"); exit 699 end if 700 ik_ibz = bz2ibz(1, ik_bz); isym_k = bz2ibz(2, ik_bz) 701 trev_k = bz2ibz(6, ik_bz); g0_k = bz2ibz(3:5, ik_bz) 702 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 703 write(msg, '(i6, 2x, 2(a,2x), 3(i4,2x), a)' ) & 704 ik_bz, trim(ktoa(bz(:, ik_bz))), trim(ktoa(ibz(:,ik_ibz))), ik_ibz, isym_k, trev_k, trim(ltoa(g0_k)) 705 call wrtout(units, msg) 706 end do 707 708 call wrtout(units, ch10) 709 710 end subroutine kpts_map_print
m_kpts/kpts_pack_in_stars [ Functions ]
[ Top ] [ m_kpts ] [ Functions ]
NAME
kpts_pack_in_stars
FUNCTION
INPUTS
OUTPUT
SOURCE
506 subroutine kpts_pack_in_stars(nkpt, kpts, kmap) 507 508 !Arguments ------------------------------------ 509 !scalars 510 integer,intent(in) :: nkpt 511 !arrays 512 real(dp),intent(inout) :: kpts(3, nkpt) 513 integer,intent(inout) :: kmap(6, nkpt) 514 515 !Local variables------------------------------- 516 !scalars 517 integer :: ikpt, seen_ibz, ik_start, ik0, nkibz 518 integer :: ik_ibz, isym_k, trev_k, tsign, g0_k(3) 519 logical :: isirr_k 520 !arrays 521 integer,allocatable :: iperm(:), ibz_ids(:), kmap_ord(:,:), star_pos(:,:) 522 real(dp) :: swap_kpt(3), swap_kmap(6) 523 real(dp),allocatable :: kpts_ord(:,:) 524 525 ! ************************************************************************* 526 527 ! Order according to ik_ibz index 528 ABI_MALLOC(ibz_ids, (nkpt)) 529 ABI_MALLOC(iperm, (nkpt)) 530 ibz_ids = kmap(1, :) 531 iperm = [(ikpt, ikpt=1, nkpt)] 532 533 call sort_int(nkpt, ibz_ids, iperm) 534 ABI_FREE(ibz_ids) 535 536 ! Rearrange items in _ord arrays. 537 ABI_MALLOC(kpts_ord, (3, nkpt)) 538 ABI_MALLOC(kmap_ord, (6, nkpt)) 539 do ikpt=1,nkpt 540 kpts_ord(:, ikpt) = kpts(:, iperm(ikpt)) 541 kmap_ord(:, ikpt) = kmap(:, iperm(ikpt)) 542 end do 543 544 ! We want each star group to start with the point in the IBZ so an extra shuffle is needed. 545 ! star_pos stores the beginning of the star group and the position of the base k0 for each star. 546 nkibz = maxval(kmap(1,:)) 547 ABI_ICALLOC(star_pos, (2, nkibz)) 548 549 seen_ibz = -1 550 do ikpt=1,nkpt 551 ik_ibz = kmap_ord(1, ikpt); isym_k = kmap_ord(2, ikpt) 552 trev_k = kmap_ord(6, ikpt); g0_k = kmap_ord(3:5, ikpt) 553 isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0)) 554 tsign = 1; if (trev_k == 1) tsign = -1 555 if (ik_ibz /= seen_ibz) then 556 star_pos(1, ik_ibz) = ikpt 557 seen_ibz = ik_ibz 558 end if 559 if (isirr_k) star_pos(2, ik_ibz) = ikpt 560 end do 561 562 ! Now put ik0 in the ik_start slot if needed. 563 do ik_ibz=1, nkibz 564 ik_start = star_pos(1, ik_ibz) 565 ik0 = star_pos(2, ik_ibz) 566 if (ik_start == ik0) cycle 567 568 swap_kpt = kpts_ord(:, ik_start) 569 swap_kmap = kmap_ord(:, ik_start) 570 571 kpts_ord(:, ik_start) = kpts_ord(:, ik0) 572 kmap_ord(:, ik_start) = kmap_ord(:, ik0) 573 kpts_ord(:, ik0) = swap_kpt 574 kmap_ord(:, ik0) = swap_kmap 575 end do 576 577 kpts = kpts_ord 578 kmap = kmap_ord 579 580 ABI_FREE(star_pos) 581 ABI_FREE(iperm) 582 ABI_FREE(kpts_ord) 583 ABI_FREE(kmap_ord) 584 585 end subroutine kpts_pack_in_stars
m_kpts/kpts_sort [ Functions ]
[ Top ] [ m_kpts ] [ Functions ]
NAME
kpts_sort
FUNCTION
Order list of k-points according to the norm.
INPUTS
OUTPUT
SOURCE
454 subroutine kpts_sort(gprimd, nkpt, kpts) 455 456 !Arguments ------------------------------------ 457 !scalars 458 integer,intent(in) :: nkpt 459 !arrays 460 real(dp),intent(in) :: gprimd(3, 3) 461 real(dp),intent(inout) :: kpts(3, nkpt) 462 463 !Local variables------------------------------- 464 !scalars 465 integer :: ikpt 466 !arrays 467 integer,allocatable :: iperm(:) 468 real(dp),allocatable :: knorm2(:), kpts_ord(:,:) 469 470 ! ************************************************************************* 471 472 ABI_MALLOC(knorm2, (nkpt)) 473 do ikpt=1,nkpt 474 knorm2(ikpt) = dot_product(kpts(:,ikpt), matmul(gprimd, kpts(:, ikpt))) 475 end do 476 477 ABI_MALLOC(iperm, (nkpt)) 478 iperm = [(ikpt, ikpt=1, nkpt)] 479 call sort_dp(nkpt, knorm2, iperm, tol12) 480 ABI_FREE(knorm2) 481 482 ABI_MALLOC(kpts_ord, (3, nkpt)) 483 do ikpt=1,nkpt 484 kpts_ord(:, ikpt) = kpts(:, iperm(ikpt)) 485 end do 486 kpts = kpts_ord 487 488 ABI_FREE(iperm) 489 ABI_FREE(kpts_ord) 490 491 end subroutine kpts_sort
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 spatial symmetries and/or time-reversal can be used)
SOURCE
87 integer pure function kpts_timrev_from_kptopt(kptopt) result(timrev) 88 89 !Arguments ------------------------------------ 90 !scalars 91 integer,intent(in) :: kptopt 92 93 ! ********************************************************************* 94 95 timrev = 1; if (any(kptopt == [3, 4])) timrev = 0 96 97 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 kpt 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 comm=MPI communicator. [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.
SOURCE
765 subroutine listkk(dksqmax, gmet, indkk, kptns1, kptns2, nkpt1, nkpt2, nsym, sppoldbl, symafm, symmat, timrev, comm, & 766 use_symrec) ! optional 767 768 !Arguments ------------------------------------ 769 !scalars 770 integer,intent(in) :: nkpt1,nkpt2,nsym,sppoldbl,timrev,comm 771 real(dp),intent(out) :: dksqmax 772 logical,optional,intent(in) :: use_symrec 773 !arrays 774 integer,intent(in) :: symafm(nsym),symmat(3,3,nsym) 775 integer,intent(out) :: indkk(nkpt2*sppoldbl,6) 776 real(dp),intent(in) :: gmet(3,3),kptns1(3,nkpt1),kptns2(3,nkpt2) 777 778 !Local variables------------------------------- 779 !scalars 780 integer,parameter :: usesym=1, limit=1 781 integer :: nprocs, my_rank, ierr, isk_start, isk_stop 782 integer :: l3,ig1,ig2,ig3,ii,ikpg1,ikpt1,ikpt2,ikpt2_done, isk 783 integer :: ilarger,ismaller,itrial 784 integer :: isppol,isym,itimrev,jkpt1,jsym,jtime 785 integer :: nsym_used,timrev_used 786 real(dp) :: dksq,dksqmn,lk2,llarger,ldiff,lsmaller,ltrial,min_l 787 !real(dp) :: cpu,wall,gflops 788 character(len=500) :: msg 789 !arrays 790 integer :: dkint(3),jdkint(3),k1int(3),k2int(3) 791 integer, allocatable :: isort(:), tmp_indkk(:,:) 792 real(dp) :: tsec(2) 793 real(dp) :: dk(3),kpg1(3),kpt1a(3),k1(3),k2(3) 794 !real(dp) :: kasq,ka(3) 795 real(dp),allocatable :: lkpg1(:),lkpg1_sorted(:) 796 797 ! ************************************************************************* 798 799 call timab(1091, 1, tsec) 800 !call cwtime(cpu, wall, gflops, "start") 801 802 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 803 804 if (sppoldbl<1 .or. sppoldbl>2) then 805 write(msg, '(a,i0,a)' )'The value of sppoldbl is: ',sppoldbl,', but it should be either 1 or 2.' 806 ABI_BUG(msg) 807 end if 808 809 ! When usesym=0, the old way of converting the wavefunctions (without using the symmetries), is recovered. 810 nsym_used=nsym 811 timrev_used=timrev 812 if(usesym==0)nsym_used=1 813 if(usesym==0)timrev_used=0 814 815 ! Precompute the length of the kpt1 vectors, also taking into account possible umklapp vectors 816 l3 = (2*limit+1)**3 817 ABI_CALLOC(lkpg1, (l3*nkpt1)) 818 ABI_CALLOC(lkpg1_sorted, (l3*nkpt1)) 819 ABI_MALLOC(isort, (l3*nkpt1)) 820 isort = 0 821 822 call xmpi_split_work(nkpt1, comm, isk_start, isk_stop) 823 !write(std_out,*)' List of kpt1 vectors'; write(std_out,*)' Length of the kpt1 vectors:' 824 825 !$OMP PARALLEL DO PRIVATE(k1, k1int, kpg1, ikpg1) 826 do ikpt1=isk_start,isk_stop 827 k1(:) = kptns1(:,ikpt1) !; write(std_out,*)ikpt1,k1(:) 828 k1int(:) = nint(k1(:) + tol12) 829 k1(:) = k1(:) - k1int(:) 830 do ig3=-limit,limit 831 kpg1(3) = k1(3) + ig3 832 do ig2=-limit,limit 833 kpg1(2) = k1(2) + ig2 834 do ig1=-limit,limit 835 kpg1(1) = k1(1) + ig1 836 837 ikpg1 = ig1 + limit + 1 + (2*limit+1)*(ig2+limit) + (2*limit+1)**2*(ig3+limit) + l3*(ikpt1-1) 838 ! Compute the norm of the vector (also taking into account possible umklapp) 839 lkpg1(ikpg1) = sqrt(gmet(1,1)*kpg1(1)**2+gmet(2,2)*kpg1(2)**2 + & 840 gmet(3,3)*kpg1(3)**2+two*(gmet(2,1)*kpg1(2)*kpg1(1) + & 841 gmet(3,2)*kpg1(3)*kpg1(2)+gmet(3,1)*kpg1(3)*kpg1(1))) 842 lkpg1_sorted(ikpg1) = lkpg1(ikpg1) 843 isort(ikpg1) = ikpg1 844 !write(std_out,*)' ikpt1,ig1,ig2,ig3,lkpg1=',ikpt1,ig1,ig2,ig3,lkpg1(ikpg1) 845 end do 846 end do 847 end do 848 end do 849 850 if (nprocs > 1) then 851 call xmpi_sum(lkpg1_sorted, comm, ierr) 852 call xmpi_sum(lkpg1, comm, ierr) 853 call xmpi_sum(isort, comm, ierr) 854 end if 855 !call cwtime_report(" listkk_loop1", cpu, wall, gflops) 856 857 call sort_dp(l3*nkpt1, lkpg1_sorted, isort, tol12) 858 ! From "precompute" to "sort_dp" represents more than 50% of the overall wall time for large meshes. 859 !call cwtime_report(" listkk_sort", cpu, wall, gflops) 860 861 !write(std_out,*)' listkk : output list of kpt1 for checking purposes ' 862 !write(std_out,*)' ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii)) ' 863 !do ii=1,l3*nkpt1 864 ! ikpt1=(isort(ii)-1)/l3+1 865 ! write(std_out,*)ii,ikpt1,isort(ii)-l3*(ikpt1-1),lkpg1_sorted(ii),lkpg1(isort(ii)) 866 !enddo 867 868 dksqmax = zero 869 indkk = 0 870 ! TODO: Should change API to use this shape. 871 ! workspace array for improved memory access. 872 ABI_MALLOC(tmp_indkk, (6, nkpt2*sppoldbl)) 873 tmp_indkk = 0 874 875 ! Split loop in contiguous blocks 876 call xmpi_split_work(sppoldbl * nkpt2, comm, isk_start, isk_stop) 877 878 do isppol=1,sppoldbl 879 do ikpt2=1,nkpt2 880 isk = ikpt2 + (isppol-1)*nkpt2 881 if (isk < isk_start .or. isk > isk_stop) cycle 882 883 ikpt2_done=0 884 ! Precompute the length of the kpt2 vector, with the Umklapp vector such that it is the closest to the Gamma point 885 k2(:)=kptns2(:,ikpt2) 886 k2int(:)=nint(k2(:)+tol12) 887 k2(:)=k2(:)-k2int(:) 888 lk2=sqrt(gmet(1,1)*k2(1)**2+gmet(2,2)*k2(2)**2+& 889 gmet(3,3)*k2(3)**2+two*(gmet(2,1)*k2(2)*k2(1)+& 890 gmet(3,2)*k2(3)*k2(2)+gmet(3,1)*k2(3)*k2(1))) 891 ! write(std_out, '(a,i4,7es16.6)' )' listkk : ikpt2,kptns2(:,ikpt2),k2(:),lk2=',ikpt2,kptns2(:,ikpt2),k2(:),lk2 892 893 ! Find the kpt1 vector whose length is the most similar to the length of lk2 up to a tolerance. 894 ! Use a bisection algorithm. 895 ismaller=0; lsmaller=zero 896 ilarger=l3*nkpt1+1; llarger=huge(one) 897 898 ! This loop should never reach l3*nkpt1, since this is a bisection algorithm 899 do ii=1,l3*nkpt1 900 if((ilarger-ismaller)<2 .or. (llarger-lsmaller)<2*tol12)exit 901 itrial=(ilarger+ismaller)/2 ; ltrial=lkpg1_sorted(itrial) 902 if((ltrial-lk2)>tol12)then 903 ilarger=itrial ; llarger=ltrial 904 else if((ltrial-lk2)<-tol12)then 905 ismaller=itrial ; lsmaller=ltrial 906 else 907 ismaller=itrial ; lsmaller=ltrial 908 ilarger=itrial ; llarger=ltrial 909 end if 910 end do 911 itrial=ismaller 912 if(abs(llarger-lk2)<abs(lsmaller-lk2)-tol12)itrial=ilarger 913 if(itrial==0)itrial=ilarger 914 ismaller=itrial ; ilarger=itrial 915 !write(std_out,*)' listkk : starting search at itrial=',itrial 916 917 dksqmn=huge(one) 918 919 ! The ii index is dummy. This avoids an infinite loop. 920 do ii=1,l3*nkpt1 921 ! If the difference in length between the trial vector and the target vector is bigger 922 ! than the already achieved distance, the search is finished ... 923 ldiff = abs(lkpg1_sorted(itrial) - lk2) 924 ! write(std_out,*)' listkk : ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,& 925 ! dksqmn=',ii,itrial,lkpg1_sorted(itrial),lk2,ldiff,dksqmn 926 927 if (ldiff**2 > dksqmn+tol8) exit 928 929 ! If this k-point has already been examined in a previous batch, skip it 930 ! First, compute the minimum of the difference of length of the sets of 931 ! associated vectors thanks to Umklapp vectors with the target vector 932 ikpt1 = (isort(itrial)-1) /l3 + 1 933 min_l = minval(abs(lkpg1((ikpt1-1)*l3+1:(ikpt1-1)*l3+l3)-lk2)) 934 935 ! Then compare with the current ldiff 936 ! write(std_out,*)' listkk : ikpt1,min_l,ldiff=',ikpt1,min_l,ldiff 937 if (min_l > ldiff-tol12) then 938 939 ! Now, will examine the trial vector, and the symmetric ones 940 ! MG FIXME: Here there's a possible problem with the order of symmetries because 941 ! in symkpt, time-reversal is the innermost loop. This can create inconsistencies in the symmetry tables. 942 ! Besides, one should use symrel^{-1 T} to keep the correspondence between isym -> R or S 943 do itimrev=0,timrev_used 944 do isym=1,nsym_used 945 946 ! Select magnetic characteristic of symmetries 947 if (isppol == 1 .and. symafm(isym) == -1) cycle 948 if (isppol == 2 .and. symafm(isym) == 1) cycle 949 950 ! Compute symmetric point to kpt1 951 if (usesym==1) then 952 ! original code only used transpose(symrel) 953 if (present(use_symrec)) then 954 if (use_symrec) then 955 kpt1a(:) = MATMUL(symmat(:,:,isym),kptns1(:,ikpt1)) 956 else 957 kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1)) 958 end if 959 else 960 kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,isym)),kptns1(:,ikpt1)) 961 end if 962 kpt1a(:)=(1-2*itimrev)*kpt1a(:) 963 else 964 kpt1a(:)=kptns1(:,ikpt1) 965 end if 966 967 ! Compute difference with respect to kpt2, modulo a lattice vector 968 dk(:)=kptns2(:,ikpt2)-kpt1a(:) 969 if (usesym==1) then 970 ! The tolerance insure similar behaviour on different platforms 971 ! XG120418: Actually, *assumes* that the closest point will have reduced 972 ! coordinates differing by less than 1/2. There might be elongated cells where this is not correct ... 973 dkint(:)=nint(dk(:)+tol12) 974 dk(:)=dk(:)-dkint(:) 975 else 976 dkint(:)=0 977 end if 978 979 ! Compute norm of the difference vector, and update kpt1 if better. 980 dksq=gmet(1,1)*dk(1)**2+gmet(2,2)*dk(2)**2+ & 981 gmet(3,3)*dk(3)**2+two*(gmet(2,1)*dk(2)*dk(1)+ & 982 gmet(3,2)*dk(3)*dk(2)+gmet(3,1)*dk(3)*dk(1)) 983 984 if (dksq < dksqmn+tol8) then 985 ! If exactly the right point (without using symmetries neither umklapp vector), will exit the search 986 ! Note that in this condition, each coordinate is tested separately, without squaring. 987 ! So, it is a much stronger condition than dksqmn < tol12 988 if (sum(abs(kptns2(:,ikpt2)-kptns1(:,ikpt1)))<3*tol12) ikpt2_done = 1 989 990 ! Update in three cases: either if succeeded to have exactly the vector, or the distance is better, 991 ! or the distance is only slightly worsened so select the lowest itimrev, isym or ikpt1, 992 ! in order to respect previous ordering 993 if (ikpt2_done==1 .or. & 994 dksq+tol12<dksqmn .or. & 995 ( abs(dksq-dksqmn)<tol12 .and. & 996 ((itimrev<jtime) .or. & 997 (itimrev==jtime .and. isym<jsym) .or. & 998 (itimrev==jtime .and. isym==jsym .and. ikpt1<jkpt1))))then 999 1000 dksqmn = dksq 1001 jkpt1 = ikpt1 1002 jsym = isym 1003 jtime = itimrev 1004 jdkint(:) = dkint(:) 1005 1006 !if (ikpt2_done == 1) then 1007 ! write(std_out,*)'Succeeded to lower dskmn,ikpt2_done=',dksqmn,ikpt2_done 1008 ! write(std_out,*)' ikpt1,ikpt2=',ikpt1, ikpt2 1009 ! write(std_out,*)' ikpt1,isym,dkint(:),itimrev=',ikpt1,isym,dkint(:),itimrev 1010 ! ka(:) = kpt1a(:) + dkint(:) 1011 ! kasq=gmet(1,1)*ka(1)**2+gmet(2,2)*ka(2)**2+& 1012 ! gmet(3,3)*ka(3)**2+two*(gmet(2,1)*ka(2)*ka(1)+& 1013 ! gmet(3,2)*ka(3)*ka(2)+gmet(3,1)*ka(3)*ka(1)) 1014 ! write(std_out,*)' k1 = ',kpt1a(:) 1015 ! write(std_out,*)' dkint = ',dkint(:) 1016 ! write(std_out,*)' Actual k1 = ',ka(:) 1017 ! write(std_out,*)' k2 = ',kptns2(:,ikpt2) 1018 ! write(std_out,*)' Actual k1sq = ',kasq 1019 !end if 1020 end if 1021 end if 1022 1023 if (ikpt2_done==1) exit 1024 end do ! isym 1025 if (ikpt2_done==1) exit 1026 end do ! itimrev 1027 if (ikpt2_done==1) exit 1028 end if 1029 1030 ! Update the interval that has been explored 1031 if (itrial < ismaller) ismaller = itrial 1032 if (itrial > ilarger) ilarger = itrial 1033 1034 ! Select the next index to be tried (preferably the smaller indices, but this is a bit arbitrary). 1035 ! write(std_out,*)' before choosing the next index :' 1036 ! write(std_out,*)' ismaller,itrial,ilarger=',ismaller,itrial,ilarger 1037 ! write(std_out,*)' lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1)=',& 1038 ! lkpg1_sorted(ismaller-1),lk2,lkpg1_sorted(ilarger+1) 1039 1040 if (ismaller>1 .and. ilarger<l3*nkpt1) then 1041 if (abs(lkpg1_sorted(ismaller-1)-lk2) < abs(lkpg1_sorted(ilarger+1)-lk2)+tol12) then 1042 itrial = ismaller-1 1043 else 1044 itrial = ilarger+1 1045 end if 1046 end if 1047 if (ismaller==1 .and. ilarger<l3*nkpt1) itrial = ilarger+1 1048 if (ismaller>1 .and. ilarger==l3*nkpt1) itrial = ismaller-1 1049 !if(ismaller==1 .and. ilarger==l3*nkpt1), we are done with the loop ! 1050 end do ! ikpt1 1051 1052 ! Store indices (lots of cache miss here) 1053 !indkk(isk, 1) = jkpt1 1054 !indkk(isk, 2) = jsym 1055 !indkk(isk, 3:5) = jdkint(:) 1056 !indkk(isk, 6) = jtime 1057 1058 tmp_indkk(1, isk) = jkpt1 1059 tmp_indkk(2, isk) = jsym 1060 tmp_indkk(3:5, isk) = jdkint(:) 1061 tmp_indkk(6, isk) = jtime 1062 1063 dksqmax = max(dksqmax, dksqmn) 1064 1065 if (dksqmn < -tol12) then 1066 write(msg, '(a,es16.6)' )'The minimum square of dk has negative norm: dksqmn= ',dksqmn 1067 ABI_BUG(msg) 1068 end if 1069 1070 ! DEBUG SECTION 1071 !if (dksqmn > tol5) then 1072 ! if (present(use_symrec)) then 1073 ! if (use_symrec) then 1074 ! kpt1a(:) = MATMUL(symmat(:,:,jsym),kptns1(:,jkpt1)) 1075 ! else 1076 ! kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,jsym)),kptns1(:,jkpt1)) 1077 ! end if 1078 ! else 1079 ! kpt1a(:) = MATMUL(TRANSPOSE(symmat(:,:,jsym)),kptns1(:,jkpt1)) 1080 ! end if 1081 ! kpt1a(:)=(1-2*jtime)*kpt1a(:) 1082 ! print *, "Cannot find k2: ", k2(:) 1083 ! print *, "Rotated TS(k1): ", kpt1a(:) 1084 ! print *, "with k1: ", kptns1(:, jkpt1) 1085 ! print *, "dksqmn: ", dksqmn 1086 !end if 1087 !END DEBUG 1088 1089 !write(std_out,'(a,i6,i2,2x,i6,5i3,es24.14)' )' listkk: ikpt2,isppol,indkk(isk,:)=',ikpt2,isppol,indkk(isk,:),dksqmn 1090 end do ! ikpt2 1091 end do ! isppol 1092 1093 ABI_FREE(isort) 1094 ABI_FREE(lkpg1) 1095 ABI_FREE(lkpg1_sorted) 1096 1097 indkk = transpose(tmp_indkk) 1098 ABI_FREE(tmp_indkk) 1099 if (nprocs > 1) then 1100 call xmpi_sum(indkk, comm, ierr) 1101 dksqmn = dksqmax 1102 call xmpi_max(dksqmn, dksqmax, comm, ierr) 1103 end if 1104 1105 call timab(1091, 2, tsec) 1106 !call cwtime_report(" listkk_end", cpu, wall, gflops) 1107 1108 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.
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.
SOURCE
3209 subroutine mknormpath(nbounds,bounds,gmet,ndiv_small,ndiv,npt_tot,path) 3210 3211 !Arguments ------------------------------------ 3212 !scalars 3213 integer,intent(in) :: nbounds,ndiv_small 3214 integer,intent(inout) :: npt_tot 3215 !arrays 3216 integer,intent(inout) :: ndiv(nbounds-1) 3217 real(dp),intent(in) :: bounds(3,nbounds),gmet(3,3) 3218 real(dp),intent(out),optional :: path(3,npt_tot) 3219 3220 !Local variables------------------------------- 3221 !scalars 3222 integer :: idx,ii,jp 3223 real(dp) :: fct 3224 character(len=500) :: msg 3225 !arrays 3226 real(dp) :: dd(3),lng(nbounds-1) 3227 3228 ! ************************************************************************* 3229 3230 if (ndiv_small<=0) then 3231 write(msg,'(3a,i0)')& 3232 & 'The argument ndiv_small should be a positive number,',ch10,& 3233 & 'however, ndiv_small=',ndiv_small 3234 ABI_ERROR(msg) 3235 end if 3236 3237 do ii=1,nbounds-1 3238 dd(:)=bounds(:,ii+1)-bounds(:,ii) 3239 lng(ii)= sqrt( dd(1)*gmet(1,1)*dd(1)+ & 3240 & dd(2)*gmet(2,2)*dd(2)+ & 3241 & dd(3)*gmet(3,3)*dd(3)+ & 3242 & 2.0d0*(dd(1)*gmet(1,2)*dd(2)+ & 3243 & dd(1)*gmet(1,3)*dd(3)+ & 3244 & dd(2)*gmet(2,3)*dd(3)) & 3245 & ) 3246 end do 3247 write(std_out,*)lng 3248 fct=minval(lng) 3249 3250 !Avoid division by zero if k(:,i+1)=k(:,i) 3251 if (abs(fct)<tol6) then 3252 write(msg,'(3a)')& 3253 & 'found two consecutive points in the path which are equal',ch10,& 3254 & 'This is not allowed, please modify the path in your input file' 3255 ABI_ERROR(msg) 3256 end if 3257 3258 fct=fct/ndiv_small 3259 ndiv(:)=nint(lng(:)/fct) 3260 !The 1 stand for the first point 3261 npt_tot=sum(ndiv)+1 3262 3263 if (.not.present(path)) then 3264 write(msg,'(2a,i8)')ch10,' mknormpath : total number of points on the path: ',npt_tot 3265 call wrtout(std_out,msg,'COLL') 3266 write(msg,'(2a)')ch10,' Number of divisions for each segment of the normalized path: ' 3267 call wrtout(std_out,msg,'COLL') 3268 do ii=1,nbounds-1 3269 write(msg,'(2(3f8.5,a),i5,a)')& 3270 bounds(:,ii),' ==> ',bounds(:,ii+1),' ( ndiv: ',ndiv(ii),' )' 3271 call wrtout(std_out,msg,'COLL') 3272 end do 3273 write(msg,'(a)')ch10 3274 call wrtout(std_out,msg,'COLL') 3275 else 3276 write(msg,'(2a)')ch10,' Normalized Path: ' 3277 call wrtout(std_out,msg,'COLL') 3278 idx=1 3279 do ii=1,nbounds-1 3280 do jp=1,ndiv(ii) 3281 path(:,idx)=bounds(:,ii)+(jp-1)*(path(:,ii+1)-path(:,ii))/ndiv(ii) 3282 write(msg,'(i4,4x,3(f8.5,1x))')idx,path(:,idx) 3283 call wrtout(std_out,msg,'COLL') 3284 idx=idx+1 3285 end do 3286 end do 3287 end if 3288 3289 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]]
SOURCE
1986 subroutine smpbz(brav,iout,kptrlatt,mkpt,nkpt,nshiftk,option,shiftk,spkpt,downsampling) 1987 1988 !Arguments ------------------------------- 1989 !scalars 1990 integer,intent(in) :: brav,iout,mkpt,nshiftk,option 1991 integer,intent(out) :: nkpt 1992 !arrays 1993 integer,intent(in) :: kptrlatt(3,3) 1994 integer,optional,intent(in) :: downsampling(3) 1995 real(dp),intent(in) :: shiftk(3,nshiftk) 1996 real(dp),intent(out) :: spkpt(3,mkpt) 1997 1998 !Local variables ------------------------- 1999 !scalars 2000 integer,parameter :: prtvol=0 2001 integer :: dividedown,ii,ikshft,jj,kk,nkpout,nkptlatt,nn,proddown 2002 real(dp) :: shift 2003 character(len=500) :: msg 2004 !arrays 2005 integer :: ads(3),boundmax(3),boundmin(3),cds(3),coord(3),ngkpt(3) 2006 integer, allocatable :: found1(:,:),found2(:,:),found3(:,:) 2007 real(dp) :: k1(3),k2(3),kcar(3),klatt(3,3),ktest(3),rlatt(3,3) 2008 2009 ! ********************************************************************* 2010 2011 !DEBUG 2012 !write(std_out,*)' smpbz : brav,iout,mkpt,nkpt,option=',brav,iout,mkpt,nkpt,option 2013 !write(std_out,*)' smpbz : kptrlatt(:,:)=',kptrlatt(:,:) 2014 !write(std_out,*)' smpbz : nshiftk=',nshiftk 2015 !write(std_out,*)' smpbz : shiftk(:,:)=',shiftk(:,:) 2016 !write(std_out,*)' smpbz : downsampling(:)=',downsampling(:) 2017 !ENDDEBUG 2018 2019 if(option/=0) call wrtout(iout,' Homogeneous q point set in the B.Z. ','COLL') 2020 2021 if(abs(brav)/=1)then 2022 ! Only generate Monkhorst-Pack lattices 2023 if(kptrlatt(1,2)/=0 .or. kptrlatt(2,1)/=0 .or. & 2024 & kptrlatt(1,3)/=0 .or. kptrlatt(3,1)/=0 .or. & 2025 & kptrlatt(2,3)/=0 .or. kptrlatt(3,2)/=0 ) then 2026 write(msg, '(2a,a,3i0,a,a,3i4,a,a,3i4)' )& 2027 & 'When abs(brav)/=1, kptrlatt must be diagonal, while it is',ch10,& 2028 & 'kptrlatt(:,1)= ',kptrlatt(:,1),ch10,& 2029 & 'kptrlatt(:,2)= ',kptrlatt(:,2),ch10,& 2030 & 'kptrlatt(:,3)= ',kptrlatt(:,3) 2031 ABI_BUG(msg) 2032 end if 2033 2034 ngkpt(1)=kptrlatt(1,1) 2035 ngkpt(2)=kptrlatt(2,2) 2036 ngkpt(3)=kptrlatt(3,3) 2037 ! 2038 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 2039 write(msg, '(5a,i4,a,a,i0,a,a,i0,a,a)' )& 2040 & 'All ngkpt (or ngqpt) must be strictly positive',ch10,& 2041 & 'or all ngk(q)pt must be zero (for Gamma sampling), but :',ch10,& 2042 & 'ngk(q)pt(1) = ',ngkpt(1),ch10,& 2043 & 'ngk(q)pt(2) = ',ngkpt(2),ch10,& 2044 & 'ngk(q)pt(3) = ',ngkpt(3),ch10,& 2045 & 'Action: correct ngkpt or ngqpt in the input file.' 2046 ABI_BUG(msg) 2047 end if 2048 end if 2049 2050 !Just in case the user wants the grid downsampled to the Gamma point, checks that it is present, and possibly exits 2051 if(present(downsampling))then 2052 if(sum(abs(downsampling(:)))==0)then 2053 do ikshft=1,nshiftk 2054 if(sum(abs(shiftk(:,ikshft)))>tol12)cycle 2055 nkpt=1 2056 spkpt(:,1)=zero 2057 return 2058 end do 2059 end if 2060 end if 2061 2062 !********************************************************************* 2063 2064 if(abs(brav)==1)then 2065 2066 ! Compute the number of k points in the G-space unit cell 2067 ! (will be multiplied by nshiftk later). 2068 nkptlatt=kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3) & 2069 & +kptrlatt(1,2)*kptrlatt(2,3)*kptrlatt(3,1) & 2070 & +kptrlatt(1,3)*kptrlatt(2,1)*kptrlatt(3,2) & 2071 & -kptrlatt(1,2)*kptrlatt(2,1)*kptrlatt(3,3) & 2072 & -kptrlatt(1,3)*kptrlatt(2,2)*kptrlatt(3,1) & 2073 & -kptrlatt(1,1)*kptrlatt(2,3)*kptrlatt(3,2) 2074 2075 if(present(downsampling))then 2076 if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then 2077 if(nshiftk>1)then 2078 write(msg, '(a,3i4,2a,i4,4a)' )& 2079 & 'Real downsampling is activated, with downsampling(1:3)=',downsampling(1:3),ch10,& 2080 & 'However, nshiftk must be 1 in this case, while the input nshiftk=',nshiftk,ch10,& 2081 & 'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,& 2082 & 'or set nshiftk=1.' 2083 ABI_ERROR(msg) 2084 end if 2085 proddown=downsampling(1)*downsampling(2)*downsampling(3) 2086 if(proddown/=0)then 2087 dividedown=abs(proddown) 2088 if(minval(downsampling(:))<0)then ! If there is at least one negative number 2089 dividedown=dividedown*2 2090 if(proddown>0)dividedown=dividedown*2 ! If there are two negative numbers 2091 end if 2092 end if 2093 if(mod(nkptlatt,dividedown)==0)then 2094 nkptlatt=nkptlatt/dividedown 2095 else 2096 write(msg, '(a,3i4,2a,i4,4a)' )& 2097 & 'The requested downsampling, with downsampling(1:3)=',downsampling(1:3),ch10,& 2098 & 'is not compatible with kptrlatt=',ch10,& 2099 & kptrlatt(:,:),ch10,& 2100 & 'that gives nkptlatt=',nkptlatt,ch10,& 2101 & 'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,& 2102 & 'or modify your k-point grid and/or your downsampling in order for them to be compatible.' 2103 ABI_ERROR(msg) 2104 end if 2105 end if 2106 end if 2107 2108 ! Simple Lattice 2109 if (prtvol > 0) call wrtout(std_out,' Simple Lattice Grid ','COLL') 2110 if (mkpt<nkptlatt*nshiftk) then 2111 write(msg, '(a,a,a,i8,a,a,a,a,a)' )& 2112 & 'The value of mkpt is not large enough. It should be',ch10,& 2113 & 'at least',nkptlatt*nshiftk,',',ch10,& 2114 & 'Action: set mkpt to that value in the main routine,',ch10,& 2115 & 'and recompile the code.' 2116 ABI_BUG(msg) 2117 end if 2118 2119 ! Build primitive vectors of the k lattice 2120 rlatt(:,:)=kptrlatt(:,:) 2121 call matr3inv(rlatt,klatt) 2122 2123 ! write(std_out,*)' First primitive vector of the k lattice :',klatt(:,1) 2124 ! write(std_out,*)' Second primitive vector of the k lattice :',klatt(:,2) 2125 ! write(std_out,*)' Third primitive vector of the k lattice :',klatt(:,3) 2126 2127 ! Now, klatt contains the three primitive vectors of the k lattice, 2128 ! in reduced coordinates. One builds all k vectors that 2129 ! are contained in the first Brillouin zone, with coordinates 2130 ! in the interval [0,1[ . First generate boundaries of a big box. 2131 2132 do jj=1,3 2133 2134 ! Mathematically, one has to find the coordinates of the corners of a 2135 ! rectangular paralleliped with integer coordinates, that multiplies the klatt primitive cell and allows 2136 ! it to incorporate completely the [0,1]^3 box. Then take the minimum and maximum 2137 ! of these coordinates, and round them negatively and positively to the next integer. 2138 ! This can be done easily using kptrlatt, considering each coordinate in turn 2139 ! and boils down to enlarging the boundaries for jj by the value of kptrlatt(:,jj), 2140 ! acting on boundmin or boundmax depending on the sign ot kptrlatt(:,jj). 2141 ! XG171020 The coding before 171020 was correct, despite being very simple. 2142 boundmin(jj)=0 ; boundmax(jj)=0 2143 do ii=1,3 2144 if(kptrlatt(ii,jj)<0)boundmin(jj)=boundmin(jj)+kptrlatt(ii,jj) 2145 if(kptrlatt(ii,jj)>0)boundmax(jj)=boundmax(jj)+kptrlatt(ii,jj) 2146 end do 2147 2148 ! To accomodate the shifts, boundmin and boundmax don't start from 0, but are enlarged by one 2149 ! positively and/or negatively. 2150 ! XG171020 Coding in v8.6.0 and before was not correct. This one is even simpler actually. 2151 boundmin(jj)=boundmin(jj)-ceiling(maxval(shiftk(jj,:))+tol14) 2152 boundmax(jj)=boundmax(jj)-floor(minval(shiftk(jj,:))-tol14) 2153 2154 end do 2155 2156 if(present(downsampling))then 2157 ABI_MALLOC(found1,(boundmin(2):boundmax(2),boundmin(3):boundmax(3))) 2158 ABI_MALLOC(found2,(boundmin(1):boundmax(1),boundmin(3):boundmax(3))) 2159 ABI_MALLOC(found3,(boundmin(1):boundmax(1),boundmin(2):boundmax(2))) 2160 found1=0 ; found2=0 ; found3=0 2161 end if 2162 2163 nn=1 2164 do kk=boundmin(3),boundmax(3) 2165 coord(3)=kk 2166 do jj=boundmin(2),boundmax(2) 2167 coord(2)=jj 2168 do ii=boundmin(1),boundmax(1) 2169 coord(1)=ii 2170 2171 ! Here, apply the downsampling : skip some of the trials 2172 if(present(downsampling))then 2173 2174 if(downsampling(1)==0 .and. found1(coord(2),coord(3))==1)cycle 2175 if(downsampling(2)==0 .and. found2(coord(1),coord(3))==1)cycle 2176 if(downsampling(3)==0 .and. found3(coord(1),coord(2))==1)cycle 2177 2178 ads(:)=abs(downsampling(:)) 2179 if(ads(1)>0 .and. mod(coord(1),ads(1))/=0)cycle 2180 if(ads(2)>0 .and. mod(coord(2),ads(2))/=0)cycle 2181 if(ads(3)>0 .and. mod(coord(3),ads(2))/=0)cycle 2182 cds(:)=coord(:)/ads(:) 2183 if(minval(downsampling(:))<0)then ! If there is at least one negative number 2184 2185 if(downsampling(1)*downsampling(2)*downsampling(3)/=0)then ! If there is no zero number 2186 ! Face-centered case 2187 if(downsampling(1)<0 .and. downsampling(2)<0 .and. downsampling(3)<0)then ! All three are negative 2188 if(mod(sum(cds(:)),2)/=0)cycle 2189 ! One-face-centered case 2190 else if(downsampling(1)*downsampling(2)*downsampling(3)<0)then ! Only one is negative 2191 if(downsampling(1)<0 .and. mod(cds(2)+cds(3),2)/=0)cycle 2192 if(downsampling(2)<0 .and. mod(cds(1)+cds(3),2)/=0)cycle 2193 if(downsampling(3)<0 .and. mod(cds(1)+cds(2),2)/=0)cycle 2194 ! Body-centered case ! What is left : two are negative 2195 else 2196 ! Either all are zero, or all are one, so skip when sum is 1 or 2. 2197 if(sum(mod(cds(:),2))==1 .or. sum(mod(cds(:),2))==2)cycle 2198 end if 2199 else 2200 if(downsampling(1)==0 .and. mod(cds(2)+cds(3),2)/=0)cycle 2201 if(downsampling(2)==0 .and. mod(cds(1)+cds(3),2)/=0)cycle 2202 if(downsampling(3)==0 .and. mod(cds(1)+cds(2),2)/=0)cycle 2203 end if 2204 end if 2205 end if 2206 2207 do ikshft=1,nshiftk 2208 2209 ! Only the first shiftk is taken into account if downsampling 2210 ! if(.false.)then 2211 if(present(downsampling))then 2212 if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then 2213 if(ikshft>1)cycle 2214 end if 2215 end if 2216 2217 ! Coordinates of the trial k point with respect to the k primitive lattice 2218 k1(1)=ii+shiftk(1,ikshft) 2219 k1(2)=jj+shiftk(2,ikshft) 2220 k1(3)=kk+shiftk(3,ikshft) 2221 ! Reduced coordinates of the trial k point 2222 k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3) 2223 ! Eliminate the point if outside [0,1[ 2224 if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle 2225 if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle 2226 if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle 2227 ! Wrap the trial values in the interval ]-1/2,1/2] . 2228 call wrap2_pmhalf(k2(1),k1(1),shift) 2229 call wrap2_pmhalf(k2(2),k1(2),shift) 2230 call wrap2_pmhalf(k2(3),k1(3),shift) 2231 spkpt(:,nn)=k1(:) 2232 nn=nn+1 2233 2234 if(present(downsampling))then 2235 found1(coord(2),coord(3))=1 2236 found2(coord(1),coord(3))=1 2237 found3(coord(1),coord(2))=1 2238 end if 2239 2240 end do 2241 end do 2242 end do 2243 end do 2244 nkpt=nn-1 2245 2246 if(present(downsampling))then 2247 ABI_FREE(found1) 2248 ABI_FREE(found2) 2249 ABI_FREE(found3) 2250 end if 2251 2252 if(nkpt/=nkptlatt*nshiftk)then 2253 write(msg, '(a,i0,3a,i0,a)' )& 2254 'The number of k points ',nkpt,' is not equal to',ch10,& 2255 'nkptlatt*nshiftk which is ',nkptlatt*nshiftk,'.' 2256 ABI_BUG(msg) 2257 end if 2258 2259 else if(brav==2)then 2260 2261 ! Face-Centered Lattice 2262 if (prtvol > 0) call wrtout(std_out,' Face-Centered Lattice Grid ','COLL') 2263 if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2) then 2264 write(msg, '(a,a,a,i0,a,a,a,a,a)' )& 2265 & 'The value of mkpt is not large enough. It should be',ch10,& 2266 & 'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,',',ch10,& 2267 & 'Action: set mkpt to that value in the main routine,',ch10,& 2268 & 'and recompile the code.' 2269 ABI_BUG(msg) 2270 end if 2271 nn=1 2272 if (ngkpt(1)/=ngkpt(2).or.ngkpt(1)/=ngkpt(3)) then 2273 write(msg, '(4a,3(a,i0,a),a)' )& 2274 & 'For face-centered lattices, the numbers ngqpt(1:3)',ch10,& 2275 & 'must be equal, while they are :',ch10,& 2276 & 'ngqpt(1) = ',ngkpt(1),ch10,& 2277 & 'ngqpt(2) = ',ngkpt(2),ch10,& 2278 & 'ngqpt(3) = ',ngkpt(3),ch10,& 2279 & 'Action: modify ngqpt(1:3) in the input file.' 2280 ABI_BUG(msg) 2281 end if 2282 if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2)) then 2283 write(msg, '(4a,3(a,i0,a),a)' )& 2284 & 'For face-centered lattices, the numbers ngqpt(1:3)*nshiftk',ch10,& 2285 & 'must be even, while they are :',ch10,& 2286 & 'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,& 2287 & 'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,& 2288 & 'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,& 2289 & 'Action: modify ngqpt(1:3)*nshiftk in the input file.' 2290 ABI_ERROR(msg) 2291 end if 2292 if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then 2293 spkpt(1,1)=0.0_dp 2294 spkpt(2,1)=0.0_dp 2295 spkpt(3,1)=0.0_dp 2296 nkpt=1 2297 else 2298 do kk=1,ngkpt(3) 2299 do jj=1,ngkpt(2) 2300 do ii=1,ngkpt(1) 2301 do ikshft=1,nshiftk 2302 k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1) 2303 k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2) 2304 k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3) 2305 ! Wrap the trial values in the interval ]-1/2,1/2] . 2306 call wrap2_pmhalf(k1(1),k2(1),shift) 2307 call wrap2_pmhalf(k1(2),k2(2),shift) 2308 call wrap2_pmhalf(k1(3),k2(3),shift) 2309 ! Test whether it is inside the FCC BZ. 2310 ktest(1)=2*k2(1)-1.0d-10 2311 ktest(2)=2*k2(2)-2.0d-10 2312 ktest(3)=2*k2(3)-5.0d-10 2313 if (abs(ktest(1))+abs(ktest(2))+abs(ktest(3))<1.5_dp) then 2314 kcar(1)=ktest(1)+1.0d-10 2315 kcar(2)=ktest(2)+2.0d-10 2316 kcar(3)=ktest(3)+5.0d-10 2317 spkpt(1,nn)=0.5_dp*kcar(2)+0.5_dp*kcar(3) 2318 spkpt(2,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(3) 2319 spkpt(3,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(2) 2320 nn=nn+1 2321 end if 2322 end do 2323 end do 2324 end do 2325 end do 2326 nkpt=nn-1 2327 if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2)then 2328 write(msg, '(a,i8,a,a,a,i8,a)' )& 2329 & 'The number of k points ',nkpt,' is not equal to',ch10,& 2330 & '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2 which is',& 2331 & (ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,'.' 2332 ABI_BUG(msg) 2333 end if 2334 end if 2335 2336 else if(brav==3)then 2337 2338 ! Body-Centered Lattice (not mandatory cubic !) 2339 if (prtvol > 0) call wrtout(std_out,' Body-Centered Lattice Grid ','COLL') 2340 if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/4) then 2341 write(msg, '(a,a,a,i8,a,a,a,a,a)' )& 2342 & 'The value of mkpt is not large enough. It should be',ch10,& 2343 & 'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,',',ch10,& 2344 & 'Action: set mkpt to that value in the main routine,',ch10,& 2345 & 'and recompile the code.' 2346 ABI_BUG(msg) 2347 end if 2348 nn=1 2349 if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2) .or.& 2350 & (ngkpt(2)*nshiftk)/=(((ngkpt(2)*nshiftk)/2)*2) .or.& 2351 & (ngkpt(3)*nshiftk)/=(((ngkpt(3)*nshiftk)/2)*2) ) then 2352 write(msg, '(4a,3(a,i6,a),a)' )& 2353 & 'For body-centered lattices, the numbers ngqpt(1:3)',ch10,& 2354 & 'must be even, while they are :',ch10,& 2355 & 'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,& 2356 & 'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,& 2357 & 'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,& 2358 & 'Action: modify ngqpt(1:3) in the input file.' 2359 ABI_ERROR(msg) 2360 end if 2361 if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then 2362 spkpt(1,1)=0.0_dp 2363 spkpt(2,1)=0.0_dp 2364 spkpt(3,1)=0.0_dp 2365 nkpt=1 2366 else 2367 do kk=1,ngkpt(3) 2368 do jj=1,ngkpt(2) 2369 do ii=1,ngkpt(1) 2370 do ikshft=1,nshiftk 2371 k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1) 2372 k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2) 2373 k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3) 2374 ! Wrap the trial values in the interval ]-1/2,1/2] . 2375 call wrap2_pmhalf(k1(1),k2(1),shift) 2376 call wrap2_pmhalf(k1(2),k2(2),shift) 2377 call wrap2_pmhalf(k1(3),k2(3),shift) 2378 ! Test whether it is inside the BCC BZ. 2379 ktest(1)=2*k2(1)-1.0d-10 2380 ktest(2)=2*k2(2)-2.0d-10 2381 ktest(3)=2*k2(3)-5.0d-10 2382 if (abs(ktest(1))+abs(ktest(2))<1._dp) then 2383 if (abs(ktest(1))+abs(ktest(3))<1._dp) then 2384 if (abs(ktest(2))+abs(ktest(3))<1._dp) then 2385 kcar(1)=ktest(1)+1.0d-10 2386 kcar(2)=ktest(2)+2.0d-10 2387 kcar(3)=ktest(3)+5.0d-10 2388 spkpt(1,nn)=-0.5*kcar(1)+0.5*kcar(2)+0.5*kcar(3) 2389 spkpt(2,nn)=0.5*kcar(1)-0.5*kcar(2)+0.5*kcar(3) 2390 spkpt(3,nn)=0.5*kcar(1)+0.5*kcar(2)-0.5*kcar(3) 2391 nn=nn+1 2392 end if 2393 end if 2394 end if 2395 end do 2396 end do 2397 end do 2398 end do 2399 nkpt=nn-1 2400 if(nkpt==0)then 2401 write(msg, '(3a)' )& 2402 & 'BCC lattice, input ngqpt=0, so no kpt is generated.',ch10,& 2403 & 'Action: modify ngqpt(1:3) in the input file.' 2404 ABI_ERROR(msg) 2405 end if 2406 if(nkpt/=(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4)then 2407 write(msg, '(a,i0,3a,i0,a)' )& 2408 & 'The number of k points ',nkpt,' is not equal to',ch10,& 2409 & '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4 which is',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,'.' 2410 ABI_BUG(msg) 2411 end if 2412 end if 2413 2414 else if(brav==4)then 2415 2416 ! Hexagonal Lattice (D6h) 2417 if (prtvol > 0) call wrtout(std_out,' Hexagonal Lattice Grid ','COLL') 2418 if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)) then 2419 write(msg, '(a,a,a,i0,a,a,a,a,a)' )& 2420 & 'The value of mkpt is not large enough. It should be',ch10,& 2421 & 'at least',ngkpt(1)*ngkpt(2)*ngkpt(3),',',ch10,& 2422 & 'Action: set mkpt to that value in the main routine,',ch10,& 2423 & 'and recompile the code.' 2424 ABI_BUG(msg) 2425 end if 2426 nn=1 2427 if (ngkpt(1)/=ngkpt(2)) then 2428 write(msg, '(4a,2(a,i0,a),a)' )& 2429 & 'For hexagonal lattices, the numbers ngqpt(1:2)',ch10,& 2430 & 'must be equal, while they are:',ch10,& 2431 & 'ngqpt(1) = ',ngkpt(1),ch10,& 2432 & 'ngqpt(2) = ',ngkpt(2),ch10,& 2433 & 'Action: modify ngqpt(1:3) in the input file.' 2434 ABI_ERROR(msg) 2435 end if 2436 if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then 2437 write(msg, '(3a)' )& 2438 & 'For hexagonal lattices, ngqpt(1:3)=0 is not permitted',ch10,& 2439 & 'Action: modify ngqpt(1:3) in the input file.' 2440 ABI_ERROR(msg) 2441 else 2442 do kk=1,ngkpt(3) 2443 do jj=1,ngkpt(2) 2444 do ii=1,ngkpt(1) 2445 do ikshft=1,nshiftk 2446 k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1) 2447 k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2) 2448 k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3) 2449 ! Wrap the trial values in the interval ]-1/2,1/2] . 2450 call wrap2_pmhalf(k1(1),k2(1),shift) 2451 call wrap2_pmhalf(k1(2),k2(2),shift) 2452 call wrap2_pmhalf(k1(3),k2(3),shift) 2453 spkpt(:,nn)=k2(:) 2454 nn=nn+1 2455 end do 2456 end do 2457 end do 2458 end do 2459 nkpt=nn-1 2460 if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)then 2461 write(msg, '(a,i0,3a,i0,a)' )& 2462 & 'The number of k points ',nkpt,' is not equal to',ch10,& 2463 & 'ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk which is',ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk,'.' 2464 ABI_BUG(msg) 2465 end if 2466 end if 2467 2468 else 2469 2470 write(msg, '(a,i0,a,a,a)' )& 2471 & 'The calling routine asks brav= ',brav,'.',ch10,& 2472 & 'but only brav=1 or -1,2,3 or 4 are allowed.' 2473 ABI_BUG(msg) 2474 end if 2475 2476 if (option/=0) then 2477 ! Put the Gamma point first 2478 if(nkpt>1)then 2479 do ii=1,nkpt 2480 if(sum(abs(spkpt(:,ii)))<tol8)then 2481 spkpt(:,ii)=spkpt(:,1) 2482 spkpt(:,1)=zero 2483 exit 2484 end if 2485 end do 2486 end if 2487 2488 write(msg,'(a,i8)')' Grid q points : ',nkpt 2489 call wrtout(iout,msg,'COLL') 2490 nkpout=nkpt 2491 if(nkpt>80)then 2492 call wrtout(iout,' greater than 80, so only write 20 of them ','COLL') 2493 nkpout=20 2494 end if 2495 do ii=1,nkpout 2496 write(msg, '(1x,i2,a2,3es16.8)' )ii,') ',spkpt(1,ii),spkpt(2,ii),spkpt(3,ii) 2497 call wrtout(iout,msg,'COLL') 2498 end do 2499 end if 2500 2501 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
SOURCE
327 integer function symkchk(kptns,nkpt,nsym,symrec,timrev,errmsg) result(ierr) 328 329 !Arguments ------------------------------- 330 !scalars 331 integer,intent(in) :: nkpt,nsym,timrev 332 character(len=*),intent(out) :: errmsg 333 !arrays 334 integer,intent(in) :: symrec(3,3,nsym) 335 real(dp),intent(in) :: kptns(3,nkpt) 336 337 !Local variables ------------------------- 338 !scalars 339 integer :: identi,ii,ikpt,ikpt2,imatch,isym,jj,tident 340 real(dp) :: difk,reduce 341 character(len=500) :: msg 342 !arrays 343 real(dp) :: ksym(3) 344 345 ! ********************************************************************* 346 ierr = 0 347 348 if(timrev/=1 .and. timrev/=0)then 349 write(errmsg, '(3a,i0,a)' )& 350 'timrev should be 0 or 1, while',ch10,& 351 'it is equal to ',timrev,'.' 352 ierr = 1; return 353 end if 354 355 if(nsym/=1)then 356 ! Find the identity symmetry operation 357 do isym=1,nsym 358 tident=1 359 do jj=1,3 360 if(symrec(jj,jj,isym)/=1)tident=0 361 do ii=1,3 362 if( ii/=jj .and.symrec(ii,jj,isym)/=0)tident=0 363 end do 364 end do 365 if(tident==1)then 366 identi=isym 367 call wrtout(std_out,sjoin(' symkchk: found identity with number:', itoa(identi))) 368 exit 369 end if 370 end do 371 if(tident==0)then 372 errmsg = 'Did not found the identity operation.' 373 ierr = 1; return 374 end if 375 end if 376 377 !Here begins the serious business 378 !The length sorting, etc. of symkpt have been dropped because the 379 !computational cost is estimated to be negligible. 380 381 if(nsym>1 .or. timrev==1)then 382 383 ! Outer loop over kpts 384 do ikpt=1,nkpt-1 385 386 ! Loop on the symmetries 387 ! For each k-point and each symmetry transformation, a matching 388 ! k-point must be found, modulo time reversal if appropriate 389 do isym=1,nsym 390 391 ! Get the symmetric of the vector 392 do ii=1,3 393 ksym(ii)= kptns(1,ikpt)*symrec(ii,1,isym)& 394 & +kptns(2,ikpt)*symrec(ii,2,isym)& 395 & +kptns(3,ikpt)*symrec(ii,3,isym) 396 end do 397 398 ! Second loop k-points 399 do ikpt2=1,nkpt 400 401 ! Test for match of symmetric and any vector (including original) 402 imatch=1 403 do ii=1,3 404 difk= ksym(ii)-kptns(ii,ikpt2) 405 reduce=difk-anint(difk) 406 if(abs(reduce)>tol8)imatch=0 407 end do 408 if(imatch==1)exit 409 410 ! Test for match with time reversal 411 if(timrev==1)then 412 imatch=1 413 do ii=1,3 414 difk= ksym(ii)+kptns(ii,ikpt2) 415 reduce=difk-anint(difk) 416 if(abs(reduce)>tol8)imatch=0 417 end do 418 if(imatch==1)exit 419 end if 420 421 end do ! End secondary loop over k-points 422 if (imatch/=1) then 423 write(errmsg, '(a,a,a,i0,a,i0,a,a,a,a)' )& 424 'k-point set must have full space-group symmetry',ch10,& 425 'there is no match for kpt: ',ikpt,' transformed by symmetry: ',isym,ch10,& 426 'Action: change kptopt to 2 or 3 and/or change or use shiftk',ch10,& 427 'shiftk = 0 0 0 is always a safe choice.' 428 ierr = 2; return 429 end if 430 431 end do ! End loop on isym 432 end do ! End primary loop over k-points 433 434 write(msg,'(a)')' symkchk : k-point set has full space-group symmetry.' 435 call wrtout([std_out, ab_out], msg, 'COLL') 436 end if 437 438 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 symmetries
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,MAX_NSHIFTK)=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.
SOURCE
2541 subroutine testkgrid(bravais,iout,kptrlatt,kptrlen,& 2542 & msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum) 2543 2544 !Arguments ------------------------------------ 2545 !scalars 2546 integer,intent(in) :: iout,msym,nsym,prtkpt 2547 integer,intent(out) :: nshiftk 2548 real(dp),intent(inout) :: kptrlen 2549 !arrays 2550 integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3) 2551 integer,intent(out) :: kptrlatt(3,3) 2552 real(dp),intent(in) :: rprimd(3,3) 2553 real(dp),intent(inout) :: shiftk(3,MAX_NSHIFTK) !vz_i 2554 2555 !Local variables------------------------------- 2556 !scalars 2557 integer,parameter :: kptopt=1,mkpt_list=100000 2558 integer :: ang90,center,dirvacuum,equal,igrid,igrid_current,iholohedry,ii,init_mult,iscale,iscf 2559 integer :: iset,mult1,mult2,mult3,ndims,nkpt,nkpt_current,nkpt_trial,nset 2560 real(dp) :: buffer_scale,determinant,fact,factor,kptrlen_current,kptrlen_max,kptrlen_target 2561 real(dp) :: kptrlen_trial,length1,length2,length3,length_axis1,length_axis2 2562 real(dp) :: length_axis3,merit_factor,mult1h,mult2h,mult3h,reduceda,reducedb 2563 real(dp) :: sca,scb,scc,surface,ucvol 2564 character(len=500) :: msg 2565 !arrays 2566 integer :: kptrlatt_current(3,3),kptrlatt_trial(3,3) 2567 integer,allocatable :: grid_list(:) 2568 real(dp) :: axes(3,3),gmet(3,3),gprimd(3,3),matrix1(3,3),matrix2(3,3) 2569 real(dp) :: metmin(3,3),minim(3,3),r2d(3,3),rmet(3,3),rsuper(3,3) 2570 real(dp) :: shiftk_current(3,MAX_NSHIFTK),shiftk_trial(3,MAX_NSHIFTK) 2571 real(dp),allocatable :: kpt(:,:),kptrlen_list(:),wtk(:) 2572 2573 ! ************************************************************************* 2574 2575 kptrlen_target=kptrlen 2576 2577 !The vacuum array must be made of 0 or 1 2578 do ii=1,3 2579 if(vacuum(ii)/=0 .and. vacuum(ii)/=1)then 2580 write(msg,'(a,a,a,i1,a,i3,a,a)')& 2581 & 'The values of vacuum must be 0 or 1.',ch10,& 2582 & 'However, the input vacuum(',ii,') is',vacuum(ii),ch10,& 2583 & 'Action: correct vacuum in your input file.' 2584 ABI_ERROR(msg) 2585 end if 2586 end do 2587 2588 !Specific preparation for 2-dimensional system 2589 if(sum(vacuum(:))==1)then 2590 2591 ! Make the non-active vector orthogonal to the active vectors, 2592 ! and take it along the z direction 2593 if(vacuum(1)==1)then 2594 r2d(1,3)=rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3) 2595 r2d(2,3)=rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3) 2596 r2d(3,3)=rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3) 2597 r2d(:,1)=rprimd(:,2) 2598 r2d(:,2)=rprimd(:,3) 2599 dirvacuum=1 2600 else if(vacuum(2)==1)then 2601 r2d(1,3)=rprimd(2,3)*rprimd(3,1)-rprimd(3,3)*rprimd(2,1) 2602 r2d(2,3)=rprimd(3,3)*rprimd(1,1)-rprimd(1,3)*rprimd(3,1) 2603 r2d(3,3)=rprimd(1,3)*rprimd(2,1)-rprimd(2,3)*rprimd(1,1) 2604 r2d(:,1)=rprimd(:,3) 2605 r2d(:,2)=rprimd(:,1) 2606 dirvacuum=2 2607 else if(vacuum(3)==1)then 2608 r2d(1,3)=rprimd(2,1)*rprimd(3,2)-rprimd(3,1)*rprimd(2,2) 2609 r2d(2,3)=rprimd(3,1)*rprimd(1,2)-rprimd(1,1)*rprimd(3,2) 2610 r2d(3,3)=rprimd(1,1)*rprimd(2,2)-rprimd(2,1)*rprimd(1,2) 2611 r2d(:,1)=rprimd(:,1) 2612 r2d(:,2)=rprimd(:,2) 2613 dirvacuum=3 2614 end if 2615 surface=sqrt(sum(r2d(:,3)**2)) 2616 ! Identify the 2-D Bravais lattice 2617 ! DEBUG 2618 ! write(std_out,*)' r2d=',r2d(:,:) 2619 ! ENDDEBUG 2620 call metric(gmet,gprimd,-1,rmet,r2d,ucvol) 2621 call smallprim(metmin,minim,r2d) 2622 ! DEBUG 2623 ! write(std_out,*)' minim=',minim(:,:) 2624 ! ENDDEBUG 2625 ang90=0 ; equal=0 ; center=0 2626 axes(:,:)=minim(:,:) 2627 if(abs(metmin(1,2))<tol8)ang90=1 2628 if(abs(metmin(1,1)-metmin(2,2))<tol8)equal=1 2629 if(ang90==1)then 2630 if(equal==1)iholohedry=4 2631 if(equal==0)iholohedry=2 2632 else if(equal==1)then 2633 reduceda=metmin(1,2)/metmin(1,1) 2634 if(abs(reduceda+0.5_dp)<tol8)then 2635 iholohedry=3 2636 else if(abs(reduceda-0.5_dp)<tol8)then 2637 iholohedry=3 2638 ! Use conventional axes 2639 axes(:,2)=minim(:,2)-minim(:,1) 2640 else 2641 iholohedry=2 ; center=1 2642 axes(:,1)=minim(:,1)+minim(:,2) 2643 axes(:,2)=minim(:,2)-minim(:,1) 2644 end if 2645 else 2646 reduceda=metmin(1,2)/metmin(1,1) 2647 reducedb=metmin(1,2)/metmin(2,2) 2648 if(abs(reduceda+0.5_dp)<tol8)then 2649 iholohedry=2 ; center=1 2650 axes(:,2)=2.0_dp*minim(:,2)+minim(:,1) 2651 else if(abs(reduceda-0.5_dp)<tol8)then 2652 iholohedry=2 ; center=1 2653 axes(:,2)=2.0_dp*minim(:,2)-minim(:,1) 2654 else if(abs(reducedb+0.5_dp)<tol8)then 2655 iholohedry=2 ; center=1 2656 axes(:,1)=2.0_dp*minim(:,1)+minim(:,2) 2657 else if(abs(reducedb-0.5_dp)<tol8)then 2658 iholohedry=2 ; center=1 2659 axes(:,1)=2.0_dp*minim(:,1)-minim(:,2) 2660 else 2661 iholohedry=1 2662 end if 2663 end if 2664 ! Make sure that axes form a right-handed coordinate system 2665 determinant=axes(1,1)*axes(2,2)*axes(3,3) & 2666 & +axes(1,2)*axes(2,3)*axes(3,1) & 2667 & +axes(1,3)*axes(3,2)*axes(2,1) & 2668 & -axes(1,1)*axes(3,2)*axes(2,3) & 2669 & -axes(1,3)*axes(2,2)*axes(3,1) & 2670 & -axes(1,2)*axes(2,1)*axes(3,3) 2671 if(determinant<zero)then 2672 axes(:,1)=-axes(:,1) 2673 end if 2674 ! Prefer symmetry axes on the same side as the primitive axes 2675 sca=axes(1,1)*r2d(1,1)+axes(2,1)*r2d(2,1)+axes(3,1)*r2d(3,1) 2676 scb=axes(1,2)*r2d(1,2)+axes(2,2)*r2d(2,2)+axes(3,2)*r2d(3,2) 2677 scc=axes(1,3)*rprimd(1,dirvacuum)& 2678 & +axes(2,3)*rprimd(2,dirvacuum)& 2679 & +axes(3,3)*rprimd(3,dirvacuum) 2680 if(sca<-tol8 .and. scb<-tol8)then 2681 axes(:,1)=-axes(:,1) ; sca=-sca 2682 axes(:,2)=-axes(:,2) ; scb=-scb 2683 end if 2684 ! Doing this might change the angle between vectors, so that 2685 ! the cell is not conventional anymore 2686 ! if(sca<-tol8 .and. scc<-tol8)then 2687 ! axes(:,1)=-axes(:,1) ; sca=-sca 2688 ! axes(:,3)=-axes(:,3) ; scc=-scc 2689 ! end if 2690 ! if(scb<-tol8 .and. scc<-tol8)then 2691 ! axes(:,2)=-axes(:,2) ; scb=-scb 2692 ! axes(:,3)=-axes(:,3) ; scc=-scc 2693 ! end if 2694 length_axis1=sqrt(axes(1,1)**2+axes(2,1)**2+axes(3,1)**2) 2695 length_axis2=sqrt(axes(1,2)**2+axes(2,2)**2+axes(3,2)**2) 2696 2697 ! DEBUG 2698 ! write(std_out,*)' testkgrid: iholohedry, center =',iholohedry,center 2699 ! write(std_out,*)' testkgrid: axis 1=',axes(:,1) 2700 ! write(std_out,*)' testkgrid: axis 2=',axes(:,2) 2701 ! write(std_out,*)' testkgrid: axis 3=',axes(:,3) 2702 ! write(std_out,*)' testkgrid: length_axis=',length_axis1,length_axis2 2703 ! ENDDEBUG 2704 2705 ! End special treatment of 2-D case 2706 end if 2707 2708 !3-dimensional system 2709 if(sum(vacuum(:))==0)then 2710 iholohedry=bravais(1) 2711 center=bravais(2) 2712 fact=1.0_dp 2713 if(center/=0)fact=0.5_dp 2714 matrix1(:,1)=bravais(3:5)*fact 2715 matrix1(:,2)=bravais(6:8)*fact 2716 matrix1(:,3)=bravais(9:11)*fact 2717 call matr3inv(matrix1,matrix2) 2718 do ii=1,3 2719 axes(:,ii)=rprimd(:,1)*matrix2(ii,1)+rprimd(:,2)*matrix2(ii,2)+rprimd(:,3)*matrix2(ii,3) 2720 end do 2721 length_axis1=sqrt(axes(1,1)**2+axes(2,1)**2+axes(3,1)**2) 2722 length_axis2=sqrt(axes(1,2)**2+axes(2,2)**2+axes(3,2)**2) 2723 length_axis3=sqrt(axes(1,3)**2+axes(2,3)**2+axes(3,3)**2) 2724 ! DEBUG 2725 ! write(std_out,*)' testkgrid: axes=',axes(:,:) 2726 ! write(std_out,*)' length_axis=',length_axis1,length_axis2,length_axis3 2727 ! ENDDEBUG 2728 end if 2729 2730 !This routine examine only primitive k lattices. 2731 nshiftk=1 2732 2733 !If prtkpt/=0, will examine more grids than strictly needed 2734 buffer_scale=one 2735 if(prtkpt/=0)buffer_scale=two 2736 2737 if(prtkpt/=0)then 2738 write(msg,'(a,a,a,a,a,a,a,a)' )ch10,& 2739 ' testkgrid : will perform the analysis of a series of k-grids.',ch10,& 2740 ' Note that kptopt=1 in this analysis, irrespective of its input value.',ch10,ch10,& 2741 ' Grid# kptrlatt shiftk kptrlen nkpt iset',ch10 2742 call wrtout(std_out,msg,'COLL') 2743 call wrtout(iout,msg,'COLL') 2744 ABI_MALLOC(grid_list,(mkpt_list)) 2745 ABI_MALLOC(kptrlen_list,(mkpt_list)) 2746 grid_list(:)=0 2747 kptrlen_list(:)=0.0_dp 2748 end if 2749 2750 if(sum(vacuum(:))==3)then 2751 2752 kptrlatt(:,:)=0 2753 kptrlatt(1,1)=1 2754 kptrlatt(2,2)=1 2755 kptrlatt(3,3)=1 2756 shiftk(:,1)=0.0_dp 2757 kptrlen=1000.0_dp 2758 nkpt_current=1 2759 igrid_current=1 2760 2761 if(prtkpt/=0)then 2762 write(msg,& 2763 & '(a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )& 2764 & ' 1 ',kptrlatt(:,1),' ',shiftk(1,1),' ',kptrlen,1,1,ch10,& 2765 & ' ',kptrlatt(:,2),' ',shiftk(2,1),ch10,& 2766 & ' ',kptrlatt(:,3),' ',shiftk(3,1),ch10 2767 call wrtout(std_out,msg,'COLL') 2768 call wrtout(iout,msg,'COLL') 2769 ! The unit cell volume is fake 2770 ucvol=kptrlen**3 2771 end if 2772 2773 else 2774 2775 nkpt=0 ; nkpt_current=0 ; iscf=1 ; iset=1 2776 kptrlen_current=0.0_dp 2777 mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1 2778 ABI_MALLOC(kpt,(3,nkpt)) 2779 ABI_MALLOC(wtk,(nkpt)) 2780 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 2781 2782 ! Loop on different grids, the upper limit is only to avoid an infinite loop 2783 do igrid=1,1000 2784 2785 kptrlatt_trial(:,:)=0 2786 kptrlatt_trial(1,1)=1 2787 kptrlatt_trial(2,2)=1 2788 kptrlatt_trial(3,3)=1 2789 shiftk_trial(:,1)=0.0_dp 2790 2791 ! 1-dimensional system 2792 if(sum(vacuum(:))==2)then 2793 if(vacuum(1)==0)then 2794 kptrlatt_trial(1,1)=2*igrid ; shiftk_trial(1,1)=0.5_dp 2795 else if(vacuum(2)==0)then 2796 kptrlatt_trial(2,2)=2*igrid ; shiftk_trial(2,1)=0.5_dp 2797 else if(vacuum(3)==0)then 2798 kptrlatt_trial(3,3)=2*igrid ; shiftk_trial(3,1)=0.5_dp 2799 end if 2800 end if 2801 2802 ! 2-dimensional system 2803 if(sum(vacuum(:))==1)then 2804 2805 ! Treat hexagonal holohedries separately 2806 if(iholohedry==3)then 2807 2808 ! write(std_out,*)' testkgrid: 2D, hexagonal' 2809 2810 mult1=mult1+1 2811 nset=4 2812 if(iset==1)then 2813 rsuper(:,1)=axes(:,1)*mult1 2814 rsuper(:,2)=axes(:,2)*mult1 2815 shiftk_trial(:,1)=0.0_dp 2816 else if(iset==2)then 2817 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 2818 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 2819 shiftk_trial(1,1)=1.0_dp/3.0_dp 2820 shiftk_trial(2,1)=1.0_dp/3.0_dp 2821 else if(iset==3)then 2822 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 2823 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 2824 shiftk_trial(:,1)=0.0_dp 2825 else if(iset==4)then 2826 rsuper(:,1)=axes(:,1)*mult1 2827 rsuper(:,2)=axes(:,2)*mult1 2828 shiftk_trial(1,1)=0.5_dp 2829 shiftk_trial(2,1)=0.5_dp 2830 end if 2831 2832 else 2833 ! Now treat all other holohedries 2834 length1=length_axis1*mult1 2835 length2=length_axis2*mult2 2836 ! write(std_out,*)' testkgrid: (2d) length=',length1,length2 2837 if(abs(length1-length2)<tol8)then 2838 mult1=mult1+1 2839 mult2=mult2+1 2840 else if(length1>length2)then 2841 mult2=mult2+1 2842 else if(length2>length1)then 2843 mult1=mult1+1 2844 end if 2845 nset=4 2846 ! iset==5 and 6 are allowed only for centered lattice 2847 if(center==1)nset=6 2848 if(iset==1 .or. iset==2)then 2849 rsuper(:,1)=axes(:,1)*mult1 2850 rsuper(:,2)=axes(:,2)*mult2 2851 else if(iset==3 .or. iset==4)then 2852 rsuper(:,1)=axes(:,1)*mult1-axes(:,2)*mult2 2853 rsuper(:,2)=axes(:,1)*mult1+axes(:,2)*mult2 2854 else if(iset==5 .or. iset==6)then 2855 rsuper(:,1)=axes(:,1)*(mult1-0.5_dp)-axes(:,2)*(mult2-0.5_dp) 2856 rsuper(:,2)=axes(:,1)*(mult1-0.5_dp)+axes(:,2)*(mult2-0.5_dp) 2857 end if 2858 ! This was the easiest way to code all even mult1 and mult2 pairs: 2859 ! make separate series for this possibility. 2860 if(iset==2 .or. iset==4 .or. iset==6)then 2861 rsuper(:,1)=2.0_dp*rsuper(:,1) 2862 rsuper(:,2)=2.0_dp*rsuper(:,2) 2863 end if 2864 shiftk_trial(1,1)=0.5_dp 2865 shiftk_trial(2,1)=0.5_dp 2866 2867 end if 2868 2869 ! Put back the inactive direction 2870 if(dirvacuum==1)then 2871 rsuper(:,3)=rsuper(:,1) 2872 shiftk_trial(3,1)=shiftk_trial(1,1) 2873 rsuper(:,1)=rprimd(:,1) 2874 shiftk_trial(1,1)=0.0_dp 2875 else if(dirvacuum==2)then 2876 rsuper(:,3)=rsuper(:,1) 2877 shiftk_trial(3,1)=shiftk_trial(1,1) 2878 rsuper(:,1)=rsuper(:,2) 2879 shiftk_trial(1,1)=shiftk_trial(2,1) 2880 rsuper(:,2)=rprimd(:,2) 2881 shiftk_trial(2,1)=0.0_dp 2882 else if(dirvacuum==3)then 2883 rsuper(:,3)=rprimd(:,3) 2884 shiftk_trial(3,1)=0.0_dp 2885 end if 2886 2887 ! The supercell and the corresponding shift have been generated ! 2888 ! Convert cartesian coordinates into kptrlatt_trial 2889 do ii=1,3 2890 kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+& 2891 & gprimd(2,:)*rsuper(2,ii)+& 2892 & gprimd(3,:)*rsuper(3,ii) ) 2893 end do 2894 2895 ! End of 2-dimensional system 2896 end if 2897 2898 ! 3-dimensional system 2899 if(sum(vacuum(:))==0)then 2900 ! Treat hexagonal holohedries separately 2901 if(iholohedry==6)then 2902 length1=length_axis1*mult1 2903 length3=length_axis3*mult3 2904 ! write(std_out,*)' testkgrid: (hex) lengths=',length1,length2 2905 if(abs(length1-length3)<tol8)then 2906 mult1=mult1+1 2907 mult3=mult3+1 2908 else if(length1>length3)then 2909 mult3=mult3+1 2910 else if(length3>length1)then 2911 mult1=mult1+1 2912 end if 2913 nset=4 2914 if(iset==1)then 2915 rsuper(:,1)=axes(:,1)*mult1 2916 rsuper(:,2)=axes(:,2)*mult1 2917 rsuper(:,3)=axes(:,3)*mult3 2918 shiftk_trial(:,1)=0.0_dp 2919 shiftk_trial(3,1)=0.5_dp 2920 else if(iset==2)then 2921 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 2922 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 2923 rsuper(:,3)=axes(:,3)*mult3 2924 shiftk_trial(1,1)=1.0_dp/3.0_dp 2925 shiftk_trial(2,1)=1.0_dp/3.0_dp 2926 shiftk_trial(3,1)=0.5_dp 2927 else if(iset==3)then 2928 rsuper(:,1)=(axes(:,1)-axes(:,2)) *mult1 2929 rsuper(:,2)=(axes(:,1)+2*axes(:,2))*mult1 2930 rsuper(:,3)=axes(:,3)*mult3 2931 shiftk_trial(:,1)=0.0_dp 2932 shiftk_trial(3,1)=0.5_dp 2933 else if(iset==4)then 2934 rsuper(:,1)=axes(:,1)*mult1 2935 rsuper(:,2)=axes(:,2)*mult1 2936 rsuper(:,3)=axes(:,3)*mult3 2937 shiftk_trial(:,1)=0.5_dp 2938 end if 2939 2940 else 2941 ! Now treat all other holohedries 2942 length1=length_axis1*mult1 2943 length2=length_axis2*mult2 2944 length3=length_axis3*mult3 2945 ! write(std_out,*)' testkgrid: length=',length1,length2,length3 2946 if(length2>length1+tol8 .and. length3>length1+tol8)then 2947 mult1=mult1+1 2948 else if(length1>length2+tol8 .and. length3>length2+tol8)then 2949 mult2=mult2+1 2950 else if(length1>length3+tol8 .and. length2>length3+tol8)then 2951 mult3=mult3+1 2952 else if(abs(length2-length3)<tol8 .and. & 2953 & abs(length1-length3)<tol8 .and. & 2954 & abs(length1-length2)<tol8 )then 2955 mult1=mult1+1 ; mult2=mult2+1 ; mult3=mult3+1 2956 else if(abs(length1-length2)<tol8)then 2957 mult1=mult1+1 ; mult2=mult2+1 2958 else if(abs(length1-length3)<tol8)then 2959 mult1=mult1+1 ; mult3=mult3+1 2960 else if(abs(length2-length3)<tol8)then 2961 mult2=mult2+1 ; mult3=mult3+1 2962 end if 2963 nset=6 2964 if(center==-1 .or. center==-3)nset=8 2965 if(iset==1 .or. iset==2)then 2966 ! Simple lattice of k points 2967 rsuper(:,1)=axes(:,1)*mult1 2968 rsuper(:,2)=axes(:,2)*mult2 2969 rsuper(:,3)=axes(:,3)*mult3 2970 shiftk_trial(:,1)=0.5_dp 2971 else if(iset==3 .or. iset==4)then 2972 ! FCC lattice of k points = BCC lattice in real space 2973 rsuper(:,1)=-axes(:,1)*mult1+axes(:,2)*mult2+axes(:,3)*mult3 2974 rsuper(:,2)= axes(:,1)*mult1-axes(:,2)*mult2+axes(:,3)*mult3 2975 rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2-axes(:,3)*mult3 2976 shiftk_trial(:,1)=0.5_dp 2977 else if(iset==5 .or. iset==6)then 2978 ! BCC lattice of k points = FCC lattice in real space 2979 rsuper(:,1)= axes(:,2)*mult2+axes(:,3)*mult3 2980 rsuper(:,2)= axes(:,1)*mult1 +axes(:,3)*mult3 2981 rsuper(:,3)= axes(:,1)*mult1+axes(:,2)*mult2 2982 ! The BCC lattice has no empty site with full symmetry 2983 shiftk_trial(:,1)=0.0_dp 2984 else if(iset==7 .or. iset==8)then 2985 ! iset==7 and 8 are allowed only for centered lattice 2986 mult1h=mult1-0.5_dp 2987 mult2h=mult2-0.5_dp 2988 mult3h=mult3-0.5_dp 2989 if(center==-1)then 2990 ! FCC lattice of k points = BCC lattice in real space 2991 rsuper(:,1)=-axes(:,1)*mult1h+axes(:,2)*mult2h+axes(:,3)*mult3h 2992 rsuper(:,2)= axes(:,1)*mult1h-axes(:,2)*mult2h+axes(:,3)*mult3h 2993 rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h-axes(:,3)*mult3h 2994 shiftk_trial(:,1)=0.5_dp 2995 else if(center==-3)then 2996 ! BCC lattice of k points = FCC lattice in real space 2997 rsuper(:,1)= axes(:,2)*mult2h+axes(:,3)*mult3h 2998 rsuper(:,2)= axes(:,1)*mult1h +axes(:,3)*mult3h 2999 rsuper(:,3)= axes(:,1)*mult1h+axes(:,2)*mult2h 3000 ! The BCC lattice has no empty site with full symmetry 3001 shiftk_trial(:,1)=0.0_dp 3002 end if 3003 end if 3004 ! This was the easiest way to code all even mult1, mult2, mult3 triplets: 3005 ! make separate series for this possibility. 3006 if(2*(iset/2)==iset)then 3007 rsuper(:,1)=2.0_dp*rsuper(:,1) 3008 rsuper(:,2)=2.0_dp*rsuper(:,2) 3009 rsuper(:,3)=2.0_dp*rsuper(:,3) 3010 end if 3011 end if 3012 3013 ! write(std_out,*)' testkgrid: gprimd=',gprimd(:,:) 3014 ! write(std_out,*)' testkgrid: rsuper=',rsuper(:,:) 3015 ! write(std_out,*)' testkgrid: iset =',iset 3016 3017 ! The supercell and the corresponding shift have been generated! 3018 ! Convert cartesian coordinates into kptrlatt_trial 3019 do ii=1,3 3020 kptrlatt_trial(:,ii)=nint( gprimd(1,:)*rsuper(1,ii)+& 3021 & gprimd(2,:)*rsuper(2,ii)+& 3022 & gprimd(3,:)*rsuper(3,ii) ) 3023 end do 3024 3025 ! End of 3-dimensional system 3026 end if 3027 3028 ! write(std_out,*)' testkgrid: before getkgrid' 3029 ! write(std_out,*)' testkgrid: rprimd=',rprimd(:,:) 3030 ! write(std_out,*)' testkgrid: kptrlatt_trial=',kptrlatt_trial(:,:) 3031 3032 call getkgrid(0,0,iscf,kpt,& 3033 & kptopt,kptrlatt_trial,kptrlen_trial,& 3034 & msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,& 3035 & shiftk_trial,symafm,symrel,vacuum,wtk) 3036 3037 ! write(std_out,*)' testkgrid: after getkgrid' 3038 3039 ! 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, 3040 ! that generates a kptrlen_trial that is just below kptrlen. 3041 if(prtkpt==0 .and. init_mult==1 .and. kptrlen_trial<(half-tol8)*kptrlen )then 3042 iscale=int((one-tol8)*kptrlen/kptrlen_trial) 3043 mult1=mult1*iscale 3044 mult2=mult2*iscale 3045 mult3=mult3*iscale 3046 init_mult=0 3047 ! write(std_out,*)' testkgrid: iscale=',iscale 3048 kptrlatt_trial(:,:)=kptrlatt_trial(:,:)*iscale 3049 call getkgrid(0,0,iscf,kpt,& 3050 & kptopt,kptrlatt_trial,kptrlen_trial,& 3051 & msym,nkpt,nkpt_trial,nshiftk,nsym,rprimd,& 3052 & shiftk_trial,symafm,symrel,vacuum,wtk) 3053 end if 3054 3055 if( (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_current==0) .or. & 3056 & (kptrlen_trial+tol8>kptrlen*(1.0_dp+tol8) .and. nkpt_trial<nkpt_current) .or. & 3057 & (nkpt_trial==nkpt_current .and. kptrlen_trial>kptrlen_current*(1.0_dp+tol8)))then 3058 3059 kptrlatt_current(:,:)=kptrlatt_trial(:,:) 3060 nkpt_current=nkpt_trial 3061 shiftk_current(:,:)=shiftk_trial(:,:) 3062 kptrlen_current=kptrlen_trial 3063 igrid_current=igrid 3064 end if 3065 3066 if(prtkpt/=0)then 3067 write(msg,'(i5,a,3i4,a,es14.4,a,es14.4,i8,i6,a,a,3i4,a,es14.4,a,a,3i4,a,es14.4,a)' )& 3068 & igrid,' ',kptrlatt_trial(:,1),' ',shiftk_trial(1,1),& 3069 & ' ',kptrlen_trial,nkpt_trial,iset,ch10,& 3070 & ' ',kptrlatt_trial(:,2),' ',shiftk_trial(2,1),ch10,& 3071 & ' ',kptrlatt_trial(:,3),' ',shiftk_trial(3,1),ch10 3072 call wrtout(std_out,msg,'COLL') 3073 call wrtout(iout,msg,'COLL') 3074 3075 ! Keep track of this grid, if it is worth 3076 if(kptrlen_trial > kptrlen_list(nkpt_trial)*(1.0_dp+tol8))then 3077 grid_list(nkpt_trial)=igrid 3078 kptrlen_list(nkpt_trial)=kptrlen_trial 3079 end if 3080 end if 3081 3082 ! Treat 1-D case 3083 if( sum(vacuum(:))==2 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )exit 3084 3085 ! Treat 2-D case or 3-D case 3086 if( sum(vacuum(:))<=1 .and. kptrlen_trial>buffer_scale*(1.0_dp+tol8)*kptrlen )then 3087 ! The present set of sets of k points is finished: 3088 ! either it was the last, or one has to go to the next one 3089 if(iset==nset)exit 3090 iset=iset+1 3091 mult1=0 ; mult2=0 ; mult3=0 ; init_mult=1 3092 end if 3093 3094 end do ! igrid=1,1000 3095 3096 ABI_FREE(kpt) 3097 ABI_FREE(wtk) 3098 3099 kptrlatt(:,:)=kptrlatt_current(:,:) 3100 shiftk(:,:)=shiftk_current(:,:) 3101 kptrlen=kptrlen_current 3102 3103 end if ! test on the number of dimensions 3104 3105 if(prtkpt/=0)then 3106 3107 ! sqrt(1/2) comes from the FCC packing, the best one 3108 factor=sqrt(0.5_dp)/ucvol/dble(nsym) 3109 ndims=3 3110 if(sum(vacuum(:))/=0)then 3111 if(sum(vacuum(:))==1)then 3112 ! sqrt(3/4) comes from the hex packing, the best one 3113 ! one multiplies by 2 because nsym is likely twice the number 3114 ! of symmetries that can be effectively used in 2D 3115 ndims=2 ; factor=sqrt(0.75_dp)/surface/dble(nsym)*2 3116 write(msg,'(2a)' )ch10,' Note that the system is bi-dimensional.' 3117 else if(sum(vacuum(:))==2)then 3118 ndims=1 ; factor=1/ucvol 3119 write(msg,'(2a)' )ch10,' Note that the system is uni-dimensional.' 3120 else if(sum(vacuum(:))==3)then 3121 ndims=0 3122 write(msg,'(2a)' )ch10,' Note that the system is zero-dimensional.' 3123 end if 3124 call wrtout(std_out,msg,'COLL') 3125 call wrtout(iout,msg,'COLL') 3126 end if 3127 3128 ! The asymptotic value of the merit factor is determined 3129 ! by the set of symmetries: in 3D, if it includes the 3130 ! inversion symmetry, the limit will be 1, if not, it 3131 ! will be two. In 2D, if it includes the inversion symmetry 3132 ! and an operation that maps z on -z, it will tend to one, 3133 ! while if only one of these operations is present, 3134 ! it will tend to two, and if none is present, it will tend to four. 3135 write(msg,'(11a)' )ch10,& 3136 & ' List of best grids, ordered by nkpt.',ch10,& 3137 & ' (stop at a value of kptrlen 20% larger than the target value).',ch10,& 3138 & ' (the merit factor will tend to one or two in 3 dimensions)',ch10,& 3139 & ' (and to one, two or four in 2 dimensions)',ch10,ch10,& 3140 & ' nkpt kptrlen grid# merit_factor' 3141 call wrtout(std_out,msg,'COLL') 3142 call wrtout(iout,msg,'COLL') 3143 3144 kptrlen_max=0.0_dp 3145 do ii=1,mkpt_list 3146 if(kptrlen_list(ii)>kptrlen_max*(1.0_dp+tol8))then 3147 kptrlen_max=kptrlen_list(ii) 3148 merit_factor=kptrlen_max**ndims/dble(ii)*factor 3149 write(msg, '(i6,es14.4,i6,f12.4)' )ii,kptrlen_max,grid_list(ii),merit_factor 3150 call wrtout(std_out,msg,'COLL') 3151 call wrtout(iout,msg,'COLL') 3152 end if 3153 if(kptrlen_max>1.2_dp*(1.0_dp-tol8)*kptrlen_target)exit 3154 end do 3155 3156 write(msg,'(a,a,es14.4,a,a,i6,a,a,a,es14.4,a,i6)' )ch10,& 3157 & ' For target kptrlen=',kptrlen_target,',',& 3158 & ' the selected grid is number',igrid_current,',',ch10,& 3159 & ' giving kptrlen=',kptrlen_current,' with nkpt=',nkpt_current 3160 call wrtout(std_out,msg,'COLL') 3161 call wrtout(iout,msg,'COLL') 3162 3163 write(msg,'(a,a,a,a)' )ch10,& 3164 & ' testkgrid : stop after analysis of a series of k-grids.',ch10,& 3165 & ' For usual production runs, set prtkpt back to 0 (the default).' 3166 call wrtout(std_out,msg,'COLL',do_flush=.True.) 3167 call wrtout(iout,msg,'COLL',do_flush=.True.) 3168 3169 call abi_abort('PERS',exit_status=0,print_config=.false.) 3170 end if 3171 3172 end subroutine testkgrid
m_kpts/tetra_from_kptrlatt [ Functions ]
[ Top ] [ m_kpts ] [ Functions ]
NAME
tetra_from_kptrlatt
FUNCTION
Helper function to to create an instance and htetra 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. comm= MPI communicator
OUTPUT
tetra<htetra_t>=Tetrahedron object, fully initialized if ierr == 0. msg=Error message if ierr /= 0 ierr=Exit status
SOURCE
216 type(htetra_t) function tetra_from_kptrlatt( & 217 cryst, kptopt, kptrlatt, nshiftk, shiftk, nkibz, kibz, comm, msg, ierr) result (htetra) 218 219 !Arguments ------------------------------------ 220 !scalars 221 integer,intent(in) :: kptopt,nshiftk,nkibz,comm 222 integer,intent(out) :: ierr 223 character(len=*),intent(out) :: msg 224 type(crystal_t),intent(in) :: cryst 225 !arrays 226 integer,intent(in) :: kptrlatt(3,3) 227 real(dp),intent(in) :: shiftk(3,nshiftk),kibz(3,nkibz) 228 229 !Local variables------------------------------- 230 !scalars 231 integer :: nkfull,my_nkibz,new_nshiftk 232 character(len=80) :: errorstring 233 !arrays 234 integer :: new_kptrlatt(3,3) 235 integer,allocatable :: indkk(:) 236 integer,allocatable :: bz2ibz(:,:) 237 real(dp) :: rlatt(3,3),klatt(3,3) 238 real(dp),allocatable :: kfull(:,:),my_kibz(:,:),my_wtk(:),new_shiftk(:,:) 239 240 ! ************************************************************************* 241 242 ierr = 0 243 244 ! Refuse only 1 kpoint: the algorithms are no longer valid. DOH! 245 if (nkibz == 1) then 246 msg = 'You need at least 2 kpoints to use the tetrahedron method.' 247 ierr = 1; goto 10 248 end if 249 if (all(kptrlatt == 0)) then 250 msg = 'Cannot generate tetrahedron because input kptrlatt == 0.' 251 ierr = 1; goto 10 252 end if 253 if (kptopt <= 0) then 254 msg = sjoin("Cannot generate tetrahedron because input kptopt:", itoa(kptopt)) 255 ierr = 1; goto 10 256 end if 257 258 call kpts_ibz_from_kptrlatt(cryst, kptrlatt, kptopt, nshiftk, shiftk, & 259 my_nkibz, my_kibz, my_wtk, nkfull, kfull, new_kptrlatt=new_kptrlatt, new_shiftk=new_shiftk, bz2ibz=bz2ibz) 260 261 ABI_FREE(my_wtk) 262 new_nshiftk = size(new_shiftk, dim=2) 263 264 if (my_nkibz /= nkibz .or. all(my_kibz /= kibz) ) then 265 msg = sjoin("Input nkibz:", itoa(nkibz), "does not agree with computed value:", itoa(my_nkibz)) 266 ierr = 1; goto 10 267 end if 268 269 ! Do not support new_nshiftk > 1: lattice must be decomposed into boxes 270 ! and this is not always possible (I think) with bizarre shifts 271 ! normally at this point we have incorporated everything into 272 ! new_kptrlatt, and only 1 shift is needed (in particular for MP grids). 273 if (new_nshiftk > 1) then 274 write(msg, "(9a)") & 275 'Cannot create tetrahedron object...',ch10, & 276 'Only simple lattices are supported. Action: use nshiftk=1.',ch10, & 277 'new_shiftk: ', trim(ltoa(reshape(new_shiftk, [3*new_nshiftk]))),ch10, & 278 'new_kptrlatt: ', trim(ltoa(reshape(new_kptrlatt, [9]))) 279 ierr = 2; goto 10 280 end if 281 282 rlatt = new_kptrlatt; call matr3inv(rlatt, klatt) 283 284 ABI_MALLOC(indkk, (nkfull)) 285 indkk(:) = bz2ibz(1, :) 286 ABI_SFREE(bz2ibz) 287 288 call htetra_init(htetra, indkk, cryst%gprimd, klatt, kfull, nkfull, my_kibz, my_nkibz, ierr, errorstring, comm) 289 if (ierr /= 0) msg = errorstring 290 291 10 continue 292 293 ABI_SFREE(my_kibz) 294 ABI_SFREE(indkk) 295 ABI_SFREE(kfull) 296 ABI_SFREE(new_shiftk) 297 298 end function tetra_from_kptrlatt