Line data Source code
1 : !!****m* ABINIT/m_nesting
2 : !! NAME
3 : !! m_nesting
4 : !!
5 : !! FUNCTION
6 : !! This module provides functions to compute the nesting factor:
7 : !! N(\qq) = \sum_{mn\kk} \delta(\ee_{\kpq m}) \delta(\ee_{\kk n})
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (MG, MJV)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_nesting
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_abicore
28 : use m_krank
29 : use m_sort
30 :
31 : use m_numeric_tools, only : wrap2_zero_one, interpol3d_0d
32 : use m_io_tools, only : open_file
33 : use m_bz_mesh, only : make_path
34 : use m_pptools, only : printxsf
35 :
36 : implicit none
37 :
38 : private
39 :
40 : public :: bfactor
41 : public :: mknesting
42 : public :: outnesting
43 : !!***
44 :
45 : !----------------------------------------------------------------------
46 :
47 : CONTAINS !=============================================================================
48 : !!***
49 :
50 : !!****f* m_nesting/bfactor
51 : !! NAME
52 : !! bfactor
53 : !!
54 : !! FUNCTION
55 : !! Calculate the nesting factor
56 : !!
57 : !! INPUTS
58 : !! nkptfull = number of k-points in full grid
59 : !! kptfull(3,nkptfull) = k-point grid
60 : !! nqpt = number of qpoints
61 : !! qpt(3,nqpt) = q-point grid (must be a subgrid of the k grid),
62 : !! the nesting factor will be calculated for each q point in this array
63 : !! nkpt = eventually reduced number of k-points
64 : !! weight(nband,nkpt) = integration weights for each k-point and band (NOT NORMALISED!!!)
65 : !! nband = number of bands
66 : !!
67 : !! OUTPUT
68 : !! nestfactor(nqpt) = array containing the nesting factor values
69 : !!
70 : !! NOTES
71 : !! Inspired to nmsq_gam_sumfs and mkqptequiv
72 : !! TODO : better use of symmetries to reduce the computational effort
73 : !! Must be called with kpt = full grid! Reduction by symmetry is not possible for q-dependent quantities (or not easy :)
74 : !!
75 : !! SOURCE
76 :
77 23 : subroutine bfactor(nkptfull,kptfull,nqpt,qpt,krank,nkpt,weight,nband,nestfactor)
78 :
79 : !Arguments ------------------------------------
80 : !scalars
81 : integer,intent(in) :: nband,nkptfull,nqpt,nkpt
82 : !arrays
83 : real(dp),intent(in) :: kptfull(3,nkptfull),qpt(3,nqpt),weight(nband,nkpt)
84 : real(dp),intent(out) :: nestfactor(nqpt)
85 : type(krank_t), intent(in) :: krank
86 :
87 : !Local variables-------------------------------
88 : !scalars
89 : integer :: ib1,ib2,ikplusq_irr,ikpt
90 : integer :: irank_kpt,ikpt_irr,iqpt,symrank_kpt
91 : real(dp) :: w1,w2
92 : !character(len=500) :: msg
93 : !arrays
94 : real(dp) :: kptpq(3)
95 : ! *************************************************************************
96 :
97 1056 : nestfactor(:)=zero
98 :
99 1056 : do iqpt=1,nqpt
100 364944 : do ikpt=1,nkptfull
101 363888 : irank_kpt = krank%get_rank(kptfull(:,ikpt))
102 363888 : ikpt_irr = krank%invrank(irank_kpt)
103 :
104 1455552 : kptpq(:) = kptfull(:,ikpt) + qpt(:,iqpt)
105 363888 : symrank_kpt = krank%get_rank(kptpq)
106 :
107 363888 : ikplusq_irr = krank%invrank(symrank_kpt)
108 363888 : if (ikplusq_irr == -1) then
109 0 : ABI_ERROR('It looks like no kpoint equiv to k+q!')
110 : end if
111 :
112 2286441 : do ib1=1,nband
113 1921520 : w1 = weight(ib1, ikpt_irr) ! weight for distance from the Fermi surface
114 1921520 : if (w1 < tol6 ) cycle
115 3912420 : do ib2=1,nband
116 2973720 : w2 = weight(ib2, ikplusq_irr) ! weight for distance from the Fermi surface
117 : if (w1 < tol6 ) cycle
118 4895240 : nestfactor(iqpt) = nestfactor(iqpt) + w1 * w2
119 : end do !ib2
120 : end do !ib1
121 :
122 : end do !ikpt
123 : end do !iqpt
124 :
125 : ! need prefactor of (1/nkptfull) for normalisation of integration
126 1056 : nestfactor(:) = (one/nkptfull) * nestfactor(:)
127 :
128 23 : end subroutine bfactor
129 : !!***
130 :
131 : !----------------------------------------------------------------------
132 :
133 : !!****f* m_nesting/mknesting
134 : !! NAME
135 : !! mknesting
136 : !!
137 : !! FUNCTION
138 : !! Calculate the nesting factor over the dense k-grid, interpolate the values along a given q path
139 : !! and write the data on file in the X-Y format or in the XCrysden format (XSF)
140 : !!
141 : !! INPUTS
142 : !! nkpt = number of k points
143 : !! kpt(3,nkpt) = k points
144 : !! nkx, nky, nkz = number of k-point along each direction
145 : !! nband = number of bands to be considered in the calculation
146 : !! weight(nband,nkpt) = integration weights for each k-point and band
147 : !! nqpath = number of points requested along the trajectory
148 : !! qpath_vertices = vertices of the reciprocal space trajectory
149 : !! base_name = prefix of the output file
150 : !! gprimd(3,3) dimensional reciprocal lattice vectors
151 : !! gmet = metric in reciprocal space
152 : !! prtnest = flags governing the format of the output file
153 : !!
154 : !! OUTPUT
155 : !! Write data to file.
156 : !!
157 : !! SOURCE
158 :
159 4 : subroutine mknesting(nkpt,kpt,kptrlatt,nband,weight,nqpath,&
160 4 : qpath_vertices,nqptfull,qptfull,base_name,gprimd,gmet,prtnest,qptrlatt,&
161 : nsym,symrec) ! optional
162 :
163 : !Arguments ------------------------------------
164 : !scalars
165 : integer,intent(in) :: nband,nkpt,nqpath,prtnest
166 : integer, intent(in) :: nqptfull
167 : integer, intent(in), optional :: nsym
168 : character(len=*),intent(in) :: base_name
169 : !arrays
170 : integer,intent(in) :: kptrlatt(3,3)
171 : integer,intent(in),optional :: symrec(3,3,*)
172 : real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt)
173 : real(dp),intent(in) :: qptfull(3,nqptfull)
174 : real(dp),intent(in) :: gmet(3,3)
175 : real(dp),intent(in) :: qpath_vertices(3,nqpath)
176 : real(dp),intent(in) :: weight(nband,nkpt)
177 : integer,intent(in) :: qptrlatt(3,3)
178 :
179 : !Local variables-------------------------------
180 : !scalars
181 : integer :: ikpt,jkpt
182 : integer :: ik1, ik2, ik3, nkptfull
183 : character(len=500) :: msg
184 4 : type(krank_t) :: krank
185 : !arrays
186 4 : integer,allocatable :: tmprank(:),ktable(:)
187 : character(len=fnlen) :: tmpname
188 4 : real(dp),allocatable :: nestfactor(:), nestordered(:), kptfull(:,:)
189 : ! *************************************************************************
190 :
191 : if (kptrlatt(1,2) /= 0 .or. kptrlatt(1,3) /= 0 .or. kptrlatt(2,1) /= 0 .or. &
192 4 : kptrlatt(2,3) /= 0 .or. kptrlatt(3,1) /= 0 .or. kptrlatt(3,2) /= 0 ) then
193 : write (msg,'(4a)')&
194 0 : 'kptrlatt should be diagonal in order to calculate the nesting factor,',ch10,&
195 0 : 'skipping the nesting factor calculation ',ch10
196 0 : ABI_WARNING(msg)
197 0 : return
198 : end if
199 :
200 4 : if (prtnest /= 1 .and. prtnest /= 2) then
201 0 : ABI_BUG('prtnest should be 1 or 2')
202 : end if
203 :
204 : !write(msg,'(a,9(i0,1x))')' mknesting : kptrlatt = ',kptrlatt
205 : !call wrtout(std_out,msg,'COLL')
206 :
207 4 : nkptfull = kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3)
208 972 : ABI_CALLOC(nestordered, (nkptfull))
209 12 : ABI_MALLOC(kptfull,(3,nkptfull))
210 :
211 28 : ikpt = 0
212 28 : do ik3 = 0, kptrlatt(3,3)-1
213 172 : do ik2 = 0, kptrlatt(2,2)-1
214 1128 : do ik1 = 0, kptrlatt(1,1)-1
215 960 : ikpt = ikpt+1
216 3984 : kptfull(:,ikpt) = [dble(ik1)/dble(kptrlatt(1,1)), dble(ik2)/dble(kptrlatt(2,2)),dble(ik3)/dble(kptrlatt(3,3))]
217 : end do
218 : end do
219 : end do
220 :
221 : ! NOTE: input weights are not normalised, the normalisation factor in introduced in bfactor
222 : ! new version now puts kptfull in correct order before bfactor, so no need to re-order...
223 4 : if (present(symrec)) then
224 1 : ABI_CHECK(present(nsym), "error - provide nsym and symrec arguments together")
225 1 : call krank%init(nkpt, kpt, nsym=nsym, symrec=symrec)
226 : else
227 3 : call krank%init(nkpt, kpt)
228 : end if
229 :
230 4 : call bfactor(nkptfull, kptfull, nkptfull, kptfull, krank, nkpt, weight, nband, nestordered)
231 :
232 : !================================================================================================
233 : !use linear interpolation to plot the bfactor along the given q-path
234 : ! 1) order the kpoints of the grid putting them in increasing x, then y, then z (FORTRAN convention)
235 : ! 2) make table from input kpts to ordered kpts
236 : ! 3) perform interpolation
237 : !================================================================================================
238 :
239 4 : call outnesting(base_name,gmet,gprimd,kptrlatt,nestordered,nkptfull,nqpath,prtnest,qpath_vertices)
240 4 : ABI_FREE(nestordered)
241 :
242 : ! Now do the same, but for the nesting factor over the phonon qpoints only
243 12 : ABI_MALLOC(nestfactor, (nqptfull))
244 4 : call bfactor(nkptfull,kptfull,nqptfull,qptfull,krank,nkpt,weight,nband,nestfactor)
245 :
246 4 : call krank%free()
247 4 : ABI_FREE(kptfull)
248 :
249 4 : call krank%init(nqptfull, qptfull)
250 :
251 12 : ABI_MALLOC(ktable,(nqptfull))
252 29 : do ikpt=1,nqptfull
253 29 : ktable(ikpt) = ikpt
254 : end do
255 :
256 8 : ABI_MALLOC(tmprank, (nqptfull))
257 29 : do ikpt=1,nqptfull
258 29 : tmprank(ikpt) = krank%get_rank(qptfull(:,ikpt))
259 : end do
260 4 : call sort_int(nqptfull, tmprank, ktable)
261 4 : ABI_FREE(tmprank)
262 4 : call krank%free()
263 :
264 : !fill the datagrid for the nesting factor using the Fortran convention and the conventional unit cell
265 : !NOTE: the Fortran convention is a must if we want to plot the data
266 : !in the BXSF format, useful for the linear interpolation since we use interpol3d_0d.F90
267 8 : ABI_MALLOC(nestordered,(nqptfull))
268 29 : nestordered(:)=zero
269 29 : do jkpt=1,nqptfull
270 25 : ikpt = ktable(jkpt)
271 29 : nestordered(ikpt)=nestfactor(jkpt)
272 : end do
273 4 : ABI_FREE(nestfactor)
274 4 : ABI_FREE(ktable)
275 :
276 4 : tmpname = trim(base_name)//"kplusq"
277 4 : call outnesting(tmpname,gmet,gprimd,qptrlatt,nestordered,nqptfull,nqpath,prtnest,qpath_vertices)
278 :
279 4 : ABI_FREE(nestordered)
280 :
281 4 : end subroutine mknesting
282 : !!***
283 :
284 : !----------------------------------------------------------------------
285 :
286 : !!****f* m_nesting/outnesting
287 : !! NAME
288 : !! outnesting
289 : !!
290 : !! FUNCTION
291 : !! Write ou the nesting factors calculated in mknesting
292 : !! Data on file in the X-Y format (prtnest 1) or
293 : !! in the XCrysden format (XSF) (prtnest 2)
294 : !!
295 : !! INPUTS
296 : !! base_name = prefix of the output file
297 : !! gmet = metric in reciprocal space
298 : !! gprimd(3,3) dimensional reciprocal lattice vectors
299 : !! kptrlatt(3,3) basis vectors for k-grid
300 : !! nestordered = nesting function on full grid, points ordered in x, then y, then z
301 : !! nkpt = number of k points
302 : !! nqpath = number of points requested along the trajectory
303 : !! prtnest = flags governing the format of the output file
304 : !! qpath_vertices = vertices of the reciprocal space trajectory
305 : !!
306 : !! OUTPUT
307 : !! only write to file
308 : !!
309 : !! SOURCE
310 :
311 8 : subroutine outnesting(base_name,gmet,gprimd,kptrlatt,nestordered,nkpt,nqpath,prtnest,qpath_vertices)
312 :
313 : !Arguments ------------------------------------
314 : integer,intent(in) :: nqpath,prtnest,nkpt
315 : character(len=*),intent(in) :: base_name
316 : !arrays
317 : integer,intent(in) :: kptrlatt(3,3)
318 : real(dp),intent(in) :: gprimd(3,3), gmet(3,3)
319 : real(dp),intent(in) :: qpath_vertices(3,nqpath), nestordered(nkpt)
320 :
321 : !Local variables-------------------------------
322 : !scalars
323 : integer :: unit_nest,nkx,nky,nkz,indx,ii,ipoint,npt_tot,realrecip
324 : character(len=fnlen) :: fname
325 : character(len=500) :: msg
326 : real(dp) :: res(3), kval
327 : !arrays
328 16 : integer :: ndiv(nqpath-1)
329 8 : real(dp),allocatable :: finepath(:,:)
330 : real(dp) :: tmpkpt(3),origin(3),qpt(3)
331 : ! dummy variables for call to printxsf
332 : integer :: natom, ntypat, typat(1)
333 : real(dp) :: xcart (3,1), znucl(1)
334 : ! *************************************************************************
335 :
336 : ! Definition of the q path along which ph linwid will be interpolated
337 8 : call make_path(nqpath,qpath_vertices,gmet,'G',20,ndiv,npt_tot,finepath)
338 :
339 8 : nkx = kptrlatt(1,1)
340 8 : nky = kptrlatt(2,2)
341 8 : nkz = kptrlatt(3,3)
342 :
343 8 : if (nkpt /= nkx*nky*nkz) then
344 0 : write(msg,'(a,9(i0,1x),2x,i0)')' Wrong input value for kptrlatt ',kptrlatt, nkpt
345 0 : ABI_BUG(msg)
346 : end if
347 :
348 : ! Open output file and write header
349 8 : if (open_file(base_name,msg,newunit=unit_nest,status="unknown",form="formatted",action="write") /= 0) then
350 0 : ABI_ERROR(msg)
351 : end if
352 :
353 8 : write(unit_nest,'(a)')'#'
354 8 : write(unit_nest,'(a)')'# ABINIT package : Nesting factor file'
355 8 : write(unit_nest,'(a)')'#'
356 8 : write(unit_nest,'(a,i10,a)')'# Nesting factor calculated on ',npt_tot,' Q-points'
357 8 : write(unit_nest,'(a)')'# Description of the Q-path :'
358 8 : write(unit_nest,'(a,i10)')'# Number of line segments = ',nqpath-1
359 8 : write(unit_nest,'(a)')'# Vertices of the Q-path and corresponding index = '
360 8 : indx=1
361 66 : do ii=1,nqpath
362 58 : write(unit_nest,'(a,3(E16.6,1x),i8)')'# ',qpath_vertices(:,ii),indx
363 66 : if(ii<nqpath) indx=indx+ndiv(ii)
364 : end do
365 8 : write(unit_nest,'(a)')'#'
366 :
367 : !Get qpoint along the q-path from finepath and interpolate the nesting factor
368 8 : indx=1
369 :
370 8 : write (unit_nest,'(a)')'# index nesting, qfrac_coords'
371 1592 : do ipoint=1, npt_tot
372 6336 : qpt(:) = finepath(:,ipoint)
373 6336 : call wrap2_zero_one(qpt, tmpkpt, res)
374 1584 : kval = interpol3d_0d(tmpkpt, nkx, nky, nkz, nestordered)
375 1584 : write(unit_nest,'(i5,e16.5,1x,3(es11.4,1x))')indx,kval,tmpkpt
376 1592 : indx = indx + 1
377 : end do
378 :
379 8 : close (unit_nest)
380 8 : ABI_FREE(finepath)
381 :
382 8 : if (prtnest==2) then
383 : ! write also the nesting factor in the XSF format
384 0 : fname = trim(base_name) // '_NEST_XSF'
385 :
386 0 : if (open_file(fname,msg,newunit=unit_nest,status="unknown",form="formatted",action="write") /= 0) then
387 0 : ABI_ERROR(msg)
388 : end if
389 :
390 0 : origin(:) = zero
391 0 : realrecip = 1 !reciprocal space
392 0 : natom = 1
393 0 : ntypat = 1
394 0 : typat = [1]
395 0 : xcart = reshape ([zero, zero, zero], [3, 1])
396 0 : znucl = [one]
397 0 : call printxsf(nkx,nky,nkz,nestordered,gprimd,origin,natom, ntypat, typat, xcart, znucl, unit_nest,realrecip)
398 :
399 0 : close (unit_nest)
400 : end if
401 :
402 8 : end subroutine outnesting
403 : !!***
404 :
405 : end module m_nesting
406 : !!***
|
