TABLE OF CONTENTS
ABINIT/smpbz [ 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.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (JCC,XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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) J.D. Pack and H.J. Monkhorst, Phys. Rev. B 16, 1748 (1977) A.H. MacDonald, Phys. Rev. B 18, 5897 (1978) R.A. Evarestov and V.P. Smirnov, Phys. Stat. Sol. (b) 119, 9 (1983)
PARENTS
ep_setupqpt,getkgrid,harmonic_thermo,initberry,initorbmag,m_fstab,m_ifc m_tdep_abitypes
CHILDREN
matr3inv,wrap2_pmhalf,wrtout
SOURCE
71 #if defined HAVE_CONFIG_H 72 #include "config.h" 73 #endif 74 75 #include "abi_common.h" 76 77 subroutine smpbz(brav,iout,kptrlatt,mkpt,nkpt,nshiftk,option,shiftk,spkpt,downsampling) 78 79 use defs_basis 80 use m_errors 81 use m_profiling_abi 82 83 use m_numeric_tools, only : wrap2_pmhalf 84 85 !This section has been created automatically by the script Abilint (TD). 86 !Do not modify the following lines by hand. 87 #undef ABI_FUNC 88 #define ABI_FUNC 'smpbz' 89 use interfaces_14_hidewrite 90 use interfaces_32_util 91 !End of the abilint section 92 93 implicit none 94 95 !Arguments ------------------------------- 96 !scalars 97 integer,intent(in) :: brav,iout,mkpt,nshiftk,option 98 integer,intent(out) :: nkpt 99 !arrays 100 integer,intent(in) :: kptrlatt(3,3) 101 integer,optional,intent(in) :: downsampling(3) 102 real(dp),intent(in) :: shiftk(3,nshiftk) 103 real(dp),intent(out) :: spkpt(3,mkpt) 104 105 !Local variables ------------------------- 106 !scalars 107 integer,parameter :: prtvol=0 108 integer :: dividedown,ii,ikshft,jj,kk,nkpout,nkptlatt,nn,proddown 109 real(dp) :: shift 110 character(len=500) :: message 111 !arrays 112 integer :: ads(3),boundmax(3),boundmin(3),cds(3),coord(3),ngkpt(3) 113 integer, allocatable :: found1(:,:),found2(:,:),found3(:,:) 114 real(dp) :: k1(3),k2(3),kcar(3),klatt(3,3),ktest(3),rlatt(3,3) 115 116 ! ********************************************************************* 117 118 !DEBUG 119 !write(std_out,*)' smpbz : brav,iout,mkpt,nkpt,option=',brav,iout,mkpt,nkpt,option 120 !write(std_out,*)' smpbz : kptrlatt(:,:)=',kptrlatt(:,:) 121 !write(std_out,*)' smpbz : nshiftk=',nshiftk 122 !write(std_out,*)' smpbz : shiftk(:,:)=',shiftk(:,:) 123 !write(std_out,*)' smpbz : downsampling(:)=',downsampling(:) 124 !ENDDEBUG 125 126 if(option/=0)then 127 call wrtout(iout,' Homogeneous q point set in the B.Z. ','COLL') 128 end if 129 130 if(abs(brav)/=1)then 131 ! Only generate Monkhorst-Pack lattices 132 if(kptrlatt(1,2)/=0 .or. kptrlatt(2,1)/=0 .or. & 133 & kptrlatt(1,3)/=0 .or. kptrlatt(3,1)/=0 .or. & 134 & kptrlatt(2,3)/=0 .or. kptrlatt(3,2)/=0 ) then 135 write(message, '(2a,a,3i4,a,a,3i4,a,a,3i4)' )& 136 & 'When abs(brav)/=1, kptrlatt must be diagonal, while it is',ch10,& 137 & 'kptrlatt(:,1)=',kptrlatt(:,1),ch10,& 138 & 'kptrlatt(:,2)=',kptrlatt(:,2),ch10,& 139 & 'kptrlatt(:,3)=',kptrlatt(:,3) 140 MSG_BUG(message) 141 end if 142 143 ngkpt(1)=kptrlatt(1,1) 144 ngkpt(2)=kptrlatt(2,2) 145 ngkpt(3)=kptrlatt(3,3) 146 ! 147 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 148 write(message, '(5a,i4,a,a,i4,a,a,i4,a,a)' )& 149 & 'All ngkpt (or ngqpt) must be strictly positive',ch10,& 150 & 'or all ngk(q)pt must be zero (for Gamma sampling), but :',ch10,& 151 & 'ngk(q)pt(1) = ',ngkpt(1),ch10,& 152 & 'ngk(q)pt(2) = ',ngkpt(2),ch10,& 153 & 'ngk(q)pt(3) = ',ngkpt(3),ch10,& 154 & 'Action: correct ngkpt or ngqpt in the input file.' 155 MSG_BUG(message) 156 end if 157 end if 158 159 !Just in case the user wants the grid downsampled to the Gamma point, checks that it is present, and possibly exits 160 if(present(downsampling))then 161 if(sum(abs(downsampling(:)))==0)then 162 do ikshft=1,nshiftk 163 if(sum(abs(shiftk(:,ikshft)))>tol12)cycle 164 nkpt=1 165 spkpt(:,1)=zero 166 return 167 end do 168 end if 169 end if 170 171 !********************************************************************* 172 173 if(abs(brav)==1)then 174 175 ! Compute the number of k points in the G-space unit cell 176 ! (will be multiplied by nshiftk later). 177 nkptlatt=kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3) & 178 & +kptrlatt(1,2)*kptrlatt(2,3)*kptrlatt(3,1) & 179 & +kptrlatt(1,3)*kptrlatt(2,1)*kptrlatt(3,2) & 180 & -kptrlatt(1,2)*kptrlatt(2,1)*kptrlatt(3,3) & 181 & -kptrlatt(1,3)*kptrlatt(2,2)*kptrlatt(3,1) & 182 & -kptrlatt(1,1)*kptrlatt(2,3)*kptrlatt(3,2) 183 184 if(present(downsampling))then 185 if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then 186 if(nshiftk>1)then 187 write(message, '(a,3i4,2a,i4,4a)' )& 188 & 'Real downsampling is activated, with downsampling(1:3)=',downsampling(1:3),ch10,& 189 & 'However, nshiftk must be 1 in this case, while the input nshiftk=',nshiftk,ch10,& 190 & 'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,& 191 & 'or set nshiftk=1.' 192 MSG_ERROR(message) 193 end if 194 proddown=downsampling(1)*downsampling(2)*downsampling(3) 195 if(proddown/=0)then 196 dividedown=abs(proddown) 197 if(minval(downsampling(:))<0)then ! If there is at least one negative number 198 dividedown=dividedown*2 199 if(proddown>0)dividedown=dividedown*2 ! If there are two negative numbers 200 end if 201 end if 202 if(mod(nkptlatt,dividedown)==0)then 203 nkptlatt=nkptlatt/dividedown 204 else 205 write(message, '(a,3i4,2a,i4,4a)' )& 206 & 'The requested downsampling, with downsampling(1:3)=',downsampling(1:3),ch10,& 207 & 'is not compatible with kptrlatt=',ch10,& 208 & kptrlatt(:,:),ch10,& 209 & 'that gives nkptlatt=',nkptlatt,ch10,& 210 & 'Action: either choose not to downsample the k point grid (e.g. fockdownsampling=1),',ch10,& 211 & 'or modify your k-point grid and/or your downsampling in order for them to be compatible.' 212 MSG_ERROR(message) 213 end if 214 end if 215 end if 216 217 ! Simple Lattice 218 if (prtvol > 0) call wrtout(std_out,' Simple Lattice Grid ','COLL') 219 if (mkpt<nkptlatt*nshiftk) then 220 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 221 & 'The value of mkpt is not large enough. It should be',ch10,& 222 & 'at least',nkptlatt*nshiftk,',',ch10,& 223 & 'Action: set mkpt to that value in the main routine,',ch10,& 224 & 'and recompile the code.' 225 MSG_BUG(message) 226 end if 227 228 ! Build primitive vectors of the k lattice 229 rlatt(:,:)=kptrlatt(:,:) 230 call matr3inv(rlatt,klatt) 231 232 !DEBUG 233 ! write(std_out,*)' First primitive vector of the k lattice :',klatt(:,1) 234 ! write(std_out,*)' Second primitive vector of the k lattice :',klatt(:,2) 235 ! write(std_out,*)' Third primitive vector of the k lattice :',klatt(:,3) 236 !ENDDEBUG 237 238 ! Now, klatt contains the three primitive vectors of the k lattice, 239 ! in reduced coordinates. One builds all k vectors that 240 ! are contained in the first Brillouin zone, with coordinates 241 ! in the interval [0,1[ . First generate boundaries of a big box. 242 243 do jj=1,3 244 245 ! Mathematically, one has to find the coordinates of the corners of a 246 ! rectangular paralleliped with integer coordinates, that multiplies the klatt primitive cell and allows 247 ! it to incorporate completely the [0,1]^3 box. Then take the minimum and maximum 248 ! of these coordinates, and round them negatively and positively to the next integer. 249 ! This can be done easily using kptrlatt, considering each coordinate in turn 250 ! and boils down to enlarging the boundaries for jj by the value of kptrlatt(:,jj), 251 ! acting on boundmin or boundmax depending on the sign ot kptrlatt(:,jj). 252 ! XG171020 The coding before 171020 was correct, despite being very simple. 253 boundmin(jj)=0 ; boundmax(jj)=0 254 do ii=1,3 255 if(kptrlatt(ii,jj)<0)boundmin(jj)=boundmin(jj)+kptrlatt(ii,jj) 256 if(kptrlatt(ii,jj)>0)boundmax(jj)=boundmax(jj)+kptrlatt(ii,jj) 257 end do 258 259 ! To accomodate the shifts, boundmin and boundmax don't start from 0, but are enlarged by one 260 ! positively and/or negatively. 261 ! XG171020 Coding in v8.6.0 and before was not correct. This one is even simpler actually. 262 boundmin(jj)=boundmin(jj)-ceiling(maxval(shiftk(jj,:))+tol14) 263 boundmax(jj)=boundmax(jj)-floor(minval(shiftk(jj,:))-tol14) 264 265 end do 266 267 if(present(downsampling))then 268 ABI_ALLOCATE(found1,(boundmin(2):boundmax(2),boundmin(3):boundmax(3))) 269 ABI_ALLOCATE(found2,(boundmin(1):boundmax(1),boundmin(3):boundmax(3))) 270 ABI_ALLOCATE(found3,(boundmin(1):boundmax(1),boundmin(2):boundmax(2))) 271 found1=0 ; found2=0 ; found3=0 272 end if 273 274 nn=1 275 do kk=boundmin(3),boundmax(3) 276 coord(3)=kk 277 do jj=boundmin(2),boundmax(2) 278 coord(2)=jj 279 do ii=boundmin(1),boundmax(1) 280 coord(1)=ii 281 282 ! Here, apply the downsampling : skip some of the trials 283 if(present(downsampling))then 284 285 if(downsampling(1)==0 .and. found1(coord(2),coord(3))==1)cycle 286 if(downsampling(2)==0 .and. found2(coord(1),coord(3))==1)cycle 287 if(downsampling(3)==0 .and. found3(coord(1),coord(2))==1)cycle 288 289 ads(:)=abs(downsampling(:)) 290 if(ads(1)>0 .and. mod(coord(1),ads(1))/=0)cycle 291 if(ads(2)>0 .and. mod(coord(2),ads(2))/=0)cycle 292 if(ads(3)>0 .and. mod(coord(3),ads(2))/=0)cycle 293 cds(:)=coord(:)/ads(:) 294 if(minval(downsampling(:))<0)then ! If there is at least one negative number 295 296 if(downsampling(1)*downsampling(2)*downsampling(3)/=0)then ! If there is no zero number 297 ! Face-centered case 298 if(downsampling(1)<0 .and. downsampling(2)<0 .and. downsampling(3)<0)then ! All three are negative 299 if(mod(sum(cds(:)),2)/=0)cycle 300 ! One-face-centered case 301 else if(downsampling(1)*downsampling(2)*downsampling(3)<0)then ! Only one is negative 302 if(downsampling(1)<0 .and. mod(cds(2)+cds(3),2)/=0)cycle 303 if(downsampling(2)<0 .and. mod(cds(1)+cds(3),2)/=0)cycle 304 if(downsampling(3)<0 .and. mod(cds(1)+cds(2),2)/=0)cycle 305 ! Body-centered case ! What is left : two are negative 306 else 307 if(sum(mod(cds(:),2))==1 .or. sum(mod(cds(:),2))==2)cycle ! Either all are zero, or all are one, so skip when sum is 1 or 2. 308 end if 309 else 310 if(downsampling(1)==0 .and. mod(cds(2)+cds(3),2)/=0)cycle 311 if(downsampling(2)==0 .and. mod(cds(1)+cds(3),2)/=0)cycle 312 if(downsampling(3)==0 .and. mod(cds(1)+cds(2),2)/=0)cycle 313 end if 314 end if 315 end if 316 317 do ikshft=1,nshiftk 318 319 ! Only the first shiftk is taken into account if downsampling 320 ! if(.false.)then 321 if(present(downsampling))then 322 if(.not.(downsampling(1)==1 .and. downsampling(2)==1 .and. downsampling(3)==1))then 323 if(ikshft>1)cycle 324 end if 325 end if 326 327 ! Coordinates of the trial k point with respect to the k primitive lattice 328 k1(1)=ii+shiftk(1,ikshft) 329 k1(2)=jj+shiftk(2,ikshft) 330 k1(3)=kk+shiftk(3,ikshft) 331 ! Reduced coordinates of the trial k point 332 k2(:)=k1(1)*klatt(:,1)+k1(2)*klatt(:,2)+k1(3)*klatt(:,3) 333 ! Eliminate the point if outside [0,1[ 334 if(k2(1)<-tol10)cycle ; if(k2(1)>one-tol10)cycle 335 if(k2(2)<-tol10)cycle ; if(k2(2)>one-tol10)cycle 336 if(k2(3)<-tol10)cycle ; if(k2(3)>one-tol10)cycle 337 ! Wrap the trial values in the interval ]-1/2,1/2] . 338 call wrap2_pmhalf(k2(1),k1(1),shift) 339 call wrap2_pmhalf(k2(2),k1(2),shift) 340 call wrap2_pmhalf(k2(3),k1(3),shift) 341 spkpt(:,nn)=k1(:) 342 nn=nn+1 343 344 if(present(downsampling))then 345 found1(coord(2),coord(3))=1 346 found2(coord(1),coord(3))=1 347 found3(coord(1),coord(2))=1 348 end if 349 350 end do 351 end do 352 end do 353 end do 354 nkpt=nn-1 355 356 if(present(downsampling))then 357 ABI_DEALLOCATE(found1) 358 ABI_DEALLOCATE(found2) 359 ABI_DEALLOCATE(found3) 360 end if 361 362 363 if(nkpt/=nkptlatt*nshiftk)then 364 write(message, '(a,i8,a,a,a,i8,a)' )& 365 & 'The number of k points ',nkpt,' is not equal to',ch10,& 366 & 'nkptlatt*nshiftk which is',nkptlatt*nshiftk,'.' 367 !DEBUG 368 ! write(std_out,*)' smpbz : brav,iout,mkpt,nkpt,option=',brav,iout,mkpt,nkpt,option 369 ! write(std_out,*)' smpbz : kptrlatt(:,:)=',kptrlatt(:,:) 370 ! write(std_out,*)' smpbz : nshiftk=',nshiftk 371 ! write(std_out,*)' smpbz : shiftk(:,:)=',shiftk(:,:) 372 ! write(std_out,*)' smpbz : downsampling(:)=',downsampling(:) 373 ! write(std_out,*)message 374 ! stop 375 !ENDDEBUG 376 MSG_BUG(message) 377 end if 378 379 else if(brav==2)then 380 381 ! Face-Centered Lattice 382 if (prtvol > 0) call wrtout(std_out,' Face-Centered Lattice Grid ','COLL') 383 if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2) then 384 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 385 & 'The value of mkpt is not large enough. It should be',ch10,& 386 & 'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,',',ch10,& 387 & 'Action: set mkpt to that value in the main routine,',ch10,& 388 & 'and recompile the code.' 389 MSG_BUG(message) 390 end if 391 nn=1 392 if (ngkpt(1)/=ngkpt(2).or.ngkpt(1)/=ngkpt(3)) then 393 write(message, '(4a,3(a,i6,a),a)' )& 394 & 'For face-centered lattices, the numbers ngqpt(1:3)',ch10,& 395 & 'must be equal, while they are :',ch10,& 396 & 'ngqpt(1) = ',ngkpt(1),ch10,& 397 & 'ngqpt(2) = ',ngkpt(2),ch10,& 398 & 'ngqpt(3) = ',ngkpt(3),ch10,& 399 & 'Action: modify ngqpt(1:3) in the input file.' 400 MSG_BUG(message) 401 end if 402 if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2)) then 403 write(message, '(4a,3(a,i6,a),a)' )& 404 & 'For face-centered lattices, the numbers ngqpt(1:3)*nshiftk',ch10,& 405 & 'must be even, while they are :',ch10,& 406 & 'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,& 407 & 'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,& 408 & 'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,& 409 & 'Action: modify ngqpt(1:3)*nshiftk in the input file.' 410 MSG_ERROR(message) 411 end if 412 if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then 413 spkpt(1,1)=0.0_dp 414 spkpt(2,1)=0.0_dp 415 spkpt(3,1)=0.0_dp 416 nkpt=1 417 else 418 do kk=1,ngkpt(3) 419 do jj=1,ngkpt(2) 420 do ii=1,ngkpt(1) 421 do ikshft=1,nshiftk 422 k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1) 423 k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2) 424 k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3) 425 ! Wrap the trial values in the interval ]-1/2,1/2] . 426 call wrap2_pmhalf(k1(1),k2(1),shift) 427 call wrap2_pmhalf(k1(2),k2(2),shift) 428 call wrap2_pmhalf(k1(3),k2(3),shift) 429 ! Test whether it is inside the FCC BZ. 430 ktest(1)=2*k2(1)-1.0d-10 431 ktest(2)=2*k2(2)-2.0d-10 432 ktest(3)=2*k2(3)-5.0d-10 433 if (abs(ktest(1))+abs(ktest(2))+abs(ktest(3))<1.5_dp) then 434 kcar(1)=ktest(1)+1.0d-10 435 kcar(2)=ktest(2)+2.0d-10 436 kcar(3)=ktest(3)+5.0d-10 437 spkpt(1,nn)=0.5_dp*kcar(2)+0.5_dp*kcar(3) 438 spkpt(2,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(3) 439 spkpt(3,nn)=0.5_dp*kcar(1)+0.5_dp*kcar(2) 440 nn=nn+1 441 end if 442 end do 443 end do 444 end do 445 end do 446 nkpt=nn-1 447 if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/2)then 448 write(message, '(a,i8,a,a,a,i8,a)' )& 449 & 'The number of k points ',nkpt,' is not equal to',ch10,& 450 & '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2 which is',& 451 & (ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/2,'.' 452 MSG_BUG(message) 453 end if 454 end if 455 456 else if(brav==3)then 457 458 ! Body-Centered Lattice (not mandatory cubic !) 459 if (prtvol > 0) call wrtout(std_out,' Body-Centered Lattice Grid ','COLL') 460 if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk/4) then 461 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 462 & 'The value of mkpt is not large enough. It should be',ch10,& 463 & 'at least',(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,',',ch10,& 464 & 'Action: set mkpt to that value in the main routine,',ch10,& 465 & 'and recompile the code.' 466 MSG_BUG(message) 467 end if 468 nn=1 469 if ((ngkpt(1)*nshiftk)/=(((ngkpt(1)*nshiftk)/2)*2) .or.& 470 & (ngkpt(2)*nshiftk)/=(((ngkpt(2)*nshiftk)/2)*2) .or.& 471 & (ngkpt(3)*nshiftk)/=(((ngkpt(3)*nshiftk)/2)*2) ) then 472 write(message, '(4a,3(a,i6,a),a)' )& 473 & 'For body-centered lattices, the numbers ngqpt(1:3)',ch10,& 474 & 'must be even, while they are :',ch10,& 475 & 'ngqpt(1)*nshiftk = ',ngkpt(1)*nshiftk,ch10,& 476 & 'ngqpt(2)*nshiftk = ',ngkpt(2)*nshiftk,ch10,& 477 & 'ngqpt(3)*nshiftk = ',ngkpt(3)*nshiftk,ch10,& 478 & 'Action: modify ngqpt(1:3) in the input file.' 479 MSG_ERROR(message) 480 end if 481 if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then 482 spkpt(1,1)=0.0_dp 483 spkpt(2,1)=0.0_dp 484 spkpt(3,1)=0.0_dp 485 nkpt=1 486 else 487 do kk=1,ngkpt(3) 488 do jj=1,ngkpt(2) 489 do ii=1,ngkpt(1) 490 do ikshft=1,nshiftk 491 k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1) 492 k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2) 493 k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3) 494 ! Wrap the trial values in the interval ]-1/2,1/2] . 495 call wrap2_pmhalf(k1(1),k2(1),shift) 496 call wrap2_pmhalf(k1(2),k2(2),shift) 497 call wrap2_pmhalf(k1(3),k2(3),shift) 498 ! Test whether it is inside the BCC BZ. 499 ktest(1)=2*k2(1)-1.0d-10 500 ktest(2)=2*k2(2)-2.0d-10 501 ktest(3)=2*k2(3)-5.0d-10 502 if (abs(ktest(1))+abs(ktest(2))<1._dp) then 503 if (abs(ktest(1))+abs(ktest(3))<1._dp) then 504 if (abs(ktest(2))+abs(ktest(3))<1._dp) then 505 kcar(1)=ktest(1)+1.0d-10 506 kcar(2)=ktest(2)+2.0d-10 507 kcar(3)=ktest(3)+5.0d-10 508 spkpt(1,nn)=-0.5*kcar(1)+0.5*kcar(2)+0.5*kcar(3) 509 spkpt(2,nn)=0.5*kcar(1)-0.5*kcar(2)+0.5*kcar(3) 510 spkpt(3,nn)=0.5*kcar(1)+0.5*kcar(2)-0.5*kcar(3) 511 nn=nn+1 512 end if 513 end if 514 end if 515 end do 516 end do 517 end do 518 end do 519 nkpt=nn-1 520 if(nkpt==0)then 521 write(message, '(3a)' )& 522 & 'BCC lattice, input ngqpt=0, so no kpt is generated.',ch10,& 523 & 'Action: modify ngqpt(1:3) in the input file.' 524 MSG_ERROR(message) 525 end if 526 if(nkpt/=(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4)then 527 write(message, '(a,i8,a,a,a,i8,a)' )& 528 & 'The number of k points ',nkpt,' is not equal to',ch10,& 529 & '(ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4 which is',& 530 & (ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)/4,'.' 531 MSG_BUG(message) 532 end if 533 end if 534 535 else if(brav==4)then 536 537 ! Hexagonal Lattice (D6h) 538 if (prtvol > 0) call wrtout(std_out,' Hexagonal Lattice Grid ','COLL') 539 if (mkpt<ngkpt(1)*ngkpt(2)*ngkpt(3)) then 540 write(message, '(a,a,a,i8,a,a,a,a,a)' )& 541 & 'The value of mkpt is not large enough. It should be',ch10,& 542 & 'at least',ngkpt(1)*ngkpt(2)*ngkpt(3),',',ch10,& 543 & 'Action: set mkpt to that value in the main routine,',ch10,& 544 & 'and recompile the code.' 545 MSG_BUG(message) 546 end if 547 nn=1 548 if (ngkpt(1)/=ngkpt(2)) then 549 write(message, '(4a,2(a,i6,a),a)' )& 550 & 'For hexagonal lattices, the numbers ngqpt(1:2)',ch10,& 551 & 'must be equal, while they are :',ch10,& 552 & 'ngqpt(1) = ',ngkpt(1),ch10,& 553 & 'ngqpt(2) = ',ngkpt(2),ch10,& 554 & 'Action: modify ngqpt(1:3) in the input file.' 555 MSG_ERROR(message) 556 end if 557 if (ngkpt(1)==0.or.ngkpt(2)==0.or.ngkpt(3)==0) then 558 write(message, '(3a)' )& 559 & 'For hexagonal lattices, ngqpt(1:3)=0 is not permitted',ch10,& 560 & 'Action: modify ngqpt(1:3) in the input file.' 561 MSG_ERROR(message) 562 else 563 do kk=1,ngkpt(3) 564 do jj=1,ngkpt(2) 565 do ii=1,ngkpt(1) 566 do ikshft=1,nshiftk 567 k1(1)=(ii-1+shiftk(1,ikshft))/ngkpt(1) 568 k1(2)=(jj-1+shiftk(2,ikshft))/ngkpt(2) 569 k1(3)=(kk-1+shiftk(3,ikshft))/ngkpt(3) 570 ! Wrap the trial values in the interval ]-1/2,1/2] . 571 call wrap2_pmhalf(k1(1),k2(1),shift) 572 call wrap2_pmhalf(k1(2),k2(2),shift) 573 call wrap2_pmhalf(k1(3),k2(3),shift) 574 spkpt(:,nn)=k2(:) 575 nn=nn+1 576 end do 577 end do 578 end do 579 end do 580 nkpt=nn-1 581 if(nkpt/=ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk)then 582 write(message, '(a,i8,a,a,a,i8,a)' )& 583 & 'The number of k points ',nkpt,' is not equal to',ch10,& 584 & 'ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk which is',& 585 & ngkpt(1)*ngkpt(2)*ngkpt(3)*nshiftk,'.' 586 MSG_BUG(message) 587 end if 588 end if 589 590 else 591 592 write(message, '(a,i6,a,a,a)' )& 593 & 'The calling routine asks brav=',brav,'.',ch10,& 594 & 'but only brav=1 or -1,2,3 or 4 are allowed.' 595 MSG_BUG(message) 596 end if 597 598 if (option/=0) then 599 ! Put the Gamma point first 600 if(nkpt>1)then 601 do ii=1,nkpt 602 if(sum(abs(spkpt(:,ii)))<tol8)then 603 spkpt(:,ii)=spkpt(:,1) 604 spkpt(:,1)=zero 605 exit 606 end if 607 end do 608 end if 609 610 write(message,'(a,i8)')' Grid q points : ',nkpt 611 call wrtout(iout,message,'COLL') 612 nkpout=nkpt 613 if(nkpt>80)then 614 call wrtout(iout,' greater than 80, so only write 20 of them ','COLL') 615 nkpout=20 616 end if 617 do ii=1,nkpout 618 write(message, '(1x,i2,a2,3es16.8)' )ii,') ',spkpt(1,ii),spkpt(2,ii),spkpt(3,ii) 619 call wrtout(iout,message,'COLL') 620 end do 621 end if 622 623 end subroutine smpbz