Line data Source code
1 : !!****m*ABINIT/m_lwf
2 : !! NAME
3 : !! m_lwf
4 : !!
5 : !! FUNCTION
6 : !! Module for the lattice Wannier function
7 : !! Container type is defined, and destruction, print subroutines
8 : !! as well as the central mkphdos
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 1999-2026 ABINIT group (HeXu)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public Licence, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! SOURCE
18 :
19 :
20 : !! Todolist:
21 : !! - LWF output filename change
22 : !! - LWF specific outputs: lwf_masses, original crystal structure
23 : !! - Units of Hamiltonian need to be changed to eV-Angstrom units
24 : !! - dipole-dipole
25 : !! - MPI
26 : !! - exclude_bands
27 :
28 : #if defined HAVE_CONFIG_H
29 : #include "config.h"
30 : #endif
31 :
32 : #include "abi_common.h"
33 :
34 : module m_lwf
35 :
36 : use, intrinsic :: iso_c_binding
37 : use defs_basis
38 : use m_errors
39 : use m_xmpi
40 : use m_abicore
41 : use m_htetra
42 : use m_numeric_tools
43 : use m_cgtools
44 : use m_crystal
45 : use m_nctk
46 : use m_atprj
47 : use m_sortph
48 : use m_ddb
49 : use netcdf
50 : use m_supercell
51 : use m_dtset
52 : use m_krank
53 :
54 : use m_fstrings, only : ltoa
55 : use m_io_tools, only : open_file
56 : use m_geometry, only : mkrdim, symredcart, normv, phdispl_cart2red
57 : use m_dynmat, only : gtdyn9, dfpt_phfrq, dfpt_prtph, pheigvec_normalize, massmult_and_breaksym, phdispl_from_eigvec
58 : use m_bz_mesh, only : make_path
59 : use m_ifc, only : ifc_type
60 : use m_anaddb_dataset, only : anaddb_dataset_type
61 : use m_kpts, only : kpts_ibz_from_kptrlatt
62 : use m_scdm_math, only : build_Rgrid
63 : use m_wannier_builder, only : WannierBuilder_witheigen_t
64 : use m_wann_netcdf, only : IOWannNC
65 :
66 :
67 : implicit none
68 :
69 : type LatticeWannier
70 : type(WannierBuilder_witheigen_t):: scdm
71 : type(ifc_type), pointer:: ifc
72 : type(crystal_t), pointer:: crystal
73 : integer:: qptrlatt(3, 3)
74 : real(dp):: shiftq(3)
75 : integer:: nqibz
76 : real(dp), allocatable:: qibz(:, :), qweights(:)
77 : integer:: nR, natom
78 : integer, allocatable:: Rlist(:, :)
79 : real(dp), allocatable:: eigenvalues(:, :) ! iband, iqpt
80 : complex(dp), allocatable:: eigenvectors(:,:, :) ! ibasis, iband, iqpt
81 : integer:: comm, nprocs, my_rank
82 : contains
83 : procedure:: initialize
84 : procedure:: finalize
85 : procedure:: init_mpi
86 : procedure:: sanity_check
87 : procedure:: prepare_qpoints
88 : procedure:: prepare_Rlist
89 : procedure:: get_ifc_eigens
90 : procedure:: run_all
91 : procedure:: write_lwf_nc
92 : procedure:: print_Rlist
93 : procedure :: write_bands
94 : end type LatticeWannier
95 :
96 : private
97 :
98 : public:: run_lattice_wannier
99 :
100 : contains
101 3 : subroutine initialize(self, ifc, crystal, dtset, comm)
102 : class(LatticeWannier), intent(inout):: self
103 : integer, intent(in):: comm
104 : type(ifc_type), target, intent(in):: ifc
105 : type(crystal_t), target, intent(in):: crystal
106 : type(anaddb_dataset_type), intent(in):: dtset
107 6 : real(dp):: symcart(3, 3, crystal%nsym)
108 : integer:: isym
109 : ! TODO: add exclude_bands
110 : integer:: exclude_bands(0)
111 :
112 3 : character(len = 500+dtset%lwf_nwann*10):: msg
113 : real(dp):: mu, sigma
114 :
115 3 : self%ifc => ifc
116 3 : self%crystal => crystal
117 :
118 : ! mpi vars
119 3 : call self%init_mpi(comm)
120 :
121 : ! sanity check
122 3 : call self%sanity_check(dtset)
123 :
124 147 : do isym = 1, crystal%nsym
125 147 : call symredcart(crystal%rprimd, crystal%gprimd, symcart(:,:,isym), crystal%symrel(:,:,isym))
126 : end do
127 :
128 : ! prepare qpoints
129 3 : call self%prepare_qpoints(crystal, dtset)
130 3 : call self%prepare_Rlist()
131 : ! prepare eigen values and eigen vectors
132 :
133 3 : call self%get_ifc_eigens(ifc, crystal)
134 :
135 : ! set mu and sigma from cm^-1 to eigenvalue
136 3 : if(dtset%lwfflag == 1 .or. dtset%lwfflag == 2) then
137 6 : mu = freq_to_eigenval(dtset%lwf_mu/Ha_cmm1)
138 3 : sigma = freq_to_eigenval(dtset%lwf_sigma/Ha_cmm1)
139 : end if
140 :
141 : ! set up scdm
142 :
143 : call self%scdm%initialize( &
144 : & kpts = self%qibz, kweights = self%qweights, Rlist = self%Rlist, &
145 : & nwann = dtset%lwf_nwann, nbasis=self%natom*3, nband=self%natom*3, disentangle_func_type = dtset%lwf_disentangle, &
146 : & mu = mu, sigma = sigma, exclude_bands = exclude_bands, &
147 3 : & project_to_anchor = (dtset%lwf_anchor_proj > 0 ), method = dtset%lwfflag)
148 3 : call self%scdm%set_eigen(evals = self%eigenvalues, psi = self%eigenvectors)
149 :
150 3 : if(dtset%lwfflag == 1) then
151 2 : write(msg, '(a)') ' Constructing LWF with SCDM-k method.'
152 6 : call wrtout([ab_out, std_out], msg )
153 2 : write(msg, '(a, i0)') ' Number of LWF: ', dtset%lwf_nwann
154 6 : call wrtout([ab_out, std_out], msg )
155 2 : write(msg, '(a, 3f8.5)') ' Anchor Points q-point: ', &
156 4 : & dtset%lwf_anchor_qpt(1), dtset%lwf_anchor_qpt(2), dtset%lwf_anchor_qpt(3)
157 : !write(msg, '(2a) ') 'Anchor points band indices: ', trim(ltoa(self%scdm%anchor_ibands))
158 6 : call wrtout([ab_out, std_out], msg )
159 :
160 2 : if(dtset%lwf_anchor_iband(1) > 0) then
161 2 : call self%scdm%set_anchor( dtset%lwf_anchor_qpt, dtset%lwf_anchor_iband)
162 : else
163 0 : call self%scdm%set_anchor(anchor_kpt = dtset%lwf_anchor_qpt)
164 : end if
165 : ! output information
166 1 : else if(dtset%lwfflag == 2) then
167 : ! TODO: projected lattice wannier function
168 1 : write(msg, '(a)') ' Constructing LWF with projected wannier function method.'
169 3 : call wrtout([ab_out, std_out], msg )
170 1 : call self%scdm%set_disp_projector(dtset%lwf_projector)
171 1 : write(msg, '(2a)') ' The projectors: ', trim(ltoa(dtset%lwf_projector))
172 3 : call wrtout([ab_out, std_out], msg )
173 : end if
174 :
175 3 : end subroutine initialize
176 :
177 3 : subroutine finalize(self)
178 : class(LatticeWannier), intent(inout):: self
179 3 : call self%scdm%finalize()
180 3 : ABI_SFREE(self%qibz)
181 3 : ABI_SFREE(self%qweights)
182 3 : ABI_SFREE(self%eigenvalues)
183 3 : ABI_SFREE(self%eigenvectors)
184 3 : ABI_SFREE(self%Rlist)
185 3 : end subroutine finalize
186 :
187 3 : subroutine init_mpi(self, comm)
188 : class(LatticeWannier), intent(inout):: self
189 : integer, intent(in):: comm
190 3 : self%comm = comm
191 3 : self%nprocs = xmpi_comm_size(comm)
192 3 : self%my_rank = xmpi_comm_rank(comm)
193 3 : end subroutine init_mpi
194 :
195 3 : subroutine sanity_check(self, dtset)
196 : class(LatticeWannier), intent(inout):: self
197 : type(anaddb_dataset_type), intent(in):: dtset
198 3 : if(self%nprocs /= 1) then
199 0 : ABI_ERROR(" MPI is not yet implemented for Lattice Wannier function.")
200 : end if
201 :
202 3 : if(dtset%dipdip > 0) then
203 0 : ABI_ERROR(" dipdip is not yet implemented for Lattice Wannier function.")
204 : end if
205 3 : end subroutine sanity_check
206 :
207 :
208 3 : subroutine prepare_qpoints(self, crystal, dtset)
209 : class(LatticeWannier), intent(inout):: self
210 : type(crystal_t), intent(in):: crystal
211 : type(anaddb_dataset_type), intent(in):: dtset
212 : integer:: nqshft = 1
213 : real(dp):: lwf_qshift(3, 1)
214 : integer:: nqbz
215 3 : real(dp), allocatable:: qbz(:, :)
216 : integer:: in_qptrlatt(3, 3), new_qptrlatt(3, 3)
217 3 : integer, allocatable:: bz2ibz_smap(:,:)!, bz2ibz(:)
218 3 : real(dp), allocatable:: new_shiftq(:,:)
219 : integer, parameter:: bcorr0 = 0, master = 0
220 : integer:: my_qptopt
221 : character(len = 500):: msg
222 : integer:: iqpt, nkpout
223 : ! Copied from m_phonons/mkphdos
224 3 : in_qptrlatt = 0
225 3 : in_qptrlatt(1, 1) = dtset%lwf_ngqpt(1)
226 3 : in_qptrlatt(2, 2) = dtset%lwf_ngqpt(2)
227 3 : in_qptrlatt(3, 3) = dtset%lwf_ngqpt(3)
228 :
229 : ! TODO: shift is now not supported.
230 3 : nqshft = 1
231 3 : lwf_qshift(:, :) = 0.0_dp
232 :
233 : ! TODO: HeXu: symmetry is not used. Check if it is available.
234 : ! In m_phonons, mkphdos, there is the comment:
235 : ! Rotate e(q) to get e(Sq) to account for symmetrical q-points in BZ.
236 : ! eigenvectors indeed are not invariant under rotation. See e.g. Eq 39-40 of PhysRevB.76.165108 [[cite:Giustino2007]].
237 : ! In principle there's a phase due to nonsymmorphic translations but we here need |e(Sq)_iatom|**2
238 3 : my_qptopt = 3
239 :
240 : call kpts_ibz_from_kptrlatt(crystal, in_qptrlatt, my_qptopt, nqshft, lwf_qshift, &
241 : & self%nqibz, self%qibz, self%qweights, nqbz, qbz, new_kptrlatt = new_qptrlatt, &
242 3 : & new_shiftk = new_shiftq, bz2ibz = bz2ibz_smap)
243 :
244 3 : ABI_FREE(bz2ibz_smap)
245 : !ABI_FREE(bz2ibz)
246 :
247 39 : self%qptrlatt = new_qptrlatt
248 12 : self%shiftq(:) = new_shiftq(:, 1) ! only one shift in output
249 3 : if (self%my_rank == master) then
250 3 : write(msg, "(3a, i0)")" LWF ngqpt: ", trim(ltoa(dtset%lwf_ngqpt)), ", qptopt: ", my_qptopt
251 9 : call wrtout([ab_out, std_out], msg)
252 3 : write(msg, "(2(a, i0))")" Number of q-points in the IBZ: ", self%nqibz, ", number of MPI processes: ", self%nprocs
253 9 : call wrtout([ab_out, std_out], msg)
254 :
255 3 : write(msg, "(a)") " List of q-points: "
256 9 : call wrtout([ab_out, std_out], msg)
257 :
258 3 : write(msg,'(a,i8)')' Grid q points for sampling in the reciprocal space : ',self%nqibz
259 9 : call wrtout([ab_out, std_out],msg,'COLL')
260 3 : nkpout=self%nqibz
261 3 : if(self%nqibz>80)then
262 0 : call wrtout([ab_out, std_out],' greater than 80, so only write 20 of them ','COLL')
263 0 : nkpout=20
264 : end if
265 195 : do iqpt=1,nkpout
266 192 : write(msg, '(1x,i2,a2,3es16.8)' )iqpt,') ',self%qibz(1,iqpt),self%qibz(2,iqpt),self%qibz(3,iqpt)
267 579 : call wrtout([ab_out, std_out], msg, 'COLL')
268 : end do
269 : end if
270 3 : ABI_SFREE(new_shiftq)
271 3 : ABI_SFREE(qbz)
272 3 : end subroutine prepare_qpoints
273 :
274 3 : subroutine print_Rlist(self, dtset)
275 : class(LatticeWannier), intent(inout):: self
276 : type(anaddb_dataset_type) :: dtset
277 : integer :: master=0
278 : integer :: nRout, ii
279 : character(len = 500):: msg
280 :
281 3 : if (self%my_rank == master) then
282 3 : write(msg, "(2a)")" LWF is transformed to the real space cells: ", trim(ltoa(dtset%lwf_ngqpt))
283 9 : call wrtout([ab_out, std_out], msg)
284 3 : write(msg, "((a, i0))")" Number of cells: ", self%nR
285 9 : call wrtout([ab_out, std_out], msg)
286 :
287 3 : write(msg,'(a,i8)')' R-vectors for cells in the real space : ',self%nR
288 9 : call wrtout([ab_out, std_out],msg,'COLL')
289 3 : nRout=self%nR
290 3 : if(self%nR>80)then
291 0 : call wrtout([ab_out, std_out],' greater than 80, so only write 20 of them ','COLL')
292 0 : nRout=20
293 : end if
294 195 : do ii=1,nRout
295 192 : write(msg, '(1x,i2,a2,3i8)' )ii,') ',self%Rlist(1,ii),self%Rlist(2,ii),self%Rlist(3,ii)
296 579 : call wrtout([ab_out, std_out], msg, 'COLL')
297 : end do
298 : end if
299 3 : end subroutine print_Rlist
300 :
301 :
302 :
303 :
304 3 : subroutine prepare_Rlist(self)
305 : class(LatticeWannier), intent(inout):: self
306 : integer:: qptrlatt(3), i
307 3 : self%nR=1
308 12 : do i = 1, 3
309 9 : qptrlatt(i) = self%qptrlatt(i, i)
310 12 : self%nR=self%nR*qptrlatt(i)
311 : end do
312 3 : call build_Rgrid(qptrlatt, self%Rlist)
313 3 : end subroutine prepare_Rlist
314 :
315 2886 : elemental function freq_to_eigenval(f) result (evalue)
316 : real(dp), intent(in):: f
317 : real(dp):: evalue
318 6 : if(f < -1d-16) then
319 144 : evalue = - f*f
320 2742 : else if (f > 1d-16) then
321 2731 : evalue = f*f
322 : else
323 : evalue = 0.0_dp
324 : end if
325 : end function freq_to_eigenval
326 :
327 3 : subroutine get_ifc_eigens(self, ifc, crystal)
328 : class(LatticeWannier), intent(inout):: self
329 : type(ifc_type), intent(in):: ifc
330 : type(crystal_t), intent(in):: crystal
331 6 : real(dp):: eigvec(2, 3, Crystal%natom, 3*Crystal%natom), phfrq(3*Crystal%natom)
332 6 : real(dp):: displ(2*3*Crystal%natom*3*Crystal%natom)
333 : integer:: iq_ibz
334 : integer:: natom, natom3
335 : integer:: iatom, iband, i3
336 : complex(dp):: phase
337 3 : natom = crystal%natom
338 3 : self%natom = natom
339 3 : natom3 = natom*3
340 12 : ABI_MALLOC(self%eigenvalues, (natom3, self%nqibz))
341 15 : ABI_MALLOC(self%eigenvectors, (natom3, natom3, self%nqibz))
342 : !print *, "trans:", ifc%trans
343 : !print *, "xred:", crystal%xred
344 195 : do iq_ibz = 1, self%nqibz
345 192 : call ifc%fourq(crystal, self%qibz(:,iq_ibz), phfrq, displ, out_eigvec = eigvec)
346 : ! frequency to eigenvalues
347 3072 : self%eigenvalues(:, iq_ibz) = freq_to_eigenval(phfrq)
348 : ! remove phases from eigenvector
349 1155 : do iatom = 1, natom
350 : ! to remove the phase factor exp(iqr)
351 : !phase = exp(-cmplx(0.0_dp, two_pi)*dot_product(crystal%xred(:, iatom), self%qibz(:, iq_ibz) ))
352 3840 : phase = exp(-cmplx(0.0_dp, two_pi)*dot_product(ifc%trans(:, iatom), self%qibz(:, iq_ibz)))
353 15552 : do iband=1, natom3
354 58560 : do i3 = 1, 3
355 : self%eigenvectors((iatom-1)*3+i3, iband, iq_ibz ) = &
356 57600 : & CMPLX(eigvec(1, i3, iatom, iband), eigvec(2, i3, iatom, iband)) * phase
357 : end do
358 : end do
359 : end do
360 : end do
361 3 : end subroutine get_ifc_eigens
362 :
363 3 : subroutine write_lwf_nc(self, prefix)
364 : class(LatticeWannier), intent(inout):: self
365 : character(len=*), intent(in) :: prefix
366 : character(len=500) :: msg
367 : type(IOWannNC):: ncfile
368 3 : call self%scdm%create_ncfile(trim(prefix)//"_lwf.nc", ncfile)
369 3 : call self%scdm%write_wann_netcdf(ncfile, wannR_unit='dimensionless', HwannR_unit='Ha')
370 : !NCF_CHECK(self%crystal%ncwrite(ncfile%ncid))
371 3 : call self%scdm%close_ncfile(ncfile)
372 3 : write(msg, '(a)') ' LWF construction finished.'
373 9 : call wrtout([ab_out, std_out], msg)
374 3 : write(msg, '(a)') ' LWF coefficients and Hamiltonian writen to file: '//trim(prefix)//"_lwf.nc ."
375 9 : call wrtout([ab_out, std_out], msg)
376 3 : end subroutine write_lwf_nc
377 :
378 :
379 3 : subroutine run_all(self, prefix, dtset)
380 : class(LatticeWannier), intent(inout):: self
381 : character(len=*), intent(in) :: prefix
382 : type(anaddb_dataset_type), intent(in):: dtset
383 3 : call self%scdm%construct_wannier()
384 3 : call self%print_Rlist(dtset)
385 3 : call self%write_lwf_nc(prefix = prefix)
386 3 : if(dtset%nqpath>0) then
387 2 : call self%write_bands(prefix, dtset)
388 : end if
389 3 : end subroutine run_all
390 :
391 3 : subroutine run_lattice_wannier(ifc, crystal, dtset, prefix, comm)
392 : integer, intent(in):: comm
393 : character(len=*), intent(in):: prefix
394 : type(ifc_type), intent(in):: ifc
395 : type(crystal_t), intent(in):: crystal
396 : type(anaddb_dataset_type), intent(in):: dtset
397 3 : type(LatticeWannier):: lwf
398 3 : call lwf%initialize(ifc, crystal, dtset, comm)
399 3 : call lwf%run_all(prefix, dtset)
400 3 : call lwf%finalize()
401 3 : end subroutine run_lattice_wannier
402 :
403 :
404 2 : subroutine write_bands(self, prefix, inp)
405 : class(LatticeWannier), intent(inout):: self
406 : character(len=*), intent(in):: prefix
407 : type(anaddb_dataset_type), intent(in):: inp
408 2 : integer,allocatable :: ndiv(:)
409 : integer :: nfineqpath
410 2 : real(dp),allocatable :: fineqpath(:,:)
411 2 : real(dp),allocatable :: phfrq(:, :)
412 2 : complex(dp), allocatable :: eigvec(:,:, :)
413 : !real(dp),allocatable :: weights(:)
414 : integer :: iq, iband
415 : real(dp) :: f
416 6 : ABI_MALLOC(ndiv,(inp%nqpath-1))
417 2 : call make_path(inp%nqpath,inp%qpath,self%Crystal%gmet,'G',inp%ndivsm,ndiv,nfineqpath,fineqpath,std_out)
418 2 : ABI_FREE(ndiv)
419 :
420 8 : ABI_MALLOC(phfrq, (inp%lwf_nwann, nfineqpath))
421 10 : ABI_MALLOC(eigvec, (inp%lwf_nwann, inp%lwf_nwann, nfineqpath))
422 :
423 2 : call self%scdm%get_wannier_eigen_klist(fineqpath, nfineqpath, phfrq, eigvec)
424 346 : do iq=1, nfineqpath
425 1378 : do iband =1, inp%lwf_nwann
426 1032 : f=phfrq(iband, iq)
427 1376 : if (f>1.0d-9) then
428 648 : phfrq(iband, iq)= sqrt(f)
429 384 : else if (f<-1.0d-9) then
430 383 : phfrq(iband, iq)= -sqrt(-f)
431 : else
432 1 : phfrq(iband, iq)= 0.0_dp
433 : end if
434 :
435 : end do
436 : end do
437 :
438 2 : call write_phfrq(trim(prefix)//"_lwf_PHFRQ", inp%lwf_nwann, nfineqpath, phfrq)
439 :
440 2 : ABI_FREE(phfrq)
441 2 : ABI_FREE(eigvec)
442 2 : ABI_FREE(fineqpath)
443 2 : end subroutine write_bands
444 :
445 2 : subroutine write_phfrq(path,nlwf,nqpts,phfreq)
446 :
447 : !Arguments ------------------------------------
448 : !scalars
449 : integer,intent(in) :: nqpts, nlwf
450 : character(len=*),intent(in) :: path
451 : !arrays
452 : !real(dp),intent(in) :: qpoints(3,nqpts)
453 : real(dp),intent(in) :: phfreq(nlwf,nqpts)
454 :
455 : !Local variables-------------------------------
456 : !scalars
457 : integer :: nphmodes, iq, iunit
458 : !real(dp) :: dummy
459 : character(len=300) :: fmt
460 : character(len=500) :: msg
461 :
462 : ! *************************************************************************
463 :
464 2 : nphmodes = nlwf
465 : !dummy = qpoints(1,1); dummy = weights(1)
466 2 : if (open_file(path, msg, newunit=iunit, form="formatted", status="unknown", action="write") /= 0) then
467 0 : ABI_ERROR(msg)
468 : end if
469 :
470 2 : write (iunit, '(a)') '# ABINIT generated LWF phonon band structure file. All in Ha atomic units'
471 2 : write (iunit, '(a)') '# '
472 2 : write (iunit, '(a,i0)') '# number_of_qpoints ', nqpts
473 2 : write (iunit, '(a,i0)') '# number_of_phonon_modes ', nphmodes
474 2 : write (iunit, '(a)') '# '
475 2 : write (fmt,'(a,i0,a)') "(I5, ", nphmodes, "E20.10)"
476 346 : do iq= 1, nqpts
477 346 : write (iunit, fmt) iq, phfreq(:,iq)
478 : end do
479 2 : close(iunit)
480 :
481 2 : end subroutine write_phfrq
482 :
483 :
484 :
485 0 : end module m_lwf
486 : !!***
|
