TABLE OF CONTENTS
ABINIT/m_berry_curvature [ Modules ]
NAME
m_berry_curvature
FUNCTION
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG, MMignolet) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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 m_crystal 28 use m_dtset 29 use m_dtfil 30 use netcdf 31 use m_nctk 32 33 use m_time, only : cwtime, cwtime_report 34 use m_fstrings, only : strcat, sjoin, ktoa 35 use defs_datatypes, 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
m_berry_curvature/berry_curvature [ Functions ]
[ Top ] [ m_berry_curvature ] [ Functions ]
NAME
berry_curvature
FUNCTION
Computes the molecular Berry curvature and writes the reulst to a new ddb Ref: see D. Saparov PRB 105, 064303 (2022)
INPUTS
gstore<type(gstore_t)> = gstore object with el-ph matrix element computed the whole BZ dtset<type(dataset_type)> even though its listed as inout, its not modified dtfil<type(datafiles_type)> ! in_ddb<type(ddb_type)> in_ifc<type(ifc_type)> dielt(3,3)=dielectric tensor. zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement qdrp_cart(3,3,3,natom)=Quadrupole tensor on each atom in cartesian cordinates
OUTPUT
Generates a new ddb file with the molecular Berry curvature: *_BERRY_DDB
SOURCE
79 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, ebands_timrev 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 type(crystal_t) :: in_ddb_crystal 107 type(ddb_type) :: in_ddb, berry_ddb 108 type(ddb_hdr_type) :: in_ddb_hdr, berry_ddb_hdr 109 !arrays 110 integer :: units(2) 111 integer,allocatable :: blkflg(:,:,:,:) 112 integer,allocatable :: my_kqmap(:,:) 113 real(dp) :: qphon(3) 114 real(dp) :: qq_ibz(3) 115 real(dp),allocatable :: tmp_mat(:,:,:,:,:) 116 complex(dp),allocatable :: gmat(:,:,:) 117 !---------------------------------------------------------------------- 118 119 comm = gstore%comm; nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 120 units = [std_out, ab_out] 121 call cwtime(cpu_all, wall_all, gflops_all, "start") 122 123 cryst => gstore%cryst; ebands => gstore%ebands 124 natom = cryst%natom; ntypat = cryst%ntypat 125 natom3 = 3 * cryst%natom; nsppol = ebands%nsppol 126 spin_occ = one; if (nsppol == 1 .and. dtset%nspinor == 1) spin_occ = two 127 ebands_timrev = kpts_timrev_from_kptopt(ebands%kptopt) 128 129 if (my_rank == master) then 130 call wrtout(std_out, " Computing berry curvature", pre_newlines=2) 131 call gstore%print(std_out, header="Gstore", prtvol=dtset%prtvol) 132 end if 133 134 ! Consistency check 135 ierr = 0 136 ABI_CHECK_NOSTOP(gstore%qzone == "ibz", "qzone == 'ibz' is required", ierr) 137 ABI_CHECK_NOSTOP(gstore%kzone == "bz", "kzone == 'bz' is required", ierr) 138 ABI_CHECK_NOSTOP(gstore%gqk(1)%cplex == 2, "cplex == 2 is required", ierr) 139 ! Need all perts on the same procs as we have to take the outer product (ipc1, ipc2) 140 ABI_CHECK_NOSTOP(gstore%gqk(1)%pert_comm%nproc == 1, "berry_curvature is not compatible with pert_parallelism", ierr) 141 ABI_CHECK(ierr == 0, "Wrong gstore object for berry_curvature. See messages above") 142 143 ABI_CALLOC(gmat, (natom3, natom3, gstore%nqibz)) 144 145 ! Loop over collinear spins (if any) 146 do spin=1,gstore%nsppol 147 my_is = gstore%spin2my_is(spin); if (my_is == 0) cycle 148 gqk => gstore%gqk(my_is) 149 nb = gqk%nb 150 ABI_MALLOC(my_kqmap, (6, gqk%my_nk)) 151 152 ! For each q-point in the IBZ treated by me. 153 do my_iq=1,gqk%my_nq 154 print_qtime = (my_iq <= LOG_MODQ .or. mod(my_iq, LOG_MODQ) == 0) 155 if (print_qtime) call cwtime(cpu_q, wall_q, gflops_q, "start") 156 iq_ibz = gqk%my_q2ibz(1, my_iq); isym_q = gqk%my_q2ibz(2, my_iq) 157 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 tsign_q = 1; if (trev_q == 1) tsign_q = -1 160 qq_ibz = gstore%qibz(:, iq_ibz) 161 iq_glob = my_iq + gqk%my_qstart - 1 162 163 ! Find k+q in the IBZ for all my k-points. 164 if (kpts_map("symrel", ebands_timrev, cryst, gstore%krank_ibz, gqk%my_nk, gqk%my_kpts, my_kqmap, qpt=qq_ibz) /= 0) then 165 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 do my_ik=1,gqk%my_nk 170 ik_glob = my_ik + gqk%my_kstart - 1 171 172 ik_ibz = gqk%my_k2ibz(1, my_ik); isym_k = gqk%my_k2ibz(2, my_ik) 173 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 tsign_k = 1; if (trev_k == 1) tsign_k = -1 176 177 ikq_ibz = my_kqmap(1, my_ik); isym_kq = my_kqmap(2, my_ik) 178 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 do ib2=1,nb 191 band2 = ib2 + gqk%bstart - 1 192 e_b2_k = ebands%eig(band2, ik_ibz, spin) 193 f_b2_k = ebands%occ(band2, ik_ibz, spin) / spin_occ 194 e_b2_kq = ebands%eig(band2, ikq_ibz, spin) 195 f_b2_kq = ebands%occ(band2, ikq_ibz, spin) / spin_occ 196 197 do ib1=1,nb 198 band1 = ib1 + gqk%bstart - 1 199 e_b1_kq = ebands%eig(band1, ikq_ibz, spin) 200 f_b1_kq = ebands%occ(band1, ikq_ibz, spin) / spin_occ 201 e_b1_k = ebands%eig(band1, ik_ibz, spin) 202 f_b1_k = ebands%occ(band1, ik_ibz, spin) / spin_occ 203 204 ! following the formula: m = b1 and mprime = b2 205 fact(1) = f_b1_k * (one - f_b2_kq) 206 fact(2) = f_b1_kq * (one - f_b2_k) 207 208 if (all(abs(fact) < tol20)) cycle 209 210 dene = e_b1_k - e_b2_kq 211 if (abs(dene) > tol12) then 212 ! the tolerance here might need some tweaking 213 fact(1) = fact(1) / dene**2 214 else 215 ! TODO: add finite delta or do Taylor series expansion of numerator + denominator? 216 fact(1) = zero 217 end if 218 219 dene = e_b2_k - e_b1_kq 220 if (abs(dene) > tol12) then 221 ! the tolerance here might need some tweaking 222 fact(2) = fact(2) / dene**2 223 else 224 ! TODO: add finite delta or do Taylor series expansion of numerator + denominator? 225 fact(2) = zero 226 end if 227 228 ! Loop over perturbations and accumulate. 229 do my_ip1=1,gqk%my_npert 230 ipc1 = gqk%my_iperts(my_ip1) 231 do my_ip2=1,gqk%my_npert 232 ipc2 = gqk%my_iperts(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 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 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 tmp = tmp1 - tmp2 244 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 if (print_qtime) then 253 write(msg,'(4x,3(a,i0),a)')"my_iq [", my_iq, "/", gqk%my_nq, "] (tot: ", gstore%nqibz, ")" 254 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 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 gmat = j_dpc * gmat / gstore%nkbz 263 if (nsppol == 1 .and. dtset%nspinor == 1) gmat = two * gmat 264 call xmpi_sum(gmat, comm, ierr) 265 call massmult_and_breaksym_cplx(cryst%natom, cryst%ntypat, cryst%typat, gstore%ifc%amu, gmat, herm_opt=0) 266 call cwtime_report(" berry_curvature:", cpu_all, wall_all, gflops_all) 267 268 if (my_rank == master) then 269 ! Print some results to ab_out for testing purposes. 270 do iq_ibz=1,gstore%nqibz 271 write(msg, "(2a,2x,2a)")ch10, ch10, "G(q) matrix for qpt: ", trim(ktoa(gstore%qibz(:,iq_ibz))) 272 call wrtout(units, msg) 273 write(msg, "(2x,4(a6,2x), 2(a12,2x))")"idir1", "ipert1", "idir2", "ipert2", "Re(gmat)", "Im(gmat)" 274 call wrtout(units, msg) 275 do ipc2=1,natom3 276 idir2 = mod(ipc2-1, 3) + 1; ipert2 = (ipc2 - idir2) / 3 + 1 277 do ipc1=1,natom3 278 idir1 = mod(ipc1-1, 3) + 1; ipert1 = (ipc1 - idir1) / 3 + 1 279 write(msg, "(2x,4(i6,2x), 2(es12.5,2x))") & 280 idir1, ipert1, idir2, ipert2, real(gmat(ipc1, ipc2, iq_ibz)), aimag(gmat(ipc1, ipc2, iq_ibz)) 281 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 call wrtout(units, sjoin("- Reading input DDB from file:", dtfil%filddbsin)) 292 call in_ddb%from_file(dtfil%filddbsin, in_ddb_hdr, in_ddb_crystal, comm, dtset%prtvol) 293 call in_ddb_crystal%free() 294 295 ! Initialize ddb header object 296 call wrtout(units, "- Initialize berry_ddb_hdr:") 297 ! dirty trick otherwise ddb_hdr_init complains about not havinig kpt weights 298 ABI_MALLOC(dtset%wtk, (gstore%nkbz)) 299 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=in_ddb_hdr%nkpt, kpt=in_ddb_hdr%kpt, & 303 occ=in_ddb_hdr%occ) 304 ABI_FREE(dtset%wtk) 305 call in_ddb_hdr%free() 306 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 nblock = gstore%nqibz 313 mpert = natom+6 ! in principle we only need natom, but some portions of anaddb 314 ! require at least natom+6 315 msize = 3*mpert * 3*mpert 316 call wrtout(units, "- Initialize berry_ddb:") 317 call berry_ddb%init(dtset, nblock, mpert=mpert, with_d2E=.true.) 318 319 ! Set the values for the 2nd order derivatives 320 ABI_CALLOC(blkflg, (3, mpert, 3, mpert)) 321 ABI_CALLOC(tmp_mat, (2, 3, mpert, 3, mpert)) 322 blkflg(:3,:natom,:3,:natom) = 1 ! all values have been computed 323 do iq_ibz=1,gstore%nqibz 324 tmp_mat(1,:3,:natom,:3,:natom) = reshape( real(gmat(:,:,iq_ibz)), (/3,natom,3,natom/)) 325 tmp_mat(2,:3,:natom,:3,:natom) = reshape(aimag(gmat(:,:,iq_ibz)), (/3,natom,3,natom/)) 326 qphon = gstore%qibz(:,iq_ibz) 327 call berry_ddb%set_qpt(iq_ibz, qphon) 328 call berry_ddb%set_typ(iq_ibz, BLKTYP_d2E_mbc) 329 call berry_ddb%set_d2matr(iq_ibz, tmp_mat, blkflg) 330 end do 331 332 ! Output Berry ddb 333 if (my_rank == master) then 334 berry_ddb_filepath = strcat(dtfil%filnam_ds(4), "_BERRY_DDB") 335 call wrtout(units, sjoin("- Writing DDB file with Berry curvature to: ", berry_ddb_filepath), pre_newlines=2) 336 call berry_ddb%write(berry_ddb_hdr, berry_ddb_filepath) 337 end if 338 339 ! Deallocate ddb object 340 call berry_ddb_hdr%free() 341 call berry_ddb%free() 342 ABI_FREE(gmat) 343 ABI_FREE(blkflg) 344 ABI_FREE(tmp_mat) 345 346 end subroutine berry_curvature