Line data Source code
1 : !!****m* ABINIT/m_inkpts
2 : !! NAME
3 : !! m_inkpts
4 : !!
5 : !! FUNCTION
6 : !! Routines to initialize k-point and q-point sampling from input file.
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_inkpts
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_xmpi
28 : use m_nctk
29 : use netcdf
30 : use m_hdr
31 :
32 : use m_time, only : timab
33 : use m_fstrings, only : sjoin, itoa
34 : use m_numeric_tools, only : isdiagmat
35 : use m_geometry, only : metric
36 : use m_symfind, only : symfind, symlatt
37 : use m_cgtools, only : set_istwfk
38 : use m_parser, only : intagm
39 : use m_kpts, only : getkgrid, testkgrid, mknormpath
40 :
41 : implicit none
42 :
43 : private
44 : !!***
45 :
46 : public :: inkpts ! Initialize k points (list of k points, weights, storage) for one particular dataset.
47 : public :: inqpt ! Initialize the q point for one particular dataset, characterized by jdtset.
48 : !!***
49 :
50 : contains
51 : !!***
52 :
53 : !!****f* m_inkpts/inkpts
54 : !! NAME
55 : !! inkpts
56 : !!
57 : !! FUNCTION
58 : !! Initialize k points (list of k points, weights, storage)
59 : !! for one particular dataset, characterized by jdtset.
60 : !! Note that nkpt (and nkpthf) can be computed by calling this routine with input value of nkpt=0, provided kptopt /= 0.
61 : !!
62 : !! INPUTS
63 : !! bravais(11): bravais(1)=iholohedry
64 : !! bravais(2)=center
65 : !! bravais(3:11)=coordinates of rprim in the axes of the conventional bravais lattice (*2 if center/=0)
66 : !! chksymbreak= if 1, will check whether the k point grid is symmetric, and stop if not.
67 : !! [impose_istwf_1]= (optional argument):
68 : !! 0: no restriction on istwfk
69 : !! 1: impose istwfk=1 for all k points
70 : !! 2: impose istwfk=1 for all k points non equal to zero
71 : !! iout=unit number for echoed output
72 : !! iscf= <= 0 => non-SCF, >0 => SCF.
73 : !! jdtset=number of the dataset looked for
74 : !! lenstr=actual length of the string
75 : !! kptopt=option for the generation of k points
76 : !! msym=default maximal number of symmetries
77 : !! getkerange_filepath= Path of KERANGE.nc file used to initialize k-point sampling if kptopt == 0 and string != ABI_NOFILE
78 : !! nqpt=number of q points (0 or 1)
79 : !! nsym=number of symetries
80 : !! occopt=option for occupation numbers
81 : !! qptn(3)=reduced coordinates of eventual q point shift (already normalized).
82 : !! response=0 if GS case, =1 if RF case.
83 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
84 : !! string*(*)=character string containing all the input data. Initialized previously in instrng.
85 : !! symafm(nsym)=(anti)ferromagnetic part of symmetry operations
86 : !! symrel(3,3,nsym)=symmetry operations in real space in terms of primitive translations
87 : !! vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
88 : !! comm= MPI communicator
89 : !!
90 : !! OUTPUT
91 : !! fockdownsampling(3)=echo of input variable fockdownsampling(3)
92 : !! kptnrm=normalisation of k points
93 : !! kptrlatt_orig(3,3)=Original value of kptrlatt as specified in the input file (if kptopt/=0)
94 : !! kptrlatt(3,3)=k-point lattice specification (if kptopt/=0)
95 : !! kptrlen=length of the smallest real space supercell vector
96 : !! nshiftk_orig=Original number of k-point shifts (0 if not read)
97 : !! nshiftk=actual number of k-point shifts in shiftk (if kptopt/=0)
98 : !! shiftk(3,MAX_NSHIFTK)=shift vectors for k point generation (if kptopt/=0)
99 : !! If nkpt/=0 the following arrays are also output:
100 : !! istwfk(nkpt)=option parameters that describes the storage of wfs
101 : !! kpt(3,nkpt)=reduced coordinates of k points.
102 : !! kpthf(3,nkpthf)=reduced coordinates of k points for Fock operator.
103 : !! wtk(nkpt)=weight assigned to each k point.
104 : !! ngkpt(3)=Number of divisions along the three reduced directions
105 : !! (0 signals that this variable has not been used.
106 : !! shiftk_orig(3,MAX_NSHIFTK)=Original shifts read from the input file
107 : !! (0 signals that this variable has not been read).
108 : !!
109 : !! SIDE EFFECTS
110 : !! Input/output:
111 : !! nkpt=number of k points
112 : !! if non-zero at input, is only an input variable
113 : !! if zero at input, its actual value will be computed
114 : !! nkpthf=number of k points for Fock operator, computed if nkpt=0 at input
115 : !!
116 : !! NOTES
117 : !! Warning: this routine can be called with nkpt=0 (in which
118 : !! case it returns the true value of nkpt), which can lead
119 : !! to strange bugs in the debugging procedure, if one tries to print wtk or istwfk, in this case!
120 : !!
121 : !! SOURCE
122 :
123 13923 : subroutine inkpts(bravais,chksymbreak,fockdownsampling,iout,iscf,istwfk,jdtset,&
124 13923 : & kpt,kpthf,kptopt,kptnrm,kptrlatt_orig,kptrlatt,kptrlen,lenstr,msym, getkerange_filepath, &
125 : & nkpt,nkpthf,nqpt,ngkpt,nshiftk,nshiftk_orig,shiftk_orig,nsym,&
126 13923 : & occopt,qptn,response,rprimd,shiftk,string,symafm,symrel,vacuum,wtk,comm,&
127 : & impose_istwf_1) ! Optional argument
128 :
129 : !Arguments ------------------------------------
130 : !scalars
131 : integer,intent(in) :: chksymbreak,iout,iscf,jdtset,kptopt,lenstr,msym,nqpt,nsym,occopt
132 : integer,intent(in) :: response, comm
133 : integer,intent(in),optional :: impose_istwf_1
134 : integer,intent(inout) :: nkpt,nkpthf
135 : integer,intent(out) :: nshiftk,nshiftk_orig
136 : integer,intent(out) :: fockdownsampling(3)
137 : real(dp),intent(out) :: kptnrm,kptrlen
138 : character(len=*),intent(in) :: string
139 : character(len=*),intent(in) :: getkerange_filepath
140 : !arrays
141 : integer,intent(in) :: bravais(11),symafm(msym),symrel(3,3,msym),vacuum(3)
142 : integer,intent(out) :: istwfk(nkpt),kptrlatt(3,3),kptrlatt_orig(3,3),ngkpt(3)
143 : real(dp),intent(in) :: rprimd(3,3),qptn(3)
144 : real(dp),intent(out) :: kpt(3,nkpt),kpthf(3,nkpthf),shiftk(3,MAX_NSHIFTK),wtk(nkpt),shiftk_orig(3,MAX_NSHIFTK)
145 :
146 : !Local variables-------------------------------
147 : !scalars
148 : integer,parameter :: master = 0
149 : integer :: dkpt,ii,ikpt,jkpt,marr,ndiv_small,nkpt_computed,my_rank,nprocs
150 : integer :: nsegment,prtkpt,tread,tread_kptrlatt,tread_ngkpt, ncid, fform, ierr
151 : logical :: use_kerange
152 : real(dp) :: fraction,norm,ucvol,wtksum
153 : character(len=500) :: msg
154 13923 : type(hdr_type) :: hdr
155 : !arrays
156 13923 : integer,allocatable :: ndivk(:),intarr(:), krange2ibz(:)
157 : real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),rmet(3,3),tsec(2)
158 13923 : real(dp),allocatable :: kptbounds(:,:),dprarr(:)
159 :
160 : ! *************************************************************************
161 :
162 13923 : call timab(192,1,tsec)
163 :
164 13923 : my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm)
165 :
166 : ! Compute the maximum size of arrays intarr and dprarr
167 13923 : marr = max(3*nkpt,3*MAX_NSHIFTK)
168 41769 : ABI_MALLOC(intarr,(marr))
169 41769 : ABI_MALLOC(dprarr,(marr))
170 :
171 : ! Use zero to signal that these values have not been read.
172 13923 : ngkpt = 0
173 13923 : shiftk_orig = zero
174 13923 : kptrlatt_orig = 0; kptrlatt = 0
175 13923 : nshiftk_orig = 1; nshiftk = 1
176 13923 : use_kerange = .False.
177 :
178 : !fockdownsampling(:)=1
179 : !kptnrm = one
180 : !kpthf = zero
181 :
182 : ! MG: FIXME These values should be initialized because they are intent(out)
183 : ! but several tests fails. So we keep this bug to avoid problems somewhere else
184 : ! The initialization of the kpoints should be rewritten in a cleaner way
185 : ! without all these side effects!
186 : !shiftk = zero
187 : ! Initializing these three variables is OK but we keep the bug to preserve the old behavior
188 : !wtk = one
189 : !kpt = zero
190 : !istwfk = 1
191 :
192 : ! Initialize kptrlen
193 13923 : kptrlen=30.0_dp
194 13923 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptrlen',tread,'DPR')
195 13923 : if(tread==1)kptrlen=dprarr(1)
196 :
197 : ! Initialize kpt, kptnrm and wtk according to kptopt.
198 13923 : if (kptopt == 0 .and. getkerange_filepath == ABI_NOFILE) then
199 : ! For kptopt==0, one must have nkpt defined.
200 48452 : kpt(:,:)=zero
201 6356 : call intagm(dprarr,intarr,jdtset,marr,3*nkpt,string(1:lenstr),'kpt',tread,'DPR')
202 10900 : if(tread==1) kpt(:,:)=reshape( dprarr(1:3*nkpt), [3,nkpt])
203 :
204 6356 : kptnrm=one
205 6356 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR')
206 6356 : if(tread==1) kptnrm=dprarr(1)
207 :
208 : ! Only read wtk when iscf >0 or iscf=-1 or iscf=-3 or (iscf=-2 and response=1)
209 : ! (this last option is for Zach Levine)
210 : ! Normalize the k-point weights when occopt/=2
211 : ! Check that k point weights add to 1 when occopt==2
212 6356 : if (iscf>0.or.iscf==-1.or.iscf==-3.or.(iscf==-2.and.response==1)) then
213 16162 : wtk = one
214 6248 : call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'wtk',tread,'DPR')
215 11524 : if(tread==1) wtk(1:nkpt)=dprarr(1:nkpt)
216 :
217 16162 : wtksum=sum(wtk(:))
218 6248 : write(msg,'(a,i0,a,f12.6)')' inkpts: Sum of ',nkpt,' k point weights is',wtksum
219 6248 : call wrtout(std_out,msg)
220 :
221 6248 : if (wtksum < tol6) then
222 : write(msg, '(3a)' )&
223 0 : 'This sum is too close to zero. ',ch10,&
224 0 : 'Action: correct the array wtk in the input file.'
225 0 : ABI_ERROR(msg)
226 : end if
227 6248 : if (abs(wtksum - one) > tol6) then
228 344 : if (occopt==2) then
229 : write(msg, '(a,1p,e18.8,a,a,a)' )&
230 0 : 'wtksum= ',wtksum,' /= 1.0 means wts do not add to 1 , while occopt=2.',ch10,&
231 0 : 'Action: correct the array wtk in input file.'
232 0 : ABI_ERROR(msg)
233 : else
234 344 : write(msg,'(a,i0,a)')' With present occopt= ',occopt,', renormalize it to one'
235 344 : call wrtout(std_out,msg)
236 344 : norm=one/wtksum
237 1334 : wtk(1:nkpt)=wtk(1:nkpt)*norm
238 : end if
239 : end if
240 : end if
241 :
242 7567 : else if (kptopt == 0 .and. getkerange_filepath /= ABI_NOFILE) then
243 : ! Initialize k-points from kerange_path file.
244 4 : use_kerange = .True.
245 12 : ABI_MALLOC(krange2ibz, (nkpt))
246 4 : if (my_rank == master) then
247 4 : NCF_CHECK(nctk_open_read(ncid, getkerange_filepath, xmpi_comm_self))
248 4 : call hdr%ncread(ncid, fform)
249 4 : ABI_CHECK(fform == fform_from_ext("KERANGE.nc"), sjoin("Error while reading:", getkerange_filepath, ", fform:", itoa(fform)))
250 : ! TODO Add code for consistency check
251 : !kptopt, nsym, occopt
252 : !ABI_CHECK(nkpt == hdr%nkpt, "nkpt from kerange != nkpt")
253 4 : NCF_CHECK(nf90_get_var(ncid, nctk_idname(ncid, "krange2ibz"), krange2ibz))
254 4 : NCF_CHECK(nf90_close(ncid))
255 : end if
256 :
257 4 : call xmpi_bcast(krange2ibz, master, comm, ierr)
258 4 : call hdr%bcast(master, my_rank, comm)
259 : ! Hdr contains the kpts in the IBZ. Extract k-points in the pockets via krange2ibz.
260 4 : nshiftk = hdr%nshiftk; nshiftk_orig = hdr%nshiftk_orig
261 38 : istwfk = hdr%istwfk(krange2ibz(:))
262 100 : kptrlatt = hdr%kptrlatt; kptrlatt_orig = hdr%kptrlatt_orig
263 4 : ABI_CHECK(isdiagmat(hdr%kptrlatt), "kptrlatt is not diagonal!")
264 4 : ngkpt(1) = hdr%kptrlatt(1, 1); ngkpt(2) = hdr%kptrlatt(2, 2); ngkpt(3) = hdr%kptrlatt(3, 3)
265 140 : kpt = hdr%kptns(:, krange2ibz(:)) !; kpthf(3,nkpthf)
266 36 : shiftk(:,1:nshiftk) = hdr%shiftk; shiftk_orig(:, 1:nshiftk_orig) = hdr%shiftk_orig
267 38 : wtk = hdr%wtk(krange2ibz(:))
268 4 : call hdr%free()
269 4 : kptnrm = one
270 4 : ABI_FREE(krange2ibz)
271 :
272 7563 : else if (kptopt < 0) then
273 : ! Band structure calculation
274 284 : nsegment=abs(kptopt)
275 :
276 284 : if (iscf /= -2)then
277 : write(msg, '(3a,i0,3a)' ) &
278 0 : 'For a negative kptopt, iscf must be -2,',ch10,&
279 0 : 'while it is found to be ',iscf,'.',ch10,&
280 0 : 'Action: change the value of iscf in your input file, or change kptopt.'
281 0 : ABI_ERROR(msg)
282 : end if
283 :
284 284 : if(marr<3*nsegment+3)then
285 0 : marr=3*nsegment+3
286 0 : ABI_FREE(dprarr)
287 0 : ABI_FREE(intarr)
288 0 : ABI_MALLOC(dprarr,(marr))
289 0 : ABI_MALLOC(intarr,(marr))
290 : end if
291 :
292 852 : ABI_MALLOC(kptbounds,(3,nsegment+1))
293 852 : ABI_MALLOC(ndivk,(nsegment))
294 :
295 284 : call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR')
296 :
297 284 : if(tread==1)then
298 852 : kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1])
299 : else
300 : write(msg,'(5a)') &
301 0 : 'When kptopt is negative, kptbounds must be initialized ',ch10,&
302 0 : 'in the input file, which is not the case.',ch10,&
303 0 : 'Action: initialize kptbounds in your input file, or change kptopt.'
304 0 : ABI_ERROR(msg)
305 : end if
306 :
307 284 : call intagm(dprarr,intarr,jdtset,marr,nsegment,string(1:lenstr),'ndivk',tread,'INT')
308 284 : if(tread==1)then
309 450 : ndivk(1:nsegment)=intarr(1:nsegment)
310 : ! The 1 stand for the first point
311 450 : nkpt_computed=1+sum(ndivk(1:nsegment))
312 :
313 : ! ndivk and ndivsm are mutually exclusive
314 214 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT')
315 214 : if (tread == 1) then
316 0 : ABI_ERROR("ndivk and ndivsm are mutually exclusive. Choose only one variable")
317 : end if
318 :
319 : else
320 : ! Calculate ndivk such as the path is normalized
321 : ! Note that if both ndivk and ndivsm are defined in in input file, only ndivk is used !
322 70 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'ndivsm',tread,'INT')
323 70 : if(tread==1)then
324 70 : ndiv_small=intarr(1)
325 70 : call metric(gmet,gprimd,std_out,rmet,rprimd,ucvol)
326 70 : call mknormpath(nsegment+1,kptbounds,gmet,ndiv_small,ndivk,nkpt_computed)
327 : else
328 : write(msg,'(5a)') &
329 0 : 'When kptopt is negative, ndivsm or ndivk must be initialized ',ch10,&
330 0 : 'in the input file, which is not the case.',ch10,&
331 0 : 'Action: initialize ndivsm or ndivk in your input file, or change kptopt.'
332 0 : ABI_ERROR(msg)
333 : end if
334 : end if
335 :
336 : ! Check that the argument nkpt is coherent with nkpt_computed
337 284 : if (nkpt/=0 .and. nkpt /= nkpt_computed) then
338 : write(msg, '(a,i0,5a,i0,7a)' ) &
339 0 : 'The argument nkpt = ',nkpt,', does not match',ch10,&
340 0 : 'the number of k points generated by kptopt, ndivk, kptbounds,',ch10,&
341 0 : 'and the eventual symmetries, that is, nkpt= ',nkpt_computed,'.',ch10,&
342 0 : 'However, note that it might due to the user,',ch10,&
343 0 : 'if nkpt is explicitely defined in the input file.',ch10,&
344 0 : 'In this case, please check your input file.'
345 0 : ABI_ERROR(msg)
346 : end if
347 :
348 284 : if (nkpt/=0) then
349 : ! The array kpt has the right dimension and we can generate the k-path
350 143 : call intagm(dprarr,intarr,jdtset,marr,3*nsegment+3,string(1:lenstr),'kptbounds',tread,'DPR')
351 143 : if(tread==1)then
352 429 : kptbounds(:,:)=reshape( dprarr(1:3*nsegment+3), [3,nsegment+1])
353 : else
354 : write(msg, '(5a)') &
355 0 : 'When kptopt is negative, kptbounds must be initialized ',ch10,&
356 0 : 'in the input file, which is not the case.',ch10,&
357 0 : 'Action: initialize kptbounds in your input file, or change kptopt.'
358 0 : ABI_ERROR(msg)
359 : end if
360 :
361 : ! First k point
362 572 : jkpt=1
363 572 : kpt(:,1)=kptbounds(:,1)
364 523 : do ii=1,nsegment
365 380 : dkpt=ndivk(ii)
366 3428 : do ikpt=1,dkpt
367 3048 : fraction=dble(ikpt)/dble(dkpt)
368 12572 : kpt(:,ikpt+jkpt)=fraction *kptbounds(:,ii+1)+(one-fraction)*kptbounds(:,ii)
369 : end do
370 523 : jkpt=jkpt+dkpt
371 : end do
372 :
373 : end if
374 :
375 284 : kptnrm=one
376 284 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'kptnrm',tread,'DPR')
377 284 : if(tread==1) kptnrm=dprarr(1)
378 :
379 284 : ABI_FREE(kptbounds)
380 284 : ABI_FREE(ndivk)
381 :
382 7279 : else if (kptopt>=1 .and. kptopt<=4) then
383 : ! Read ngkpt
384 7279 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngkpt',tread_ngkpt,'INT')
385 7279 : if(tread_ngkpt==1)then
386 26080 : ngkpt(1:3)=intarr(1:3)
387 26080 : do ii=1,3
388 26080 : if(ngkpt(ii)<1)then
389 : write(msg,'(a,i0,3a,i0,3a)') &
390 0 : 'The input variable ngkpt(',ii,') must be strictly positive,',ch10,&
391 0 : 'while it is found to be ',ngkpt(ii),'.',ch10,&
392 0 : 'Action: change it in your input file, or change kptopt.'
393 0 : ABI_ERROR(msg)
394 : end if
395 : end do
396 : end if
397 :
398 7279 : call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'kptrlatt',tread_kptrlatt,'INT')
399 7279 : if(tread_kptrlatt==1) kptrlatt(:,:)=reshape(intarr(1:9), [3,3])
400 :
401 7279 : if(tread_ngkpt==1 .and. tread_kptrlatt==1)then
402 : write(msg, '(5a)') &
403 0 : 'The input variables ngkpt and kptrlatt cannot both ',ch10,&
404 0 : 'be defined in the input file.',ch10,&
405 0 : 'Action: change one of ngkpt or kptrlatt in your input file.'
406 0 : ABI_ERROR(msg)
407 7279 : else if(tread_ngkpt==1)then
408 6520 : kptrlatt(:,:)=0
409 6520 : kptrlatt(1,1)=ngkpt(1)
410 6520 : kptrlatt(2,2)=ngkpt(2)
411 6520 : kptrlatt(3,3)=ngkpt(3)
412 : ! Save kptrlatt for reference.
413 6520 : kptrlatt_orig = kptrlatt
414 : end if
415 :
416 7279 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftk',tread,'INT')
417 7279 : if(tread==1)nshiftk=intarr(1)
418 :
419 7279 : if (nshiftk < 1 .or. nshiftk > MAX_NSHIFTK )then
420 : write(msg, '(a,i0,2a,i0,3a)' )&
421 0 : 'The only allowed values of nshiftk are between 1 and ',MAX_NSHIFTK,ch10,&
422 0 : 'while it is found to be',nshiftk,'.',ch10,&
423 0 : 'Action: change the value of nshiftk in your input file, or change kptopt.'
424 0 : ABI_ERROR(msg)
425 : end if
426 :
427 7279 : call intagm(dprarr,intarr,jdtset,marr,3*nshiftk,string(1:lenstr),'shiftk',tread,'DPR')
428 7279 : if(tread==1)then
429 19416 : shiftk(:,1:nshiftk)=reshape( dprarr(1:3*nshiftk), [3,nshiftk])
430 : ! Save input shifts as they will be changes in getkgrid.
431 6472 : nshiftk_orig = nshiftk
432 65992 : shiftk_orig(:,1:nshiftk) = shiftk(:,1:nshiftk)
433 : else
434 807 : if(nshiftk/=1)then
435 : write(msg, '(3a,i0,2a)' )&
436 0 : 'When nshiftk is not equal to 1, shiftk must be defined in the input file.',ch10,&
437 0 : 'However, shiftk is not defined, while nshiftk=',nshiftk,ch10,&
438 0 : 'Action: change the value of nshiftk in your input file, or define shiftk.'
439 0 : ABI_ERROR(msg)
440 : end if
441 : ! Default values used in indefo
442 807 : nshiftk_orig = 1
443 4035 : shiftk_orig(:,1:nshiftk) = half
444 : end if
445 :
446 7279 : prtkpt=0
447 7279 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'prtkpt',tread,'INT')
448 7279 : if(tread==1)prtkpt=intarr(1)
449 :
450 94627 : if(sum(abs(kptrlatt(:,:)))==0)then
451 197 : kptrlatt_orig = 0
452 788 : do ii=1,3
453 788 : kptrlatt_orig(ii,ii) = ngkpt(ii)
454 : end do
455 : ! The parameters of the k lattice are not known, compute kptrlatt, nshiftk, shiftk.
456 197 : call testkgrid(bravais,iout,kptrlatt,kptrlen, msym,nshiftk,nsym,prtkpt,rprimd,shiftk,symafm,symrel,vacuum)
457 : end if
458 :
459 : ! TODO: Avoid call to getkgrid if eph
460 : !eph_task = -1
461 : !call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'eph_task',tread,'INT')
462 : !if(tread==1) eph_task=intarr(1)
463 :
464 29064 : fockdownsampling(:)=1
465 7266 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'fockdownsampling',tread,'INT')
466 7782 : if(tread==1)fockdownsampling=intarr(1:3)
467 :
468 : call getkgrid(chksymbreak,0,iscf,kpt,kptopt,kptrlatt,kptrlen,&
469 : msym,nkpt,nkpt_computed,nshiftk,nsym,rprimd,&
470 7266 : shiftk,symafm,symrel,vacuum,wtk,nkpthf=nkpthf,kpthf=kpthf,downsampling=fockdownsampling)
471 :
472 7266 : kptnrm=one
473 :
474 : else
475 : write(msg,'(3a,i0,3a)' ) &
476 0 : 'The only values of kptopt allowed are smaller than 4.',ch10,&
477 0 : 'The input value of kptopt is: ',kptopt,'.',ch10,&
478 0 : 'Action: change kptopt in your input file.'
479 0 : ABI_ERROR(msg)
480 : end if
481 :
482 13910 : if (kptnrm < tol10) then
483 : write(msg, '(5a)' )&
484 0 : 'The input variable kptnrm is lower than 1.0d-10,',ch10,&
485 0 : 'while it must be a positive, non-zero number. ',ch10,&
486 0 : 'Action: correct the kptnrm in the input file.'
487 0 : ABI_ERROR(msg)
488 : end if
489 :
490 : ! The k point number has been computed, and, if nkpt/=0, also the list of k points.
491 : ! Also nkpthf has been computed, and, if nkpt/=0, also the list kpthf.
492 : ! Now, determine istwfk, and eventually shift the k points by the value of qptn.
493 13910 : if (nkpt /= 0) then
494 138147 : istwfk(1:nkpt)=0
495 10450 : call intagm(dprarr,intarr,jdtset,marr,nkpt,string(1:lenstr),'istwfk',tread,'INT')
496 25427 : if(tread==1) istwfk(1:nkpt)=intarr(1:nkpt)
497 :
498 : ! Impose istwfk=1 for RF calculations or NSCF calculation with kpts from kerange.
499 75827 : if (response == 1 .or. use_kerange) istwfk(1:nkpt)=1
500 :
501 : ! Also impose istwfk=1 for spinor calculations
502 10450 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nspinor',tread,'INT')
503 12820 : if(tread/=0 .and. intarr(1)/=1)istwfk(1:nkpt)=1
504 :
505 10450 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'pawspnorb',tread,'INT')
506 10591 : if(tread/=0 .and. intarr(1)/=0)istwfk(1:nkpt)=1
507 :
508 10450 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'zora',tread,'INT')
509 10718 : if( tread/=0 .and. intarr(1)/=0 .and. intarr(1)/=1)istwfk(1:nkpt)=1
510 :
511 138147 : do ikpt=1,nkpt
512 127697 : if(istwfk(ikpt)==0)then
513 209202 : kpoint=kpt(:,ikpt)/kptnrm; if (nqpt/=0.and.response==0) kpoint=kpoint+qptn
514 48495 : istwfk(ikpt) = set_istwfk(kpoint)
515 : end if
516 138147 : if (present(impose_istwf_1)) then
517 113194 : if (impose_istwf_1==1) then
518 1452 : istwfk(ikpt)=1
519 232259 : else if (impose_istwf_1==2.and.any(kpt(:,ikpt)>tol10)) then
520 1175 : istwfk(ikpt)=1
521 : end if
522 : end if
523 : end do
524 : end if
525 :
526 : ! If nkpt was to be computed, transfer it from nkpt_computed
527 13910 : if (nkpt == 0) nkpt = nkpt_computed
528 :
529 13910 : ABI_FREE(intarr)
530 13910 : ABI_FREE(dprarr)
531 :
532 13910 : call timab(192,2,tsec)
533 :
534 13910 : end subroutine inkpts
535 : !!***
536 :
537 : !!****f* m_inkpts/inqpt
538 : !! NAME
539 : !! inqpt
540 : !!
541 : !! FUNCTION
542 : !! Initialize the q point for one particular dataset, characterized by jdtset.
543 : !!
544 : !! INPUTS
545 : !! chksymbreak= if 1, will check whether the k point grid is symmetric, and stop if not.
546 : !! iout=unit number for echoed output
547 : !! jdtset=number of the dataset looked for
548 : !! lenstr=actual length of the string
549 : !! msym=default maximal number of symmetries
550 : !! natom=number of atoms
551 : !! rprimd(3,3)=dimensional real space primitive translations (bohr)
552 : !! spinat(3,1:natom)=spin-magnetization of the atoms
553 : !! string*(*)=character string containing all the input data. Initialized previously in instrng.
554 : !! typat(natom)=type for each atom
555 : !! vacuum(3)=for each direction, 0 if no vacuum, 1 if vacuum
556 : !! xred(3,natom,nimage) =reduced coordinates of atoms
557 : !!
558 : !! OUTPUT
559 : !! qptn(3)=reduced coordinates of eventual q point (normalisation is already included)
560 : !! kptrlatt(3,3)=q-point lattice specification (if kptopt/=0)
561 : !! wtqc=weigth of the eventual current q point
562 : !!
563 : !! SOURCE
564 :
565 713 : subroutine inqpt(chksymbreak,iout,jdtset,lenstr,msym,natom,qptn,wtqc,rprimd,spinat,string,typat,vacuum,xred,qptrlatt)
566 :
567 : !Arguments ------------------------------------
568 : !scalars
569 : integer,intent(in) :: chksymbreak,iout,jdtset,lenstr,msym,natom
570 : real(dp),intent(inout) :: wtqc
571 : character(len=*),intent(in) :: string
572 : !arrays
573 : integer,intent(in) :: typat(natom),vacuum(3)
574 : real(dp),intent(out) :: qptn(3)
575 : integer,intent(inout) :: qptrlatt(3,3) !vz_i
576 : real(dp),intent(in) :: rprimd(3,3)
577 : real(dp),intent(in) :: spinat(3,natom)
578 : real(dp),intent(in) :: xred(3,natom)
579 :
580 : !Local variables-------------------------------
581 : !scalars
582 : integer :: ii,iqpt,iscf_fake,marr,nptsym,nqpt_max,nqpt_computed,nshiftq,nsym_new,qptopt
583 : integer :: tread,tread_q_sum,tread_qptrlatt,tread_ngqpt,use_inversion
584 : real(dp) :: qptnrm,qptrlen,tolsym,ucvol
585 : character(len=500) :: msg
586 : !arrays
587 : integer :: bravais(11), ngqpt(3)
588 713 : integer, allocatable :: symafm_new(:), ptsymrel(:,:,:),symrel_new(:,:,:), intarr(:)
589 : real(dp) :: gmet(3,3),gprimd(3,3),qpt(3),rmet(3,3),shiftq(3,MAX_NSHIFTK)
590 713 : real(dp),allocatable :: qpts(:,:),tnons_new(:,:),wtq(:), dprarr(:)
591 :
592 : ! *************************************************************************
593 :
594 : ! Compute the maximum size of arrays intarr and dprarr (nshiftq is MAX_NSHIFTK at maximum)
595 713 : marr=630
596 713 : ABI_MALLOC(intarr,(marr))
597 713 : ABI_MALLOC(dprarr,(marr))
598 713 : tread_q_sum=0
599 :
600 : ! Find the method to generate the q-points
601 713 : qptopt=0
602 713 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptopt',tread,'INT')
603 713 : tread_q_sum=tread_q_sum+tread
604 713 : if(tread==1)qptopt=intarr(1)
605 :
606 713 : if(qptopt==0)then
607 : ! Read qpt and qptnrm
608 615 : qpt=zero
609 615 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'qpt',tread,'DPR')
610 615 : tread_q_sum=tread_q_sum+tread
611 3019 : if(tread==1) qpt(1:3)=dprarr(1:3)
612 :
613 615 : qptnrm=1.0_dp
614 615 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'qptnrm',tread,'DPR')
615 615 : tread_q_sum=tread_q_sum+tread
616 :
617 615 : if(tread==1) qptnrm=dprarr(1)
618 4 : if(qptnrm<tol10)then
619 : write(msg, '(5a)' )&
620 0 : 'The input variable qptnrm is lower than 1.0d-10,',ch10,&
621 0 : 'while it must be a positive, non-zero number. ',ch10,&
622 0 : 'Action: correct the qptnrm in the input file.'
623 0 : ABI_ERROR(msg)
624 : end if
625 :
626 2460 : qptn(:)=qpt(:)/qptnrm
627 :
628 : ! DBSP: one could want ot define wtq in order to reproduce what is obtained
629 : ! with ngqpt but without having to do initialize the qgrid (extremly slow in case of large grid > 50x50x50
630 615 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'wtq',tread,'DPR')
631 615 : tread_q_sum=tread_q_sum+tread
632 615 : if(tread==1) wtqc=dprarr(1)
633 :
634 98 : else if (qptopt>=1 .and. qptopt<=4) then
635 98 : ngqpt(:)=0
636 98 : call intagm(dprarr,intarr,jdtset,marr,3,string(1:lenstr),'ngqpt',tread_ngqpt,'INT')
637 98 : tread_q_sum=tread_q_sum+tread_ngqpt
638 :
639 98 : if(tread_ngqpt==1)then
640 348 : ngqpt(1:3)=intarr(1:3)
641 348 : do ii=1,3
642 348 : if(ngqpt(ii)<1)then
643 : write(msg, '(a,i0,a,a,a,i0,a,a,a)' ) &
644 0 : 'The input variable ngqpt(',ii,') must be strictly positive,',ch10,&
645 0 : 'while it is found to be',ngqpt(ii),'.',ch10,&
646 0 : 'Action: change it in your input file, or change qptopt.'
647 0 : ABI_ERROR(msg)
648 : end if
649 : end do
650 : end if
651 :
652 98 : call intagm(dprarr,intarr,jdtset,marr,9,string(1:lenstr),'qptrlatt',tread_qptrlatt,'INT')
653 98 : if(tread_qptrlatt==1) qptrlatt(:,:)=reshape(intarr(1:9), (/3,3/) )
654 98 : tread_q_sum=tread_q_sum+tread_qptrlatt
655 :
656 98 : if(tread_ngqpt==1 .and. tread_qptrlatt==1)then
657 : write(msg, '(5a)' ) &
658 0 : 'The input variables ngqpt and qptrlatt cannot both ',ch10,&
659 0 : 'be defined in the input file.',ch10,&
660 0 : 'Action: change one of ngqpt or qptrlatt in your input file.'
661 0 : ABI_ERROR(msg)
662 98 : else if(tread_ngqpt==1)then
663 87 : qptrlatt(:,:)=0
664 87 : qptrlatt(1,1)=ngqpt(1)
665 87 : qptrlatt(2,2)=ngqpt(2)
666 87 : qptrlatt(3,3)=ngqpt(3)
667 : end if
668 :
669 98 : nshiftq=1
670 98 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'nshiftq',tread,'INT')
671 98 : tread_q_sum=tread_q_sum+tread
672 98 : if(tread==1)nshiftq=intarr(1)
673 :
674 98 : if (nshiftq<1 .or. nshiftq>MAX_NSHIFTK) then
675 : write(msg, '(a,i0,2a,i0,3a)' )&
676 0 : 'The only allowed values of nshiftq are between 1 and,',MAX_NSHIFTK,ch10,&
677 0 : 'while it is found to be',nshiftq,'.',ch10,&
678 0 : 'Action: change the value of nshiftq in your input file, or change qptopt.'
679 0 : ABI_ERROR(msg)
680 : end if
681 :
682 98 : shiftq=zero
683 98 : call intagm(dprarr,intarr,jdtset,marr,3*nshiftq,string(1:lenstr),'shiftq',tread,'DPR')
684 98 : tread_q_sum=tread_q_sum+tread
685 :
686 98 : if(tread==1)then
687 270 : shiftq(:,1:nshiftq)=reshape( dprarr(1:3*nshiftq), (/3,nshiftq/) )
688 : else
689 8 : if(nshiftq/=1)then
690 : write(msg, '(3a,i0,2a)' )&
691 0 : 'When nshiftq is not equal to 1, shiftq must be defined in the input file.',ch10,&
692 0 : 'However, shiftq is not defined, while nshiftq=',nshiftq,ch10,&
693 0 : 'Action: change the value of nshiftq in your input file, or define shiftq.'
694 0 : ABI_ERROR(msg)
695 : end if
696 : end if
697 :
698 : !write(std_out,'(a)')' m_inkpts%inqpt : before symlatt '
699 :
700 : ! Re-generate symmetry operations from the lattice and atomic coordinates
701 : ! This is a fundamental difference with respect to the k point generation.
702 98 : tolsym=tol8
703 294 : ABI_MALLOC(ptsymrel,(3,3,msym))
704 294 : ABI_MALLOC(symafm_new,(msym))
705 196 : ABI_MALLOC(symrel_new,(3,3,msym))
706 294 : ABI_MALLOC(tnons_new,(3,msym))
707 98 : call symlatt(bravais,dev_null,msym,nptsym,ptsymrel,rprimd,tolsym)
708 98 : use_inversion=1
709 98 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
710 :
711 : !write(std_out,'(a)')' m_inkpts%inqpt : before symfind '
712 :
713 : call symfind(gprimd,msym,natom,nptsym,1,nsym_new,0,&
714 98 : ptsymrel,spinat,symafm_new,symrel_new,tnons_new,tolsym,typat,use_inversion,xred)
715 :
716 : !write(std_out,'(a)')' m_inkpts%inqpt : after symfind '
717 :
718 : ! Prepare to compute the q-point grid in the ZB or IZB
719 98 : iscf_fake=0 ! Do not need the weights
720 :
721 : ! Compute the maximum number of q points
722 98 : nqpt_max=0
723 98 : ABI_MALLOC(qpts,(3,nqpt_max))
724 98 : ABI_MALLOC(wtq,(nqpt_max))
725 : call getkgrid(chksymbreak,0,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,&
726 : msym,nqpt_max,nqpt_computed,nshiftq,nsym_new,rprimd,&
727 98 : shiftq,symafm_new,symrel_new,vacuum,wtq)
728 :
729 98 : nqpt_max=nqpt_computed
730 98 : ABI_FREE(qpts)
731 98 : ABI_FREE(wtq)
732 :
733 : ! Find the index of the q point within the set of q points that will be generated
734 98 : iqpt=0
735 98 : call intagm(dprarr,intarr,jdtset,marr,1,string(1:lenstr),'iqpt',tread,'INT')
736 98 : tread_q_sum=tread_q_sum+tread
737 98 : if(tread==1)iqpt=intarr(1)
738 :
739 : ! Checks that iqpt is among the computed q points
740 98 : if(iqpt<0)then
741 : write(msg, '(a,i0,3a)' )&
742 0 : 'The input variable iqpt,',iqpt,' is negative, while it should be 0 or positive.',ch10,&
743 0 : 'Action: correct iqpt in the input file.'
744 0 : ABI_ERROR(msg)
745 : end if
746 :
747 98 : if (iqpt > nqpt_computed) then
748 : write(msg, '(a,i0,a,i0,7a)' )&
749 0 : 'The input variable iqpt:',iqpt,' is bigger than the computed number of q-points in the grid which is ',nqpt_max,'.',ch10,&
750 0 : 'The latter has been computed from the input variables qptrlatt, ngqpt, nshiftq,',ch10,&
751 0 : 'shiftq, as well as qptopt, the symmetries of the lattice, and spinat.',ch10,&
752 0 : 'Action: correct iqpt in the input file, or correct the computed q-point grid.'
753 0 : ABI_ERROR(msg)
754 : end if
755 :
756 : ! Compute the q-point grid in the BZ or the IBZ
757 294 : ABI_MALLOC(qpts,(3,nqpt_max))
758 294 : ABI_MALLOC(wtq,(nqpt_max))
759 :
760 : call getkgrid(chksymbreak,iout,iscf_fake,qpts,qptopt,qptrlatt,qptrlen,&
761 : msym,nqpt_max,nqpt_computed,nshiftq,nsym_new,rprimd,&
762 98 : shiftq,symafm_new,symrel_new,vacuum,wtq)
763 :
764 : ! Transfer to qptn, and deallocate
765 98 : qptn(:)=zero
766 98 : if(iqpt/=0)then
767 392 : qptn(:)=qpts(:,iqpt)
768 98 : wtqc = wtq(iqpt)
769 : end if
770 :
771 98 : ABI_FREE(ptsymrel)
772 98 : ABI_FREE(symafm_new)
773 98 : ABI_FREE(symrel_new)
774 98 : ABI_FREE(tnons_new)
775 98 : ABI_FREE(qpts)
776 196 : ABI_FREE(wtq)
777 :
778 : else
779 : write(msg, '(3a,i0,3a)' ) &
780 0 : 'The only values of qptopt allowed are smaller than 4.',ch10,&
781 0 : 'The input value of qptopt is',qptopt,'.',ch10,&
782 0 : 'Action: change qptopt in your input file.'
783 0 : ABI_ERROR(msg)
784 : end if
785 :
786 : ! See issue #31 on gitlab. Not really a good idea.
787 : !if(nqpt==0 .and. tread_q_sum/=0)then
788 : ! write(msg, '(5a)' ) &
789 : ! 'When nqpt is zero, the following input variables cannot be defined :',ch10, &
790 : ! ' iqpt, ngqpt, nshiftq, qptopt, qpt, qptnrm, qptrlatt, shiftq, wtq . ',ch10, &
791 : ! 'Action: change nqpt to 1, or un-define all the variables above.'
792 : ! ABI_ERROR(msg)
793 : !endif
794 :
795 713 : ABI_FREE(intarr)
796 713 : ABI_FREE(dprarr)
797 :
798 713 : end subroutine inqpt
799 : !!***
800 :
801 : end module m_inkpts
802 : !!***
|
