Line data Source code
1 : !!****m* ABINIT/m_berry_curvature
2 : !! NAME
3 : !! m_berry_curvature
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2008-2026 ABINIT group (MG, MMignolet)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_berry_curvature
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_xmpi
26 : use m_errors
27 : use netcdf
28 : use m_nctk
29 :
30 : use m_crystal, only : crystal_t
31 : use m_dtset, only : dataset_type
32 : use m_dtfil, only : datafiles_type
33 : use m_time, only : cwtime, cwtime_report
34 : use m_fstrings, only : strcat, sjoin, ktoa, itoa
35 : use m_ebands, only : ebands_t
36 : use m_kpts, only : kpts_timrev_from_kptopt, kpts_map, smpbz
37 : use m_ddb_hdr, only : ddb_hdr_type, BLKTYP_d2E_mbc
38 : use m_ddb, only : ddb_type
39 : use m_ifc, only : ifc_type
40 : use m_gstore, only : gstore_t, gqk_t
41 : use m_dynmat, only : massmult_and_breaksym_cplx
42 :
43 : implicit none
44 :
45 : private
46 :
47 : public :: berry_curvature
48 : !!***
49 :
50 : contains
51 : !!***
52 :
53 : !----------------------------------------------------------------------
54 :
55 : !!****f* m_berry_curvature/berry_curvature
56 : !! NAME
57 : !! berry_curvature
58 : !!
59 : !! FUNCTION
60 : !! Computes the molecular Berry curvature and writes the reulst to a new ddb
61 : !! Ref: see D. Saparov PRB 105, 064303 (2022)
62 : !!
63 : !! INPUTS
64 : !! gstore<type(gstore_t)> = gstore object with el-ph matrix element computed
65 : !! the whole BZ
66 : !! dtset<type(dataset_type)> even though its listed as inout, its not modified
67 : !! dtfil<type(datafiles_type)>
68 : !! ! in_ddb<type(ddb_type)>
69 : !! in_ifc<type(ifc_type)>
70 : !! dielt(3,3)=dielectric tensor.
71 : !! zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
72 : !! qdrp_cart(3,3,3,natom)=Quadrupole tensor on each atom in cartesian cordinates
73 : !!
74 : !! OUTPUT
75 : !! Generates a new ddb file with the molecular Berry curvature: *_BERRY_DDB
76 : !!
77 : !! SOURCE
78 :
79 0 : subroutine berry_curvature(gstore, dtset, dtfil)
80 :
81 : !Arguments ------------------------------------
82 : !scalars
83 : type(gstore_t),target,intent(inout) :: gstore
84 : type(dataset_type),intent(inout) :: dtset
85 : type(datafiles_type),intent(in) :: dtfil
86 :
87 : !Local variables-------------------------------
88 : !scalars
89 : integer,parameter :: master = 0, LOG_MODQ = 5
90 : integer :: nproc, my_rank, ierr, comm, mpert, msize
91 : integer :: my_is, spin, nsppol, ntypat, natom, natom3, ib1, ib2, band1, band2, nb
92 : integer :: ik_ibz, isym_k, trev_k, tsign_k, g0_k(3)
93 : integer :: ikq_ibz, isym_kq, trev_kq, tsign_kq, g0_kq(3)
94 : integer :: iq_ibz, isym_q, trev_q, tsign_q, g0_q(3)
95 : integer :: my_iq, iq_glob, my_ik, ik_glob
96 : integer :: my_ip1, my_ip2, ipc1, ipc2, ipert1, ipert2, nblock, idir1, idir2
97 : logical :: isirr_k, isirr_q, isirr_kq, print_qtime
98 : real(dp) :: e_b1_k, e_b2_k, e_b1_kq, e_b2_kq, f_b1_k, f_b1_kq, f_b2_k, f_b2_kq, dene, spin_occ, fact(2)
99 : real(dp) :: cpu_all, wall_all, gflops_all, cpu_q, wall_q, gflops_q
100 : complex(dp) :: tmp, tmp1, tmp2
101 : character(len=5000) :: msg
102 : character(len=fnlen) :: berry_ddb_filepath
103 : class(crystal_t),pointer :: cryst
104 : class(ebands_t),pointer :: ebands
105 : type(gqk_t),pointer :: gqk
106 0 : type(crystal_t) :: in_ddb_crystal
107 0 : type(ddb_type) :: in_ddb, berry_ddb
108 0 : type(ddb_hdr_type) :: in_ddb_hdr, berry_ddb_hdr
109 : !arrays
110 : integer :: units(2)
111 0 : integer,allocatable :: blkflg(:,:,:,:)
112 0 : integer,allocatable :: my_kqmap(:,:)
113 : real(dp) :: qphon(3)
114 : real(dp) :: qq_ibz(3)
115 0 : real(dp),allocatable :: tmp_mat(:,:,:,:,:)
116 0 : complex(dp),allocatable :: gmat(:,:,:)
117 : !----------------------------------------------------------------------
118 :
119 0 : comm = gstore%comm; nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
120 0 : units = [std_out, ab_out]
121 0 : call cwtime(cpu_all, wall_all, gflops_all, "start")
122 :
123 0 : cryst => gstore%cryst; ebands => gstore%ebands
124 0 : natom = cryst%natom; ntypat = cryst%ntypat
125 0 : natom3 = 3 * cryst%natom; nsppol = ebands%nsppol
126 0 : spin_occ = one; if (nsppol == 1 .and. dtset%nspinor == 1) spin_occ = two
127 :
128 0 : if (my_rank == master) then
129 0 : call wrtout(std_out, " Computing berry curvature", pre_newlines=2)
130 0 : call gstore%print([std_out], header="Gstore", prtvol=dtset%prtvol)
131 : end if
132 :
133 : ! Consistency check
134 0 : ierr = 0
135 0 : ABI_CHECK_NOSTOP(gstore%qzone == "ibz", "qzone == 'ibz' is required", ierr)
136 0 : ABI_CHECK_NOSTOP(gstore%kzone == "bz", "kzone == 'bz' is required", ierr)
137 0 : ABI_CHECK_NOSTOP(gstore%gqk(1)%cplex == 2, "cplex == 2 is required", ierr)
138 : ! Need all perts on the same procs as we have to take the outer product (ipc1, ipc2)
139 0 : ABI_CHECK_NOSTOP(gstore%gqk(1)%pert_comm%nproc == 1, "berry_curvature is not compatible with pert_parallelism", ierr)
140 0 : ABI_CHECK(ierr == 0, "Wrong gstore object for berry_curvature. See messages above")
141 :
142 0 : ABI_CALLOC(gmat, (natom3, natom3, gstore%nqibz))
143 :
144 : ! Loop over collinear spins (if any)
145 0 : do spin=1,gstore%nsppol
146 0 : my_is = gstore%spin2my_is(spin); if (my_is == 0) cycle
147 0 : gqk => gstore%gqk(my_is)
148 0 : ABI_CHECK_IEQ(gqk%nb_kq, gqk%nb_k, "nb_kq /= nb_k not tested")
149 0 : nb = gqk%nb_k
150 0 : ABI_MALLOC(my_kqmap, (6, gqk%my_nk))
151 :
152 : ! For each q-point in the IBZ treated by me.
153 0 : do my_iq=1,gqk%my_nq
154 0 : print_qtime = (my_iq <= LOG_MODQ .or. mod(my_iq, LOG_MODQ) == 0)
155 0 : if (print_qtime) call cwtime(cpu_q, wall_q, gflops_q, "start")
156 0 : iq_ibz = gqk%my_q2ibz(1, my_iq); isym_q = gqk%my_q2ibz(2, my_iq)
157 0 : trev_q = gqk%my_q2ibz(6, my_iq); g0_q = gqk%my_q2ibz(3:5, my_iq)
158 : isirr_q = (isym_q == 1 .and. trev_q == 0 .and. all(g0_q == 0))
159 0 : tsign_q = 1; if (trev_q == 1) tsign_q = -1
160 0 : qq_ibz = gstore%qibz(:, iq_ibz)
161 0 : iq_glob = my_iq + gqk%my_qstart - 1
162 :
163 : ! Find k+q in the IBZ for all my k-points.
164 0 : if (kpts_map("symrel", ebands%kptopt, cryst, gstore%krank_ibz, gqk%my_nk, gqk%my_kpts, my_kqmap, qpt=qq_ibz) /= 0) then
165 0 : ABI_ERROR(sjoin("Cannot map k+q to IBZ with qpt:", ktoa(qq_ibz)))
166 : end if
167 :
168 : ! Integration over my k-points in the BZ
169 0 : do my_ik=1,gqk%my_nk
170 0 : ik_glob = my_ik + gqk%my_kstart - 1
171 :
172 0 : ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik)
173 0 : trev_k = gqk%my_k2ibz(6, my_ik); g0_k = gqk%my_k2ibz(3:5, my_ik)
174 : isirr_k = (isym_k == 1 .and. trev_k == 0 .and. all(g0_k == 0))
175 0 : tsign_k = 1; if (trev_k == 1) tsign_k = -1
176 :
177 0 : ikq_ibz = my_kqmap(1, my_ik); isym_kq = my_kqmap(2, my_ik)
178 0 : trev_kq = my_kqmap(6, my_ik); g0_kq = my_kqmap(3:5, my_ik)
179 : isirr_kq = (isym_kq == 1 .and. trev_kq == 0 .and. all(g0_kq == 0))
180 : tsign_kq = 1; if (trev_kq == 1) tsign_kq = -1
181 :
182 : ! Summation over the two band indices. NB: occupancies f are rescaled in [0, 1] when nsppol 2.
183 : ! Here we accumulate:
184 : !
185 : ! + <b1,k|D_{-q,p1}|b2,k+q> <b2,k+q|D_{q,p2}|b1,k> / (e_{b1,k} - e_{b2,k+q})^2 <<< term 1
186 : ! - <b2,k|D_{-q,p1}|b1,k+q> <b1,k+q|D_{q,p2}|b2,k> / (e_{b2,k} - e_{b1,k+q})^2 <<< term 2
187 : !
188 : ! where b1 is the initial band and b2 the final band, p1, p2 are atomic perturbations in reduced coords
189 : ! and we're summing over k in the BZ at fixed q.
190 0 : do ib2=1,nb
191 0 : band2 = ib2 + gqk%bstart_k - 1
192 0 : e_b2_k = ebands%eig(band2, ik_ibz, spin)
193 0 : f_b2_k = ebands%occ(band2, ik_ibz, spin) / spin_occ
194 0 : e_b2_kq = ebands%eig(band2, ikq_ibz, spin)
195 0 : f_b2_kq = ebands%occ(band2, ikq_ibz, spin) / spin_occ
196 :
197 0 : do ib1=1,nb
198 0 : band1 = ib1 + gqk%bstart_k - 1
199 0 : e_b1_kq = ebands%eig(band1, ikq_ibz, spin)
200 0 : f_b1_kq = ebands%occ(band1, ikq_ibz, spin) / spin_occ
201 0 : e_b1_k = ebands%eig(band1, ik_ibz, spin)
202 0 : f_b1_k = ebands%occ(band1, ik_ibz, spin) / spin_occ
203 :
204 : ! following the formula: m = b1 and mprime = b2
205 0 : fact(1) = f_b1_k * (one - f_b2_kq)
206 0 : fact(2) = f_b1_kq * (one - f_b2_k)
207 :
208 0 : if (all(abs(fact) < tol20)) cycle
209 :
210 0 : dene = e_b1_k - e_b2_kq
211 0 : if (abs(dene) > tol12) then
212 : ! the tolerance here might need some tweaking
213 0 : fact(1) = fact(1) / dene**2
214 : else
215 : ! TODO: add finite delta or do Taylor series expansion of numerator + denominator?
216 0 : fact(1) = zero
217 : end if
218 :
219 0 : dene = e_b2_k - e_b1_kq
220 0 : if (abs(dene) > tol12) then
221 : ! the tolerance here might need some tweaking
222 0 : fact(2) = fact(2) / dene**2
223 : else
224 : ! TODO: add finite delta or do Taylor series expansion of numerator + denominator?
225 0 : fact(2) = zero
226 : end if
227 :
228 : ! Loop over perturbations and accumulate.
229 0 : do my_ip1=1,gqk%my_npert
230 0 : ipc1 = gqk%my_pertcases(my_ip1)
231 0 : do my_ip2=1,gqk%my_npert
232 0 : ipc2 = gqk%my_pertcases(my_ip2)
233 : ! my_g(my_npert, nb, my_nq, nb, my_nk)
234 :
235 : ! 1st term
236 : ! <k, b1| D_{-q,p1}H |k+q, b2> * <k+q, b2| D_{q,p2}H |k, b1>
237 0 : tmp1 = fact(1) * conjg(gqk%my_g(my_ip1,ib2,my_iq,ib1,my_ik)) * gqk%my_g(my_ip2,ib2,my_iq,ib1,my_ik)
238 :
239 : ! 2nd term
240 : ! <k, b2| D_{-q,p1}H |k+q, b1> * <k+q, b1| D_{q,p2}H |k, b2>
241 0 : tmp2 = fact(2) * conjg(gqk%my_g(my_ip1,ib1,my_iq,ib2,my_ik)) * gqk%my_g(my_ip2,ib1,my_iq,ib2,my_ik)
242 :
243 0 : tmp = tmp1 - tmp2
244 0 : gmat(ipc1, ipc2, iq_ibz) = gmat(ipc1, ipc2, iq_ibz) + tmp
245 : end do
246 : end do
247 :
248 : end do ! ib1
249 : end do ! ib2
250 : end do ! my_ik
251 :
252 0 : if (print_qtime) then
253 0 : write(msg,'(4x,3(a,i0),a)')"my_iq [", my_iq, "/", gqk%my_nq, "] (tot: ", gstore%nqibz, ")"
254 0 : call cwtime_report(msg, cpu_q, wall_q, gflops_q); if (my_iq == LOG_MODQ) call wrtout(std_out, " ...")
255 : end if
256 : end do ! my_iq
257 0 : ABI_FREE(my_kqmap)
258 : end do ! spin
259 :
260 : ! Here we ALL_REDUCE all partial integrals (sum over MPI-distributed dims i.e. spins and k-points in BZ).
261 : ! Also, account for spin degeneracy as f in [0,1] if collinear.
262 0 : gmat = j_dpc * gmat / gstore%nkbz
263 0 : if (nsppol == 1 .and. dtset%nspinor == 1) gmat = two * gmat
264 0 : call xmpi_sum(gmat, comm, ierr)
265 0 : call massmult_and_breaksym_cplx(cryst%natom, cryst%ntypat, cryst%typat, gstore%ifc%amu, gmat, herm_opt=0)
266 0 : call cwtime_report(" berry_curvature:", cpu_all, wall_all, gflops_all)
267 :
268 0 : if (my_rank == master) then
269 : ! Print some results to ab_out for testing purposes.
270 0 : do iq_ibz=1,gstore%nqibz
271 0 : write(msg, "(2a,2x,2a)")ch10, ch10, "G(q) matrix for qpt: ", trim(ktoa(gstore%qibz(:,iq_ibz)))
272 0 : call wrtout(units, msg)
273 0 : write(msg, "(2x,4(a6,2x), 2(a12,2x))")"idir1", "ipert1", "idir2", "ipert2", "Re(gmat)", "Im(gmat)"
274 0 : call wrtout(units, msg)
275 0 : do ipc2=1,natom3
276 0 : idir2 = mod(ipc2-1, 3) + 1; ipert2 = (ipc2 - idir2) / 3 + 1
277 0 : do ipc1=1,natom3
278 0 : idir1 = mod(ipc1-1, 3) + 1; ipert1 = (ipc1 - idir1) / 3 + 1
279 : write(msg, "(2x,4(i6,2x), 2(es12.5,2x))") &
280 0 : idir1, ipert1, idir2, ipert2, real(gmat(ipc1, ipc2, iq_ibz)), aimag(gmat(ipc1, ipc2, iq_ibz))
281 0 : call wrtout(units, msg)
282 : end do
283 : end do
284 : end do ! iq_ibz
285 : end if
286 : !stop
287 :
288 : ! Symmetry properties:
289 : ! 1) G(-q) = G(q)^*
290 : ! 2) In the presence of time reversal symmetry, the Berry curvature would be zero
291 0 : call wrtout(units, sjoin("- Reading input DDB from file:", dtfil%filddbsin))
292 0 : call in_ddb%from_file(dtfil%filddbsin, in_ddb_hdr, in_ddb_crystal, comm, dtset%prtvol)
293 0 : call in_ddb_crystal%free()
294 :
295 : ! Initialize ddb header object
296 0 : call wrtout(units, "- Initialize berry_ddb_hdr:")
297 : ! dirty trick otherwise ddb_hdr_init complains about not havinig kpt weights
298 0 : ABI_MALLOC(dtset%wtk, (gstore%nkbz))
299 0 : dtset%wtk(:) = one / (one * gstore%nkbz)
300 : call berry_ddb_hdr%init(dtset,in_ddb_hdr%psps, in_ddb_hdr%pawtab, &
301 : dscrpt=' Molecular Berry curvature ', &
302 : nblok=gstore%nqibz, nkpt=gstore%ebands%nkpt, kpt=gstore%ebands%kptns, &
303 0 : occ=gstore%ebands%occ)
304 0 : ABI_FREE(dtset%wtk)
305 0 : call in_ddb_hdr%free()
306 0 : call in_ddb%free()
307 :
308 : ! Initialize ddb object
309 : ! nblock = number of block -> number of qpt here
310 : ! mpert = maximum number of perturbations (atom displacements + electric field + ...)
311 : ! msize = maximum size of one block of the ddb e.g. 3*mpert * 3*mpert.
312 0 : nblock = gstore%nqibz
313 0 : mpert = natom+6 ! in principle we only need natom, but some portions of anaddb
314 : ! require at least natom+6
315 0 : msize = 3*mpert * 3*mpert
316 0 : call wrtout(units, "- Initialize berry_ddb:")
317 0 : call berry_ddb%init(dtset, nblock, mpert=mpert, with_d2E=.true.)
318 :
319 : ! Set the values for the 2nd order derivatives
320 0 : ABI_CALLOC(blkflg, (3, mpert, 3, mpert))
321 0 : ABI_CALLOC(tmp_mat, (2, 3, mpert, 3, mpert))
322 0 : blkflg(:3,:natom,:3,:natom) = 1 ! all values have been computed
323 0 : do iq_ibz=1,gstore%nqibz
324 0 : tmp_mat(1,:3,:natom,:3,:natom) = reshape( real(gmat(:,:,iq_ibz)), [3,natom,3,natom])
325 0 : tmp_mat(2,:3,:natom,:3,:natom) = reshape(aimag(gmat(:,:,iq_ibz)), [3,natom,3,natom])
326 0 : qphon = gstore%qibz(:,iq_ibz)
327 0 : call berry_ddb%set_qpt(iq_ibz, qphon)
328 0 : call berry_ddb%set_typ(iq_ibz, BLKTYP_d2E_mbc)
329 0 : call berry_ddb%set_d2matr(iq_ibz, tmp_mat, blkflg)
330 : end do
331 :
332 : ! Output Berry ddb
333 0 : if (my_rank == master) then
334 0 : berry_ddb_filepath = strcat(dtfil%filnam_ds(4), "_BERRY_DDB")
335 0 : call wrtout(units, sjoin("- Writing DDB file with Berry curvature to: ", berry_ddb_filepath), pre_newlines=2)
336 0 : call berry_ddb%write(berry_ddb_hdr, berry_ddb_filepath)
337 : end if
338 :
339 : ! Deallocate ddb object
340 0 : call berry_ddb_hdr%free(); call berry_ddb%free()
341 0 : ABI_FREE(gmat)
342 0 : ABI_FREE(blkflg)
343 0 : ABI_FREE(tmp_mat)
344 :
345 0 : end subroutine berry_curvature
346 : !!***
347 :
348 : end module m_berry_curvature
349 : !!***
|
