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