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 : : !! nsppol = number of spin channels. For nsppol 2, the nesting factor is the average of the two spin channels
147 : : !! so that a non-magnetic system gives the same result with nsppol 1 and nsppol 2.
148 : : !! weight(nband,nkpt,nsppol) = integration weights for each k-point, band and spin
149 : : !! nqpath = number of points requested along the trajectory
150 : : !! qpath_vertices = vertices of the reciprocal space trajectory
151 : : !! base_name = prefix of the output file
152 : : !! gprimd(3,3) dimensional reciprocal lattice vectors
153 : : !! gmet = metric in reciprocal space
154 : : !! prtnest = flags governing the format of the output file
155 : : !! [ndivsm] = number of divisions for the smallest segment of the q-path. Default: 20
156 : : !!
157 : : !! OUTPUT
158 : : !! Write data to file.
159 : : !!
160 : : !! SOURCE
161 : :
162 : 4 : subroutine mknesting(nkpt,kpt,kptrlatt,nband,nsppol,weight,nqpath,&
163 : 4 : qpath_vertices,nqptfull,qptfull,base_name,gprimd,gmet,prtnest,qptrlatt,&
164 : : nsym,symrec,ndivsm) ! optional
165 : :
166 : : !Arguments ------------------------------------
167 : : !scalars
168 : : integer,intent(in) :: nband,nsppol,nkpt,nqpath,prtnest
169 : : integer, intent(in) :: nqptfull
170 : : integer, intent(in), optional :: nsym, ndivsm
171 : : character(len=*),intent(in) :: base_name
172 : : !arrays
173 : : integer,intent(in) :: kptrlatt(3,3)
174 : : integer,intent(in),optional :: symrec(3,3,*)
175 : : real(dp),intent(in) :: gprimd(3,3),kpt(3,nkpt)
176 : : real(dp),intent(in) :: qptfull(3,nqptfull)
177 : : real(dp),intent(in) :: gmet(3,3)
178 : : real(dp),intent(in) :: qpath_vertices(3,nqpath)
179 : : real(dp),intent(in) :: weight(nband,nkpt,nsppol)
180 : : integer,intent(in) :: qptrlatt(3,3)
181 : :
182 : : !Local variables-------------------------------
183 : : !scalars
184 : : integer :: ikpt,jkpt,spin
185 : : integer :: ik1, ik2, ik3, nkptfull
186 : : character(len=500) :: msg
187 : 4 : type(krank_t) :: krank
188 : : !arrays
189 : 4 : integer,allocatable :: tmprank(:),ktable(:)
190 : : character(len=fnlen) :: tmpname
191 : 4 : real(dp),allocatable :: nestfactor(:), nestordered(:), kptfull(:,:), nest_spin(:)
192 : : ! *************************************************************************
193 : :
194 : : if (kptrlatt(1,2) /= 0 .or. kptrlatt(1,3) /= 0 .or. kptrlatt(2,1) /= 0 .or. &
195 [ + - + - : 4 : kptrlatt(2,3) /= 0 .or. kptrlatt(3,1) /= 0 .or. kptrlatt(3,2) /= 0 ) then
+ - + - +
- - + ]
196 : : write (msg,'(4a)')&
197 : 0 : 'kptrlatt should be diagonal in order to calculate the nesting factor,',ch10,&
198 : 0 : 'skipping the nesting factor calculation ',ch10
199 : 0 : ABI_WARNING(msg)
200 : 0 : return
201 : : end if
202 : :
203 [ - + ]: 4 : if (prtnest /= 1 .and. prtnest /= 2) then
204 : 0 : ABI_BUG('prtnest should be 1 or 2')
205 : : end if
206 : :
207 : : !write(msg,'(a,9(i0,1x))')' mknesting : kptrlatt = ',kptrlatt
208 : : !call wrtout(std_out,msg,'COLL')
209 : :
210 : 4 : nkptfull = kptrlatt(1,1)*kptrlatt(2,2)*kptrlatt(3,3)
211 [ - + + - : 972 : ABI_CALLOC(nestordered, (nkptfull))
+ - + - -
+ - + +
+ ]
212 [ - + + - : 12 : ABI_MALLOC(kptfull,(3,nkptfull))
+ - + - -
+ - + ]
213 : :
214 : 28 : ikpt = 0
215 [ + + ]: 28 : do ik3 = 0, kptrlatt(3,3)-1
216 [ + + ]: 172 : do ik2 = 0, kptrlatt(2,2)-1
217 [ + + ]: 1128 : do ik1 = 0, kptrlatt(1,1)-1
218 : 960 : ikpt = ikpt+1
219 [ + + ]: 3984 : kptfull(:,ikpt) = [dble(ik1)/dble(kptrlatt(1,1)), dble(ik2)/dble(kptrlatt(2,2)),dble(ik3)/dble(kptrlatt(3,3))]
220 : : end do
221 : : end do
222 : : end do
223 : :
224 : : ! NOTE: input weights are not normalised, the normalisation factor in introduced in bfactor
225 : : ! new version now puts kptfull in correct order before bfactor, so no need to re-order...
226 [ + + ]: 4 : if (present(symrec)) then
227 [ - + ]: 1 : ABI_CHECK(present(nsym), "error - provide nsym and symrec arguments together")
228 : 1 : call krank%init(nkpt, kpt, nsym=nsym, symrec=symrec)
229 : : else
230 : 3 : call krank%init(nkpt, kpt)
231 : : end if
232 : :
233 : : ! Average over spins.
234 [ - + + - : 8 : ABI_MALLOC(nest_spin, (nkptfull))
- + - + ]
235 [ + + ]: 8 : do spin=1,nsppol
236 : 4 : call bfactor(nkptfull, kptfull, nkptfull, kptfull, krank, nkpt, weight(:,:,spin), nband, nest_spin)
237 [ - + + + ]: 972 : nestordered = nestordered + nest_spin / nsppol
238 : : end do
239 [ - + ]: 4 : ABI_FREE(nest_spin)
240 : :
241 : : !================================================================================================
242 : : !use linear interpolation to plot the bfactor along the given q-path
243 : : ! 1) order the kpoints of the grid putting them in increasing x, then y, then z (FORTRAN convention)
244 : : ! 2) make table from input kpts to ordered kpts
245 : : ! 3) perform interpolation
246 : : !================================================================================================
247 : :
248 : 4 : call outnesting(base_name,gmet,gprimd,kptrlatt,nestordered,nkptfull,nqpath,prtnest,qpath_vertices,ndivsm=ndivsm)
249 [ - + ]: 4 : ABI_FREE(nestordered)
250 : :
251 : : ! Now do the same, but for the nesting factor over the phonon qpoints only
252 [ - + + - : 37 : ABI_CALLOC(nestfactor, (nqptfull))
+ - + - -
+ - + +
+ ]
253 [ - + + - : 8 : ABI_MALLOC(nest_spin, (nqptfull))
- + - + ]
254 [ + + ]: 8 : do spin=1,nsppol
255 : 4 : call bfactor(nkptfull,kptfull,nqptfull,qptfull,krank,nkpt,weight(:,:,spin),nband,nest_spin)
256 [ - + + + ]: 37 : nestfactor = nestfactor + nest_spin / nsppol
257 : : end do
258 [ - + ]: 4 : ABI_FREE(nest_spin)
259 : :
260 : 4 : call krank%free()
261 [ - + ]: 4 : ABI_FREE(kptfull)
262 : :
263 : 4 : call krank%init(nqptfull, qptfull)
264 : :
265 [ - + + - : 12 : ABI_MALLOC(ktable,(nqptfull))
+ - + - -
+ - + ]
266 [ + + ]: 29 : do ikpt=1,nqptfull
267 : 29 : ktable(ikpt) = ikpt
268 : : end do
269 : :
270 [ - + + - : 8 : ABI_MALLOC(tmprank, (nqptfull))
- + - + ]
271 [ + + ]: 29 : do ikpt=1,nqptfull
272 : 29 : tmprank(ikpt) = krank%get_rank(qptfull(:,ikpt))
273 : : end do
274 : 4 : call sort_int(nqptfull, tmprank, ktable)
275 [ - + ]: 4 : ABI_FREE(tmprank)
276 : 4 : call krank%free()
277 : :
278 : : !fill the datagrid for the nesting factor using the Fortran convention and the conventional unit cell
279 : : !NOTE: the Fortran convention is a must if we want to plot the data
280 : : !in the BXSF format, useful for the linear interpolation since we use interpol3d_0d.F90
281 [ - + + - : 8 : ABI_MALLOC(nestordered,(nqptfull))
- + - + ]
282 [ + + ]: 29 : nestordered(:)=zero
283 [ + + ]: 29 : do jkpt=1,nqptfull
284 : 25 : ikpt = ktable(jkpt)
285 : 29 : nestordered(ikpt)=nestfactor(jkpt)
286 : : end do
287 [ - + ]: 4 : ABI_FREE(nestfactor)
288 [ - + ]: 4 : ABI_FREE(ktable)
289 : :
290 [ + - ]: 4 : tmpname = trim(base_name)//"kplusq"
291 : 4 : call outnesting(tmpname,gmet,gprimd,qptrlatt,nestordered,nqptfull,nqpath,prtnest,qpath_vertices,ndivsm=ndivsm)
292 : :
293 [ - + ]: 4 : ABI_FREE(nestordered)
294 : :
295 [ - + - + : 4 : end subroutine mknesting
- + ]
296 : : !!***
297 : :
298 : : !----------------------------------------------------------------------
299 : :
300 : : !!****f* m_nesting/outnesting
301 : : !! NAME
302 : : !! outnesting
303 : : !!
304 : : !! FUNCTION
305 : : !! Write ou the nesting factors calculated in mknesting
306 : : !! Data on file in the X-Y format (prtnest 1) or
307 : : !! in the XCrysden format (XSF) (prtnest 2)
308 : : !!
309 : : !! INPUTS
310 : : !! base_name = prefix of the output file
311 : : !! gmet = metric in reciprocal space
312 : : !! gprimd(3,3) dimensional reciprocal lattice vectors
313 : : !! kptrlatt(3,3) basis vectors for k-grid
314 : : !! nestordered = nesting function on full grid, points ordered in x, then y, then z
315 : : !! nkpt = number of k points
316 : : !! nqpath = number of points requested along the trajectory
317 : : !! prtnest = flags governing the format of the output file
318 : : !! qpath_vertices = vertices of the reciprocal space trajectory
319 : : !! [ndivsm] = number of divisions for the smallest segment of the q-path. Default: 20
320 : : !!
321 : : !! OUTPUT
322 : : !! only write to file
323 : : !!
324 : : !! SOURCE
325 : :
326 : 8 : subroutine outnesting(base_name,gmet,gprimd,kptrlatt,nestordered,nkpt,nqpath,prtnest,qpath_vertices,ndivsm)
327 : :
328 : : !Arguments ------------------------------------
329 : : integer,intent(in) :: nqpath,prtnest,nkpt
330 : : integer,optional,intent(in) :: ndivsm
331 : : character(len=*),intent(in) :: base_name
332 : : !arrays
333 : : integer,intent(in) :: kptrlatt(3,3)
334 : : real(dp),intent(in) :: gprimd(3,3), gmet(3,3)
335 : : real(dp),intent(in) :: qpath_vertices(3,nqpath), nestordered(nkpt)
336 : :
337 : : !Local variables-------------------------------
338 : : !scalars
339 : : integer :: unit_nest,nkx,nky,nkz,indx,ii,ipoint,npt_tot,realrecip,my_ndivsm
340 : : character(len=fnlen) :: fname
341 : : character(len=500) :: msg
342 : : real(dp) :: res(3), kval
343 : : !arrays
344 : 16 : integer :: ndiv(nqpath-1)
345 : 8 : real(dp),allocatable :: finepath(:,:)
346 : : real(dp) :: tmpkpt(3),origin(3),qpt(3)
347 : : ! dummy variables for call to printxsf
348 : : integer :: natom, ntypat, typat(1)
349 : : real(dp) :: xcart (3,1), znucl(1)
350 : : ! *************************************************************************
351 : :
352 : : ! Definition of the q path along which ph linwid will be interpolated
353 [ - + ]: 8 : my_ndivsm = 20; if (present(ndivsm)) my_ndivsm = ndivsm
354 : 8 : call make_path(nqpath,qpath_vertices,gmet,'G',my_ndivsm,ndiv,npt_tot,finepath)
355 : :
356 : 8 : nkx = kptrlatt(1,1)
357 : 8 : nky = kptrlatt(2,2)
358 : 8 : nkz = kptrlatt(3,3)
359 : :
360 [ - + ]: 8 : if (nkpt /= nkx*nky*nkz) then
361 : 0 : write(msg,'(a,9(i0,1x),2x,i0)')' Wrong input value for kptrlatt ',kptrlatt, nkpt
362 : 0 : ABI_BUG(msg)
363 : : end if
364 : :
365 : : ! Open output file and write header
366 [ - + ]: 8 : if (open_file(base_name,msg,newunit=unit_nest,status="unknown",form="formatted",action="write") /= 0) then
367 : 0 : ABI_ERROR(msg)
368 : : end if
369 : :
370 : 8 : write(unit_nest,'(a)')'#'
371 : 8 : write(unit_nest,'(a)')'# ABINIT package : Nesting factor file'
372 : 8 : write(unit_nest,'(a)')'#'
373 : 8 : write(unit_nest,'(a,i10,a)')'# Nesting factor calculated on ',npt_tot,' Q-points'
374 : 8 : write(unit_nest,'(a)')'# Description of the Q-path :'
375 : 8 : write(unit_nest,'(a,i10)')'# Number of line segments = ',nqpath-1
376 : 8 : write(unit_nest,'(a)')'# Vertices of the Q-path and corresponding index = '
377 : 8 : indx=1
378 [ + + ]: 66 : do ii=1,nqpath
379 : 58 : write(unit_nest,'(a,3(E16.6,1x),i8)')'# ',qpath_vertices(:,ii),indx
380 [ + + ]: 66 : if(ii<nqpath) indx=indx+ndiv(ii)
381 : : end do
382 : 8 : write(unit_nest,'(a)')'#'
383 : :
384 : : !Get qpoint along the q-path from finepath and interpolate the nesting factor
385 : 8 : indx=1
386 : :
387 : 8 : write (unit_nest,'(a)')'# index nesting, qfrac_coords'
388 [ + + ]: 1592 : do ipoint=1, npt_tot
389 [ + + ]: 6336 : qpt(:) = finepath(:,ipoint)
390 [ + + ]: 6336 : call wrap2_zero_one(qpt, tmpkpt, res)
391 : 1584 : kval = interpol3d_0d(tmpkpt, nkx, nky, nkz, nestordered)
392 : 1584 : write(unit_nest,'(i5,e16.5,1x,3(es11.4,1x))')indx,kval,tmpkpt
393 : 1592 : indx = indx + 1
394 : : end do
395 : :
396 : 8 : close (unit_nest)
397 [ - + ]: 8 : ABI_FREE(finepath)
398 : :
399 [ - + ]: 8 : if (prtnest==2) then
400 : : ! write also the nesting factor in the XSF format
401 [ # # ]: 0 : fname = trim(base_name) // '_NEST_XSF'
402 : :
403 [ # # ]: 0 : if (open_file(fname,msg,newunit=unit_nest,status="unknown",form="formatted",action="write") /= 0) then
404 : 0 : ABI_ERROR(msg)
405 : : end if
406 : :
407 : 0 : origin(:) = zero
408 : 0 : realrecip = 1 !reciprocal space
409 : 0 : natom = 1
410 : 0 : ntypat = 1
411 : 0 : typat = [1]
412 : 0 : xcart = reshape ([zero, zero, zero], [3, 1])
413 : 0 : znucl = [one]
414 : 0 : call printxsf(nkx,nky,nkz,nestordered,gprimd,origin,natom, ntypat, typat, xcart, znucl, unit_nest,realrecip)
415 : :
416 : 0 : close (unit_nest)
417 : : end if
418 : :
419 : 8 : end subroutine outnesting
420 : : !!***
421 : :
422 : : end module m_nesting
423 : : !!***
|
