Branch data Line data Source code
1 : : !!****m*ABINIT/m_wannier_builder
2 : : !! NAME
3 : : !! m_wannier_builder
4 : : !!
5 : : !! FUNCTION
6 : : !! Algorithms for building Wannier functions
7 : : !! Methods:
8 : : !! SCDM (select columns of density matrix method) and
9 : : !! projected wannier function (PWF)
10 : : !! COPYRIGHT
11 : : !! Copyright (C) 2005-2026 ABINIT group (hexu)
12 : : !! This file is distributed under the terms of the
13 : : !! GNU General Public License, see ~abinit/COPYING
14 : : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : : !!
16 : : !! SOURCE
17 : : !!
18 : : !! Todolist:
19 : : !! - output of Wannier90 format Amn, Mmn, eigenvalue files
20 : : !! - output header in netcdf file
21 : : !! - allow units in Hamiltonian
22 : : !! - projected wannier functions
23 : :
24 : : #if defined HAVE_CONFIG_H
25 : : #include "config.h"
26 : : #endif
27 : :
28 : :
29 : : #include "abi_common.h"
30 : :
31 : :
32 : : !===============================================================
33 : : ! SCDM-k
34 : : !> @description: Select Column Density matrix method for
35 : : ! generating Wannier functions.
36 : : !===============================================================
37 : : module m_wannier_builder
38 : : use defs_basis
39 : : use m_abicore
40 : : use m_errors
41 : : use m_io_tools, only : open_file
42 : : use m_fstrings, only : ltoa
43 : : use m_scdm_math, only: complex_QRCP_piv_only, complex_svd, tpi_im, &
44 : : & gaussian, fermi, insertion_sort_double, eigensolver
45 : : use m_wann_netcdf, only: IOWannNC
46 : : implicit none
47 : : public:: WannierBuilder_t
48 : : public :: WannierBuilder_witheigen_t
49 : : public:: Amn_to_H
50 : : private
51 : :
52 : :
53 : : !===============================================================
54 : : ! WannierBuilder_t type:
55 : : !> @ description: the class for scdmk method.
56 : : !===============================================================
57 : : type:: WannierBuilder_t
58 : :
59 : : ! Inputs
60 : : real(dp), allocatable:: kpts(:, :) !(idim, ikpt)
61 : : real(dp), allocatable:: kweights(:) !(ikpt)
62 : : !real(dp), allocatable:: weight(:, :) !(iband, ikpt)
63 : : integer, allocatable:: exclude_bands(:)
64 : :
65 : : integer:: method = 0 ! 1: SCDM-k 2: projected wf
66 : : ! Disentanglement
67 : : integer:: disentangle_func_type ! 1:unity 2: erfc function 3: gauss function
68 : : real(dp):: mu, sigma
69 : :
70 : : ! Projected Wannier function related
71 : : complex(dp), allocatable:: projectors(:, :)
72 : :
73 : : ! SCDM related
74 : : integer, allocatable:: cols(:)
75 : : real(dp), allocatable:: anchor_kpt(:)
76 : : integer:: anchor_ikpt
77 : : integer, allocatable:: anchor_ibands(:)
78 : : logical:: project_to_anchor = .False.
79 : :
80 : : ! Output Wannier function and Hamiltonian
81 : : integer:: nwann, nkpt, nband, nbasis, nkdim
82 : : integer:: dim ! dimension of position
83 : : integer:: nR
84 : : integer, allocatable:: Rlist(:,:) !(idim, iRpt)
85 : : complex(dp), allocatable:: Amnk(:, :, :) !(nband, nwann, nkpt)
86 : : complex(dp), allocatable:: Hwannk(:, :, :) !(nwann, nwann, nkpt)
87 : :
88 : : complex(dp), allocatable:: psi_wann_k(:, :, :) !(nbasis, nwann, nkpt)
89 : :
90 : : complex(dp), allocatable:: HwannR(:, :, :) !(nwann, nwann, nR)
91 : : complex(dp), allocatable:: wannR(:,:,:) !(nbasis, nwann, nR)
92 : :
93 : : contains
94 : : procedure:: initialize
95 : : procedure:: finalize
96 : : procedure:: get_psi_k
97 : : procedure:: get_evals_k
98 : : procedure:: run_all
99 : : procedure:: find_kpoint
100 : : !procedure:: remove_phase
101 : : procedure:: auto_find_anchors
102 : : procedure:: get_columns
103 : : procedure:: get_scdm_Amnk ! Amnk for all kpoints
104 : : procedure:: get_projected_Amnk ! Amnk for all kpoints
105 : : procedure:: get_Amnk
106 : : procedure:: get_weight
107 : : procedure:: set_anchor
108 : : procedure:: set_disp_projector
109 : : procedure:: set_mode_projector
110 : : procedure:: get_wannR_and_HwannR
111 : : procedure:: select_columns
112 : : procedure:: construct_wannier
113 : : procedure:: write_Amnk_w90
114 : : procedure:: write_Hwann_w90
115 : : procedure:: create_ncfile
116 : : procedure:: close_ncfile
117 : : procedure:: write_wann_netcdf
118 : : procedure:: get_wannier_eigen
119 : : procedure:: get_wannier_eigen_klist
120 : : end type WannierBuilder_t
121 : :
122 : :
123 : : ! Extends WannierBuilder_t with pointer to eigens.
124 : : type, extends(WannierBuilder_t):: WannierBuilder_witheigen_t
125 : : real(dp), pointer:: evals(:, :) => null() !(iband, ikpt)
126 : : complex(dp), pointer:: psi(:, :, :) => null() ! (ibasis, iband, ikpt)
127 : : contains
128 : : procedure :: set_eigen => set_eigen
129 : : procedure :: get_psi_k => get_psi_k_from_eigen
130 : : procedure :: get_evals_k => get_evals_k_from_eigen
131 : : end type WannierBuilder_witheigen_t
132 : :
133 : : contains
134 : :
135 : : !===============================================================
136 : : !
137 : : !> @
138 : : !> kpts: kpoints. indices: (idim, ikpt)
139 : : !> kweights: kweights. indices: (ikpt)
140 : : !> nwann: number of Wannier functions to be calcualted.
141 : : !> nbasis: number of basis in the original wavefunction.
142 : : !> disentangle_func_type: the type of the disentanglement function. 1: unity function. 2. Fermi. 3. Gauss
143 : : !> project_to_anchor: whether to multiply the weight function by the projection to the anchor states.
144 : : !===============================================================
145 [ - + - + : 3 : subroutine initialize(self,kpts, kweights, Rlist, nwann, nbasis, nband, &
- + ]
146 [ + - + - ]: 3 : & disentangle_func_type, mu, sigma, exclude_bands, project_to_anchor, method)
147 : : class(WannierBuilder_t), intent(inout):: self
148 : : integer, intent(in):: nwann, nbasis, nband
149 : : real(dp), intent(in):: kpts(:, :) !(idim, ikpt)
150 : : real(dp), intent(in):: kweights(:) !(ikpt)
151 : : integer, intent(in):: Rlist(:, :) !(idim, iRpt)
152 : : integer, optional, intent(in):: disentangle_func_type
153 : : real(dp), optional, intent(in):: mu, sigma
154 : : integer, optional, intent(in):: exclude_bands(:)
155 : : logical, optional, intent(in):: project_to_anchor
156 : : integer, intent(in):: method
157 : :
158 : 3 : self%nkdim = size(kpts, 1)
159 : 3 : self%nkpt = size(kpts, 2)
160 : 3 : self%nwann = nwann
161 : 3 : self%nbasis = nbasis
162 : 3 : self%nband = nband
163 : 3 : self%method = method
164 : : !if (present(psi_phase)) then
165 : : ! if (psi_phase) then
166 : : ! ABI_MALLOC(self%psi, (self%nbasis, self%nband, self%nkpt))
167 : : ! call self%remove_phase(psi)
168 : : ! else
169 : : ! self%psi => psi
170 : : ! end if
171 : : !else
172 : : ! self%psi => psi
173 : : !end if
174 : :
175 [ - + + - : 9 : ABI_MALLOC(self%cols, (self%nwann))
+ - + - -
+ - + -
+ ]
176 [ + + ]: 12 : self%cols(:) = 0
177 [ - + + - : 12 : ABI_MALLOC(self%kpts, (self%nkdim, self%nkpt))
- + + - +
- + - - +
- + - + ]
178 [ + - + - : 774 : self%kpts = kpts
- + - - +
+ + + ]
179 [ - + + - : 9 : ABI_MALLOC(self%kweights, (self%nkpt))
+ - + - -
+ - + -
+ ]
180 [ + - - + : 198 : self%kweights = kweights
- - + + ]
181 : :
182 : 3 : self%nR = size(Rlist, 2)
183 [ - + + - : 12 : ABI_MALLOC(self%Rlist, (self%nkdim, self%nR))
- + + - +
- + - - +
- + - + ]
184 [ + - + - : 774 : self%Rlist = Rlist
- + - - +
+ + + ]
185 : :
186 : : !ABI_MALLOC(self%weight, (self%nband, self%nkpt))
187 : : !self%weight(:, :) = 0.0_dp
188 : :
189 [ + - ]: 3 : if (present(disentangle_func_type)) then
190 : 3 : self%disentangle_func_type = disentangle_func_type
191 : : else
192 : 0 : self%disentangle_func_type = 0
193 : : end if
194 : :
195 [ + - ]: 3 : if (present(mu)) then
196 : 3 : self%mu = mu
197 : : else
198 : 0 : self%mu = 0
199 : : end if
200 : :
201 [ + - ]: 3 : if (present(sigma)) then
202 : 3 : self%sigma = sigma
203 : : else
204 : 0 : self%sigma = sigma
205 : : end if
206 : :
207 [ - + + - : 15 : ABI_MALLOC(self%Amnk, (self%nband, self%nwann, self%nkpt))
- + + - -
+ + - + -
+ - - + -
+ - + ]
208 [ - + + - : 15 : ABI_MALLOC(self%psi_wann_k, (self%nbasis, self%nwann, self%nkpt))
- + + - -
+ + - + -
+ - - + -
+ - + ]
209 [ - + + - : 15 : ABI_MALLOC(self%Hwannk, (self%nwann, self%nwann, self%nkpt))
- + + - -
+ + - + -
+ - - + -
+ - + ]
210 : :
211 [ + - ]: 3 : if (present(exclude_bands)) then
212 [ + - - - : 9 : ABI_MALLOC(self%exclude_bands, (size(exclude_bands, 1)))
+ - - + -
+ - + -
+ ]
213 : : end if
214 : :
215 [ + - ]: 3 : if (present(project_to_anchor)) self%project_to_anchor = project_to_anchor
216 : 3 : end subroutine initialize
217 : :
218 : :
219 : :
220 : 3 : subroutine finalize(self)
221 : : class(WannierBuilder_t), intent(inout):: self
222 [ + - ]: 3 : ABI_SFREE(self%cols)
223 [ + - ]: 3 : ABI_SFREE(self%kpts)
224 [ + - ]: 3 : ABI_SFREE(self%kweights)
225 [ + - ]: 3 : ABI_SFREE(self%Amnk)
226 [ + - ]: 3 : ABI_SFREE(self%psi_wann_k)
227 [ + - ]: 3 : ABI_SFREE(self%Hwannk)
228 [ + - ]: 3 : ABI_SFREE(self%Rlist)
229 [ + - ]: 3 : ABI_SFREE(self%wannR)
230 [ + - ]: 3 : ABI_SFREE(self%HwannR)
231 [ + - ]: 3 : ABI_SFREE(self%exclude_bands)
232 : 5 : select case (self%method)
233 : : case(1)
234 [ + - ]: 2 : ABI_SFREE(self%anchor_kpt)
235 [ + - ]: 2 : ABI_SFREE(self%anchor_ibands)
236 [ + - ]: 2 : ABI_SFREE(self%projectors)
237 : : case(2)
238 [ + + - + : 3 : ABI_SFREE(self%projectors)
- ]
239 : : end select
240 : 3 : end subroutine finalize
241 : :
242 : 0 : function get_psi_k(self, ikpt) result(psik)
243 : : class(WannierBuilder_t), intent(inout):: self
244 : : integer, intent(in):: ikpt
245 : : complex(dp), pointer:: psik(:, :)
246 : 0 : ABI_UNUSED_A(self)
247 : 0 : ABI_UNUSED_A(ikpt)
248 : 0 : ABI_UNUSED_A(psik)
249 : 0 : ABI_ERROR("WannierBuilder_t%get_psi_k should be overrided!")
250 : 0 : end function get_psi_k
251 : :
252 : :
253 : 0 : function get_evals_k(self, ikpt) result(ek)
254 : : class(WannierBuilder_t), intent(inout):: self
255 : : integer, intent(in):: ikpt
256 : : real(dp), pointer:: ek(:)
257 : 0 : ABI_UNUSED_A(self)
258 : 0 : ABI_UNUSED_A(ikpt)
259 : 0 : ABI_UNUSED_A(ek)
260 : 0 : ABI_ERROR("WannierBuilder_t%get_evals_k should be overrided")
261 : 0 : end function get_evals_k
262 : :
263 : :
264 : : ! automatically set the anchor points using the weight functions.
265 : : ! The bands with the largest weights are selected as the anchor points.
266 : 0 : subroutine auto_find_anchors(self, ianchors)
267 : : class(WannierBuilder_t), intent(inout):: self
268 : : integer, intent(out):: ianchors(self%nwann)
269 : : integer:: i
270 : 0 : real(dp):: weights(self%nband)
271 : 0 : integer:: order(self%nband)
272 : :
273 : : call self%get_weight(self%anchor_ikpt, self%disentangle_func_type, &
274 : 0 : & self%mu, self%sigma, weights, project_to_anchor=.False.)
275 : 0 : call insertion_sort_double(weights, order)
276 : :
277 [ # # ]: 0 : do i = 1, self%nwann
278 : 0 : ianchors(i)= order(self%nband-i+1)
279 : : end do
280 : 0 : end subroutine auto_find_anchors
281 : :
282 [ - + + - : 2 : subroutine set_anchor(self, anchor_kpt, anchor_ibands)
+ - - + ]
283 : : !> anchor_kpt: anchor kpoint (optional).
284 : : !> anchor_ibands: the indices of mode used as anchor points (optional).
285 : : class(WannierBuilder_t), intent(inout):: self
286 : : real(dp), intent(in) :: anchor_kpt(:)
287 : : integer, optional, intent(in):: anchor_ibands(:)
288 : 2 : complex(dp), pointer :: psik(:,:)
289 : : character(len = 500):: msg
290 : : integer :: i
291 [ + - + - : 6 : ABI_MALLOC(self%anchor_ibands, (self%nwann))
+ - - + -
+ - + ]
292 [ - + + - : 6 : ABI_MALLOC(self%anchor_kpt, (size(anchor_kpt)))
+ - + - -
+ - + -
+ ]
293 [ - + + - : 8 : ABI_MALLOC(self%projectors, (self%nbasis, self%nwann))
- + + - +
- + - - +
- + - + ]
294 [ + - - + : 10 : self%anchor_kpt = anchor_kpt
- - + + ]
295 : 2 : self%anchor_ikpt = self%find_kpoint(anchor_kpt)
296 : :
297 [ - + ]: 2 : if (.not. present(anchor_ibands)) then
298 : 0 : call wrtout( std_out, "Anchor points not specified, finding atomatically")
299 : 0 : call self%auto_find_anchors( self%anchor_ibands)
300 [ - + ]: 2 : else if (.not. size(anchor_ibands) == self%nwann) then
301 : 0 : ABI_ERROR("The number of anchor points should be equal to the number of Wannier functions.")
302 : : else
303 [ + - - + : 10 : self%anchor_ibands = anchor_ibands
- - + + ]
304 : : end if
305 [ + - ]: 2 : write(msg, "(2a)") "Anchor point band indices set to ", trim(ltoa(self%anchor_ibands))
306 [ + + ]: 6 : call wrtout([ab_out, std_out], msg )
307 : 2 : psik=> self%get_psi_k(self%anchor_ikpt)
308 [ + + ]: 8 : do i = 1, self%nwann
309 [ + + ]: 98 : self%projectors(:, i)=psik(:, self%anchor_ibands(i))
310 : : end do
311 [ - + ]: 2 : end subroutine set_anchor
312 : :
313 : :
314 [ - + ]: 1 : subroutine set_disp_projector(self, id_projectors)
315 : : class(WannierBuilder_t), intent(inout):: self
316 : : integer, intent(in):: id_projectors(:)
317 : : integer:: i
318 [ - + + - : 4 : ABI_MALLOC(self%projectors, (self%nbasis, self%nwann))
- + + - +
- + - - +
- + - + ]
319 [ + + + + ]: 49 : self%projectors(:,:)=zero
320 [ + + ]: 4 : do i = 1, self%nwann
321 : 4 : self%projectors(id_projectors(i), i)= cmplx(1.0, 0.0, dp)
322 : : end do
323 : 1 : end subroutine set_disp_projector
324 : :
325 : 0 : subroutine set_mode_projector(self, kpoint, iband)
326 : : class(WannierBuilder_t), intent(inout):: self
327 : : real(dp), intent(in):: kpoint(:)
328 : : integer:: iband(:)
329 : 0 : ABI_UNUSED_A(self)
330 : : ABI_UNUSED(kpoint)
331 : : ABI_UNUSED(iband)
332 : 0 : end subroutine set_mode_projector
333 : :
334 : :
335 : 2 : subroutine select_columns(self)
336 : : class(WannierBuilder_t), intent(inout):: self
337 : : integer:: iband
338 : 4 : complex(dp):: psi_dagger(self%nband, self%nbasis)
339 : 4 : complex(dp):: psi_dagger_copy(self%nband, self%nbasis)
340 : 4 : real(dp):: weight(self%nband)
341 : : !real(dp), pointer:: evals_anchor(:)
342 : : !type(eigensolver):: esolver
343 : : character(len = 500):: msg
344 : : ! find anchor points, by default gamma
345 : 2 : self%anchor_ikpt = self%find_kpoint(self%anchor_kpt)
346 : : ! TODO: add qpt if anchor_ikpt is not found
347 : : !
348 : : !if (size(self%anchor_ibands) /= 0) then
349 : : !
350 : : !end if
351 : : ! calculate weight matrix for each kpoint
352 : : call self%get_weight(self%anchor_ikpt, self%disentangle_func_type, self%mu, self%sigma, weight, &
353 : 2 : &project_to_anchor = .True.)
354 : : ! &project_to_anchor = self%project_to_anchor)
355 : :
356 : : ! at anchor-kpoint, find cols
357 : : ! psi: (ibasis, iband)
358 [ + + + + ]: 482 : psi_dagger = transpose(conjg(self%get_psi_k(self%anchor_ikpt)))
359 [ + + ]: 32 : do iband = 1, self%nband
360 [ + + ]: 482 : psi_dagger(iband, :) = psi_dagger(iband, :)*weight(iband)
361 : : end do
362 [ + + + + ]: 482 : psi_dagger_copy(:,:) = psi_dagger(:,:)
363 : 2 : call self%get_columns(psi_dagger_copy, self%cols)
364 [ + - ]: 2 : write(msg, '(2a) ') 'Columns selected: ', trim(ltoa(self%cols))
365 [ + + ]: 6 : call wrtout([ab_out, std_out], msg )
366 : :
367 : : !psi_dagger_copy(:,:) = psi_dagger(:,:)
368 : : !! check the eigen values:
369 : : !call self%get_Amnk(self%anchor_ikpt, Amn)
370 : : !! print the anchor point eigen values:
371 : : !evals_anchor => self%get_evals_k(self%anchor_ikpt)
372 : : !call Amn_to_H_from_evals(Amn, evals_anchor, &
373 : : ! & self%nwann, self%nband, Hwann)
374 : : !evals = evals_anchor(self%anchor_ibands)
375 : : !write(msg, '(2a)') "The eigen values of the anchor points: ", &
376 : : ! & trim(ltoa(evals))
377 : : !call wrtout([ab_out, std_out], msg )
378 : : !! calculate the eigen values of the Hwannk at anchor point
379 : : !call esolver%run(evals, Hwann)
380 : : !write(msg, '(2a)') "The eigen values of Hwann(k_anchor): ", &
381 : : ! & trim(ltoa(evals))
382 : : !call wrtout([ab_out, std_out], msg )
383 : : !call esolver%finalize()
384 : 2 : end subroutine select_columns
385 : :
386 : 3 : subroutine construct_wannier(self)
387 : : class(WannierBuilder_t), intent(inout):: self
388 : : integer:: ikpt
389 : 3 : complex(dp), pointer:: p(:, :)
390 : : !real(dp):: weight(self%nband)
391 : :
392 : : !complex(dp):: tmp(self%nwann, self%nwann)
393 : : !real(dp):: evals(self%nwann)
394 : : !type(eigensolver):: esolver
395 : :
396 : 3 : if(self%method == 1)then
397 : 2 : call self%select_columns()
398 : : end if
399 : :
400 : : ! For each kpoint, calculate Amn matrix, wannier function, and Hk at k
401 [ + + ]: 195 : do ikpt = 1, self%nkpt
402 : : !Amnk (nband, nwann, nkpt)
403 : 192 : p => self%get_psi_k(ikpt)
404 : 192 : call self%get_Amnk(ikpt, self%Amnk(:, :, ikpt))
405 : : ! psik*Amnk
406 [ + - + + : 796416 : self%psi_wann_k(:, :, ikpt) = matmul(p, self%Amnk(:, :, ikpt))
+ + + + +
+ + + + -
+ - + - +
- + - + -
- - - - ]
407 : : call Amn_to_H_from_evals(self%Amnk(:, :, ikpt), self%get_evals_k(ikpt), &
408 [ - + ]: 195 : self%nwann, self%nband, self%Hwannk(:, :, ikpt))
409 : : !tmp = self%Hwannk(:,:,ikpt)
410 : : !call esolver%run(evals, tmp)
411 : : !call esolver%finalize()
412 : : !print *, "ev1:", self%get_evals_k(ikpt)
413 : : ! print *, "ev2:",evals
414 : :
415 : : end do
416 : : ! Fourier transform of wannier function to real space
417 : 3 : call self%get_wannR_and_HwannR(self%Rlist)
418 [ + + ]: 3 : end subroutine construct_wannier
419 : :
420 : :
421 : :
422 : 0 : subroutine run_all(self, ncfilename, Amnkfilename)
423 : : class(WannierBuilder_t), intent(inout):: self
424 : : character(*), intent(in):: ncfilename
425 : : character(*), intent(in):: Amnkfilename
426 : : type(IOWannNC):: ncfile
427 : 0 : call self%construct_wannier()
428 : 0 : call self%create_ncfile(ncfilename, ncfile)
429 : : call self%write_wann_netcdf( ncfile, &
430 : 0 : &wannR_unit='dimensionless', HwannR_unit='eV')
431 : 0 : call self%close_ncfile(ncfile)
432 [ # # ]: 0 : call self%write_Amnk_w90(trim(Amnkfilename))
433 : 0 : end subroutine run_all
434 : :
435 : : !subroutine remove_phase(self, psip)
436 : : ! class(WannierBuilder_t), intent(inout):: self
437 : : ! complex(dp), intent(in):: psip(:, :, :) ! (ibasis, iband, ikpt)
438 : : ! !complex(dp), intent(out):: psi(:,:,:) ! (ibasis, iband, ikpt)
439 : : ! integer:: ikpt, ibasis
440 : : ! complex(dp):: phase
441 : : ! do ikpt = 1, self%nkpt
442 : : ! do ibasis = 1, self%nbasis
443 : : ! phase = exp(-tpi_im*dot_product(self%kpts(:, ikpt), self%positions_red(:, ibasis)))
444 : : ! self%psi(ibasis, :, ikpt) = psip(ibasis, :, ikpt)*phase
445 : : ! end do
446 : : ! end do
447 : : !end subroutine remove_phase
448 : :
449 : : !===============================================================
450 : : ! Find one kpoint in a list of kpoints.
451 : : !> @
452 : : !===============================================================
453 [ - + ]: 4 : function find_kpoint(self, kpoint) result(ik)
454 : : class(WannierBuilder_t), intent(inout):: self
455 : : real(dp), intent(in):: kpoint(:)
456 : : integer:: ik, nk
457 : : integer:: i
458 : 8 : real(dp):: a(size(self%kpts, 2))
459 : 4 : nk = size(self%kpts, 2)
460 : : ! should transfer back to 1st BZ?
461 [ + + ]: 260 : do i = 1, nk
462 [ + + ]: 1028 : a(i) = sum((self%kpts(:, i) - kpoint)**2)
463 : : end do
464 : :
465 [ + - + - : 264 : ik = minloc(a, dim = 1)
- - + + -
+ ]
466 [ - + ]: 4 : if (a(ik) > 0.001) then
467 : 0 : ABI_ERROR("Error in finding kpoint from kpoint list. ")
468 : : end if
469 : 4 : end function find_kpoint
470 : :
471 : : !===============================================================
472 : : ! Calculate the weight function for each mode described by iband and ikpt
473 : : ! The
474 : : !> @
475 : : !===============================================================
476 : 194 : subroutine get_weight(self, ikpt, disentanglement, mu, sigma, weight, project_to_anchor)
477 : : class(WannierBuilder_t), intent(inout):: self
478 : : integer, intent(in):: ikpt, disentanglement
479 : : real(dp), intent(in):: mu, sigma
480 : : real(dp), intent(inout):: weight(self%nband)
481 : : logical, intent(in):: project_to_anchor
482 : :
483 : : integer:: iband, ianchor
484 : : real(dp):: proj
485 : : complex(dp):: p
486 : 194 : real(dp), pointer:: ek(:)
487 : 194 : complex(dp), pointer:: psik(:,:)
488 : :
489 : 388 : ek => self%get_evals_k(ikpt)
490 : 194 : select case (disentanglement)
491 : : case (1)
492 [ + + ]: 1040 : weight(:) = 1.0
493 : : case (2)
494 [ + + ]: 2064 : do iband = 1, self%nband
495 : 2064 : weight(iband) = fermi(ek(iband), mu, sigma)
496 : : end do
497 : : case (3)
498 [ # # ]: 0 : do iband = 1, self%nband
499 : 0 : weight(iband) = gaussian(ek(iband), mu, sigma)
500 : : end do
501 : : case default
502 [ + + - - ]: 194 : ABI_ERROR("The disentanglement function type can only be 1:unity, 2: fermi, 3: gaussian")
503 : : end select
504 : :
505 : : ! weight_mk *=\sum_anchor < psi_anchor | psi mk>
506 : :
507 [ + + ]: 194 : if( project_to_anchor) then
508 [ + - ]: 130 : if (size(self%anchor_ibands) /= 0) then
509 : 130 : psik=>self%get_psi_k(ikpt)
510 [ + + ]: 2080 : do iband = 1, self%nband
511 : 1950 : proj = 0.0_dp
512 [ + + ]: 7800 : do ianchor = 1, size(self%anchor_ibands)
513 [ + + ]: 93600 : p = dot_product(self%projectors(:, ianchor), psik(:, iband))
514 : 7800 : proj = proj+real(conjg(p)*p)
515 : : end do
516 : 2080 : weight(iband) = weight(iband)*proj
517 : : end do
518 : : end if
519 : : end if
520 : 194 : end subroutine get_weight
521 : :
522 [ - + ]: 2 : subroutine get_columns(self, psi_dagger, cols)
523 : : class(WannierBuilder_t), intent(inout):: self
524 : : complex(dp), intent(in):: psi_dagger(:, :)
525 : 4 : integer:: piv(size(psi_dagger, 2))
526 : : integer, intent(inout):: cols(self%nwann)
527 : 2 : call complex_QRCP_piv_only(psi_dagger, piv)
528 [ + + ]: 8 : cols = piv(:self%nwann)
529 : 2 : end subroutine get_columns
530 : :
531 : 192 : subroutine get_Amnk(self, ikpt, Amnk)
532 : : class(WannierBuilder_t), intent(inout):: self
533 : : integer, intent(in):: ikpt
534 : : complex(dp), intent(inout):: Amnk(self%nband, self%nwann)
535 [ + + ]: 192 : if(self%method == 1) then
536 : 128 : call self%get_scdm_Amnk(ikpt, Amnk)
537 [ + - ]: 64 : else if(self%method == 2) then
538 : 64 : call self%get_projected_Amnk(ikpt, Amnk)
539 : : end if
540 : 192 : end subroutine get_Amnk
541 : :
542 : 128 : subroutine get_scdm_Amnk(self, ikpt, Amnk)
543 : : class(WannierBuilder_t), intent(inout):: self
544 : : integer, intent(in):: ikpt
545 : 256 : complex(dp):: psi_dagger(self%nband, self%nbasis)
546 : : complex(dp), intent(inout):: Amnk(self%nband, self%nwann)
547 : 256 : real(dp):: weight(self%nband)
548 : 256 : complex(dp):: U(self%nband, self%nwann), VT(self%nwann, self%nwann)
549 : 256 : real(dp):: S(self%nband)
550 : : !real(dp):: weights(self%nband)
551 : : integer:: iband
552 : 128 : complex(dp), pointer:: p(:, :)
553 : :
554 : : call self%get_weight(ikpt, self%disentangle_func_type, self%mu, self%sigma, weight, &
555 : 128 : &project_to_anchor = self%project_to_anchor)
556 : :
557 : 128 : p => self%get_psi_k(ikpt)
558 [ + + ]: 2048 : do iband = 1, self%nband
559 [ + + ]: 30848 : psi_dagger(iband, :) = conjg(p(:, iband)) *weight(iband)
560 : : end do
561 : :
562 : : ! orthogonalize selected columns
563 [ + - + + : 6272 : call complex_svd(psi_dagger(:, self%cols), U, S, VT, 'S')
+ + ]
564 [ + - + + : 25344 : Amnk(:, :) = matmul(U, VT)
+ + + + +
+ + + ]
565 : 128 : end subroutine get_scdm_Amnk
566 : :
567 : :
568 : 64 : subroutine get_projected_Amnk(self, ikpt, Amnk)
569 : : class(WannierBuilder_t), intent(inout):: self
570 : : integer, intent(in):: ikpt
571 : : complex(dp), intent(inout):: Amnk(self%nband, self%nwann)
572 : 128 : complex(dp):: U(self%nband, self%nwann), VT(self%nwann, self%nwann)
573 : 128 : real(dp):: S(self%nband), weights(self%nband)
574 : : integer:: iband, iwann
575 : 64 : complex(dp), pointer:: psi(:, :)
576 : :
577 : 128 : psi => self%get_psi_k(ikpt)
578 : : call self%get_weight(ikpt, self%disentangle_func_type, &
579 : 64 : & self%mu, self%sigma, weights, project_to_anchor=.False.)
580 : :
581 : : ! <proj||psi> * weight
582 [ + + ]: 1024 : do iband = 1, self%nband
583 [ + + ]: 3904 : do iwann = 1, self%nwann
584 : : !Amnk(iband, iwann)= dot_product(conjg(self%projectors(:, iwann)), psi(:, iband )) * weights(iband)
585 : : !Amnk(iband, iwann)= dot_product(self%projectors(:, iwann), conjg(psi(:, iband ))) * weights(iband)
586 [ + + ]: 47040 : Amnk(iband, iwann)= dot_product(self%projectors(:, iwann), psi(:, iband )) * weights(iband)
587 : : end do
588 : : end do
589 : 64 : call complex_svd(Amnk, U, S, VT, 'S')
590 [ + - + + : 12672 : Amnk(:, :) = matmul(U, VT)
+ + + + +
+ + + ]
591 : 64 : end subroutine get_projected_Amnk
592 : :
593 : :
594 [ - + ]: 3 : subroutine get_wannR_and_HwannR(self, Rlist)
595 : : class(WannierBuilder_t), intent(inout):: self
596 : : integer, intent(in):: Rlist(:, :)
597 : : !complex(dp), intent(out):: HR(self%nwann, self%nwann, size(Rlist, 2))
598 : : !-- H(R)= \sum_k H(k) * exp(i2pi k.R)
599 : : integer:: ik, iR, nR
600 : : complex(dp):: factor
601 : 3 : nR = size(Rlist, 2)
602 [ - + + - : 15 : ABI_MALLOC(self%HwannR, (self%nwann, self%nwann, nR))
- + + - -
+ + - + -
+ - - + -
+ - + ]
603 [ - + + - : 15 : ABI_MALLOC(self%WannR, (self%nbasis, self%nwann, nR))
- + + - -
+ + - + -
+ - - + -
+ - + ]
604 : :
605 [ + + + + : 2499 : self%HwannR(:, :, :) = cmplx(0.0, 0.0, dp)
+ + ]
606 [ + + + + : 9411 : self%wannR(:, :, :) = cmplx(0.0, 0.0, dp)
+ + ]
607 [ + + ]: 195 : do ik = 1, self%nkpt
608 [ + + ]: 12483 : do iR = 1, nR
609 [ + + ]: 49152 : factor = exp(-tpi_Im*dot_product(self%kpts(:, ik), Rlist(:, iR))) * self%kweights(ik)
610 [ + + + + ]: 159744 : self%HwannR(:, :, iR) = self%HwannR(:,:, iR) + self%Hwannk(:, :, ik)*factor
611 [ + + + + ]: 602304 : self%wannR(:, :, iR) = self%wannR(:,:, iR) + self%psi_wann_k(:, :, ik)*factor
612 : : end do
613 : : end do
614 : 3 : end subroutine get_wannR_and_HwannR
615 : :
616 : 192 : subroutine Amn_to_H_from_evals(Amn, evals, nwann, nband, Hwann)
617 : : integer, intent(in):: nwann, nband
618 : : complex(dp), intent(in):: Amn(nband, nwann)
619 : : real(dp), intent(in):: evals(nband)
620 : : complex(dp), intent(inout):: Hwann(nwann, nwann)
621 : : integer:: i
622 : 384 : complex(dp):: tmp(nwann, nband)
623 : :
624 : : ! A\dagger E
625 [ + + ]: 3072 : do i = 1, nband
626 [ + + ]: 11712 : tmp(:, i) = conjg(Amn(i, :))*evals(i)
627 : : end do
628 : : ! Hwann = A\dagger @ E @ A
629 [ + - + + : 38016 : Hwann = matmul(tmp, Amn)
+ + + + +
+ + + ]
630 : 192 : end subroutine Amn_to_H_from_evals
631 : :
632 [ # # # # : 0 : subroutine Amn_to_H(Amn, psi, H0, nbasis, nwann, Hwann)
# # # # ]
633 : : ! Hwann = psi\dagger A
634 : : complex(dp), intent(in):: Amn(:, :), psi(:, :), H0(:, :)
635 : : complex(dp), intent(inout):: Hwann(:, :)
636 : : integer, intent(in):: nbasis, nwann
637 : 0 : complex(dp):: tmp(nbasis, nwann)
638 : : !Hwann = A_dagger@psi_dagger@H0@psi@A
639 [ # # # # : 0 : tmp(:, :) = matmul(psi, Amn)
# # # # #
# # # ]
640 [ # # # # : 0 : Hwann = matmul(matmul(transpose(conjg(tmp)), H0), tmp)
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# ]
641 : 0 : end subroutine Amn_to_H
642 : :
643 : 0 : subroutine write_Amnk_w90(self, fname)
644 : : ! write to Amnk file
645 : : class(WannierBuilder_t), intent(inout):: self
646 : : character(len=*), intent(in):: fname
647 : : integer:: iwann, iband, ikpt, locibnd
648 : : integer:: iun_amn
649 : : character(len=500):: msg
650 : :
651 [ # # ]: 0 : if (open_file(trim(fname)//".amn", msg, newunit=iun_amn, &
652 : : & form="formatted", status="unknown", action="write") /= 0) then
653 : 0 : ABI_ERROR(msg)
654 : : end if
655 : :
656 : :
657 : : !IF (wan_mode=='standalone') THEN
658 : : ! iun_amn = find_free_unit()
659 : : ! IF (ionode) OPEN (unit = iun_amn, file = trim(seedname)//".amn",form='formatted')
660 : : !ENDIF
661 : :
662 : : !TODO: re-enable this.
663 : : !WRITE(stdout, '(a, i8)') ' AMN: iknum = ',iknum
664 : : !
665 : : !IF (wan_mode=='standalone') THEN
666 : : ! CALL date_and_tim( cdate, ctime )
667 : : ! header='Created on '//cdate//' at '//ctime//' with SCDM '
668 : : ! IF (ionode) THEN
669 : : ! WRITE (iun_amn, *) header
670 : : ! WRITE (iun_amn, '(3i8, xxx, 2f10.6)') numbands, iknum, n_wannier, scdm_mu, scdm_sigma
671 : : ! ENDIF
672 : : !ENDIF
673 : :
674 [ # # ]: 0 : do ikpt = 1, self%nkpt
675 [ # # ]: 0 : do iwann = 1, self%nwann
676 : 0 : locibnd = 0
677 [ # # ]: 0 : do iband = 1, self%nband
678 : : !IF (excluded_band(iband)) CYCLE
679 : 0 : locibnd = locibnd+1
680 : 0 : WRITE (iun_amn, '(3i5, 2f18.12)') locibnd, iwann, ikpt, &
681 : 0 : & REAL(self%Amnk(locibnd, iwann, ikpt)), &
682 : 0 : & AIMAG(self%Amnk(locibnd, iwann, ikpt))
683 : : end do
684 : : end do
685 : : end do
686 : 0 : close (iun_amn)
687 : 0 : end subroutine write_Amnk_w90
688 : :
689 [ # # # # ]: 0 : subroutine write_Hwann_w90(self, HR, Rlist, fname)
690 : : class(WannierBuilder_t), intent(inout):: self
691 : : complex(dp), intent(in):: HR(:, :, :)
692 : : integer, intent(in):: Rlist(:, :)
693 : : character(len=*), intent(in):: fname
694 : : integer:: iR, ifile, iwann1, iwann2
695 : : character(len=500):: msg
696 : :
697 [ # # ]: 0 : if (open_file(trim(fname)//".hr", msg, newunit=ifile, &
698 : : & form="formatted", status="unknown", action="write") /= 0) then
699 : 0 : ABI_ERROR(msg)
700 : : end if
701 : :
702 : :
703 [ # # ]: 0 : do iR = 1, size(Rlist, 2)
704 : 0 : WRITE (ifile, '(3i5)') Rlist(:, iR)
705 [ # # ]: 0 : do iwann1 = 1, self%nwann
706 [ # # ]: 0 : do iwann2 = 1, self%nwann
707 : 0 : WRITE (ifile, '(f18.12)') HR(iwann1, iwann2, iR)
708 : : end do
709 : : end do
710 : : end do
711 : 0 : close (ifile)
712 : 0 : end subroutine write_Hwann_w90
713 : :
714 : 3 : subroutine create_ncfile(self, fname, ncfile)
715 : : class(WannierBuilder_t), intent(inout):: self
716 : : type(IOWannNC):: ncfile
717 : : character(len=*), intent(in):: fname
718 : 3 : ABI_UNUSED_A(self)
719 : 3 : call ncfile%initialize(filename = fname)
720 : 3 : end subroutine create_ncfile
721 : :
722 : 3 : subroutine close_ncfile(self, ncfile)
723 : : class(WannierBuilder_t), intent(inout):: self
724 : : type(IOWannNC):: ncfile
725 : 3 : ABI_UNUSED_A(self)
726 : 3 : call ncfile%close_file()
727 : 3 : end subroutine close_ncfile
728 : :
729 : :
730 : 3 : subroutine write_wann_netcdf(self, ncfile, wannR_unit, HwannR_unit)
731 : : class(WannierBuilder_t), intent(inout):: self
732 : : character(*), intent(in):: HwannR_unit, wannR_unit
733 : : type(IOWannNC), intent(inout):: ncfile
734 : : call ncfile%write_wann( nR = self%nR, ndim = self%nkdim, &
735 : : & nwann = self%nwann, nbasis = self%nbasis, Rlist = self%Rlist, &
736 : : & wannR = self%wannR, HwannR = self%HwannR, &
737 : 3 : & wannR_unit = wannR_unit, HwannR_unit = HwannR_unit)
738 : : call ncfile%write_Amnk(nkpt = self%nkpt, nband = self%nband, nwann = self%nwann, &
739 : 3 : & kpoints = self%kpts, Amnk = self%Amnk)
740 : : !TODO: add write_eigenvalues
741 : 3 : end subroutine write_wann_netcdf
742 : :
743 : 344 : subroutine get_wannier_eigen(self, kpoint, evals, evecs)
744 : : class(WannierBuilder_t), intent(inout):: self
745 : : real(dp), intent(in) :: kpoint(3)
746 : : real(dp), intent(inout) :: evals(self%nwann)
747 : : complex(dp), optional, intent(inout) :: evecs(self%nwann, self%nwann)
748 : 688 : complex(dp) :: Hk(self%nwann, self%nwann), phase
749 : 344 : type(eigensolver):: esolver
750 : : integer :: iR
751 [ + + + + ]: 4472 : Hk(:,:)=0.0_dp
752 [ + + ]: 22360 : do iR=1, self%nR
753 [ + + ]: 88064 : phase = exp(tpi_im * dot_product(kpoint, self%Rlist(:, iR)))
754 [ + + + + ]: 286552 : Hk = Hk + self%HwannR(:, :, iR) * phase
755 : : end do
756 : 344 : call esolver%run(evals, Hk)
757 : : ! Hk is overwritten as evecs
758 [ + - ]: 344 : if (present(evecs)) then
759 [ + + + + ]: 4472 : evecs=Hk
760 : : end if
761 : 344 : call esolver%finalize()
762 [ + - - + : 688 : end subroutine get_wannier_eigen
- + ]
763 : :
764 : 2 : subroutine get_wannier_eigen_klist(self, kpoints, nk, evals_nk, evecs_nk)
765 : : class(WannierBuilder_t), intent(inout):: self
766 : : integer, intent(in) :: nk
767 : : real(dp), intent(in) :: kpoints(3, nk)
768 : : real(dp), intent(inout) :: evals_nk(self%nwann, nk)
769 : : complex(dp), optional, intent(inout) :: evecs_nk(self%nwann, self%nwann, nk)
770 : : integer :: ik
771 [ + + ]: 346 : do ik =1, nk
772 [ + - ]: 346 : if (present(evecs_nk)) then
773 : : call self%get_wannier_eigen(kpoints(:, ik), &
774 : 344 : & evals_nk(:, ik), evecs_nk(:,:, ik))
775 : : else
776 : 0 : call self%get_wannier_eigen(kpoints(:, ik), evals_nk(:, ik))
777 : : end if
778 : : end do
779 [ + - ]: 2 : end subroutine get_wannier_eigen_klist
780 : :
781 : :
782 : : !============================ WannierBuilder_witheigen_t =========================
783 : :
784 [ - + - + ]: 3 : subroutine set_eigen(self, evals, psi)
785 : : class(WannierBuilder_witheigen_t) :: self
786 : : real(dp), intent(in), target:: evals(:, :) !(iband, ikpt)
787 : : complex(dp), intent(in), target:: psi(:, :, :) ! (ibasis, iband, ikpt)
788 : 3 : self%psi => psi
789 : 3 : self%evals => evals
790 : 3 : end subroutine set_eigen
791 : :
792 : 386 : function get_evals_k_from_eigen(self, ikpt) result(ek)
793 : : class(WannierBuilder_witheigen_t), intent(inout):: self
794 : : integer, intent(in):: ikpt
795 : : real(dp), pointer:: ek(:)
796 : 386 : ek => self%evals( :, ikpt)
797 : 386 : end function get_evals_k_from_eigen
798 : :
799 : :
800 : 518 : function get_psi_k_from_eigen(self, ikpt) result(psik)
801 : : class(WannierBuilder_witheigen_t), intent(inout):: self
802 : : integer, intent(in):: ikpt
803 : : complex(dp), pointer:: psik(:, :)
804 : 518 : psik => self%psi(:, :, ikpt)
805 : 518 : end function get_psi_k_from_eigen
806 : :
807 : :
808 [ # # # # : 0 : end module m_wannier_builder
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # ]
809 : : !!***
|
