Branch data 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 : : !!***
|
