TABLE OF CONTENTS
ABINIT/make_grad_berry [ Functions ]
NAME
make_grad_berry
FUNCTION
compute gradient contribution from berry phase in finite electric field case
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
cg(2,mcg)=input wavefunctions cgq(2,mcgq) = wavefunctions at neighboring k points cprj_k(natom,nband_k*usepaw)=cprj at this k point dimlmn(natom)=lmn_size for each atom in input order dimlmn_srt(natom)=lmn_size for each atom sorted by type direc(2,npw*nspinor)=gradient vector gs_hamk <type(gs_hamiltonian_type)>=all data for the Hamiltonian at k iband=index of band currently being treated icg=shift to be applied on the location of data in the array cg ikpt=number of the k-point currently being treated isppol=spin polarization currently treated natom=number of atoms in cell. mband =maximum number of bands mpw=maximum dimensioned size of npw mcg=second dimension of the cg array mcgq=second dimension of the cgq array mkgq = second dimension of pwnsfacq nkpt=number of k points mpi_enreg=information about MPI parallelization npw=number of planewaves in basis sphere at given k. nspinor=number of spinorial components of the wavefunctions (on current proc) nsppol=number of spin polarizations pwind(pwind_alloc,2,3) = array used to compute the overlap matrix smat between k-points (see initberry.f) pwind_alloc = first dimension of pwind pwnsfac(2,pwind_alloc) = phase factors for non-symmorphic translations (see initberry.f) pwnsfacq(2,mkgq) = phase factors for the nearest neighbours of the current k-point (electric field, MPI //)
OUTPUT
grad_berry(2,npw*nspinor) :: contribution to gradient in finite electric field case
SIDE EFFECTS
dtefield <type(efield_type)> = variables related to Berry phase calculations (see initberry.f)
NOTES
PARENTS
cgwf
CHILDREN
nonlop,pawcprj_alloc,pawcprj_copy,pawcprj_free,pawcprj_get pawcprj_symkn,smatrix,smatrix_k_paw
SOURCE
65 #if defined HAVE_CONFIG_H 66 #include "config.h" 67 #endif 68 69 #include "abi_common.h" 70 71 subroutine make_grad_berry(cg,cgq,cprj_k,detovc,dimlmn,dimlmn_srt,direc,dtefield,grad_berry,& 72 & gs_hamk,iband,icg,ikpt,isppol,mband,mcg,mcgq,mkgq,mpi_enreg,mpw,natom,nkpt,& 73 & npw,npwarr,nspinor,nsppol,pwind,pwind_alloc,pwnsfac,pwnsfacq) 74 75 use defs_abitypes 76 use defs_basis 77 use m_profiling_abi 78 use m_errors 79 use m_xmpi 80 use m_efield 81 82 use m_pawcprj, only : pawcprj_type, pawcprj_get, pawcprj_alloc, pawcprj_free, pawcprj_copy, pawcprj_symkn 83 use m_hamiltonian, only : gs_hamiltonian_type 84 85 !This section has been created automatically by the script Abilint (TD). 86 !Do not modify the following lines by hand. 87 #undef ABI_FUNC 88 #define ABI_FUNC 'make_grad_berry' 89 use interfaces_32_util 90 use interfaces_65_paw 91 use interfaces_66_nonlocal 92 !End of the abilint section 93 94 implicit none 95 96 !Arguments ------------------------------------ 97 98 !scalars 99 integer,intent(in) :: iband,icg,ikpt,isppol,mband,mcg,mcgq 100 integer,intent(in) :: mkgq,mpw,natom,nkpt,npw,nspinor,nsppol,pwind_alloc 101 type(gs_hamiltonian_type),intent(in) :: gs_hamk 102 type(efield_type),intent(inout) :: dtefield 103 type(MPI_type),intent(in) :: mpi_enreg 104 105 !arrays 106 integer,intent(in) :: dimlmn(natom),dimlmn_srt(natom) 107 integer,intent(in) :: npwarr(nkpt),pwind(pwind_alloc,2,3) 108 real(dp),intent(in) :: cg(2,mcg),cgq(2,mcgq) 109 real(dp),intent(inout) :: direc(2,npw*nspinor) 110 real(dp),intent(in) :: pwnsfac(2,pwind_alloc),pwnsfacq(2,mkgq) 111 real(dp),intent(out) :: detovc(2,2,3),grad_berry(2,npw*nspinor) 112 type(pawcprj_type),intent(in) :: cprj_k(natom,dtefield%mband_occ*gs_hamk%usepaw*dtefield%nspinor) 113 114 !Local variables------------------------------- 115 !scalars 116 integer :: choice,cpopt,ddkflag,dimenlc1,dimenlr1,iatom,icg1,icp2,idum1 117 integer :: idir,ifor,ikgf,ikptf,ikpt2,ikpt2f,ipw,i_paw_band,ispinor,itrs,itypat,job 118 integer :: klmn,mcg1_k,mcg_q,nbo,npw_k2,nspinortot,paw_opt,shiftbd,signs 119 real(dp) :: fac 120 character(len=500) :: message 121 !arrays 122 integer :: pwind_k(npw),sflag_k(dtefield%mband_occ) 123 real(dp) :: cg1_k(2,npw*nspinor),dtm_k(2),pwnsfac_k(4,mpw) 124 real(dp) :: smat_k(2,dtefield%mband_occ,dtefield%mband_occ) 125 real(dp) :: smat_inv(2,dtefield%mband_occ,dtefield%mband_occ),svectout_dum(2,0) 126 real(dp) :: dummy_enlout(0) 127 real(dp),allocatable :: cgq_k(:,:),enl_rij(:,:,:),grad_berry_ev(:,:) 128 real(dp),allocatable :: qijbkk(:,:,:),smat_k_paw(:,:,:) 129 ! type(pawcprj_type) :: cprj_dum(1,1) ! was used in on-site dipole, now suppressed 130 ! 15 June 2012 J Zwanziger 131 type(pawcprj_type),allocatable :: cprj_kb(:,:),cprj_band_srt(:,:) 132 type(pawcprj_type),allocatable :: cprj_fkn(:,:),cprj_ikn(:,:) 133 134 135 ! ********************************************************************* 136 137 !DBG_ENTER("COLL") 138 139 nbo = dtefield%mband_occ 140 141 !allocations 142 143 !Electric field: compute the gradient of the Berry phase part of the energy functional. 144 !See PRL 89, 117602 (2002), grad_berry(:,:) is the second term of Eq. (4) 145 grad_berry(:,:) = zero 146 job = 11 ; shiftbd = 1 147 mcg_q = mpw*mband*nspinor 148 mcg1_k = npw*nspinor 149 150 if (gs_hamk%usepaw /= 0) then 151 dimenlr1 = gs_hamk%lmnmax*(gs_hamk%lmnmax+1)/2 152 dimenlc1 = 2*dimenlr1 153 ABI_ALLOCATE(qijbkk,(dimenlc1,natom,nspinor**2)) 154 ABI_ALLOCATE(enl_rij,(nspinor*dimenlr1,natom,nspinor**2)) 155 ABI_ALLOCATE(smat_k_paw,(2,nbo,nbo)) 156 ABI_ALLOCATE(grad_berry_ev,(2,npw*nspinor)) 157 enl_rij = zero 158 qijbkk = zero 159 smat_k_paw = zero 160 ABI_DATATYPE_ALLOCATE(cprj_kb,(natom,nbo*nspinor)) 161 call pawcprj_alloc(cprj_kb,0,dimlmn) 162 ABI_DATATYPE_ALLOCATE(cprj_band_srt,(natom,nspinor)) 163 call pawcprj_alloc(cprj_band_srt,0,dimlmn_srt) 164 if (nkpt /= dtefield%fnkpt) then 165 ABI_DATATYPE_ALLOCATE(cprj_fkn,(natom,nbo*nspinor)) 166 ABI_DATATYPE_ALLOCATE(cprj_ikn,(natom,nbo*nspinor)) 167 call pawcprj_alloc(cprj_fkn,0,dimlmn) 168 call pawcprj_alloc(cprj_ikn,0,dimlmn) 169 else 170 ABI_DATATYPE_ALLOCATE(cprj_fkn,(0,0)) 171 ABI_DATATYPE_ALLOCATE(cprj_ikn,(0,0)) 172 end if 173 else 174 ABI_ALLOCATE(qijbkk,(0,0,0)) 175 ABI_ALLOCATE(enl_rij,(0,0,0)) 176 ABI_ALLOCATE(smat_k_paw,(0,0,0)) 177 ABI_ALLOCATE(grad_berry_ev,(0,0)) 178 ABI_DATATYPE_ALLOCATE(cprj_kb,(0,0)) 179 ABI_DATATYPE_ALLOCATE(cprj_band_srt,(0,0)) 180 ABI_DATATYPE_ALLOCATE(cprj_fkn,(0,0)) 181 ABI_DATATYPE_ALLOCATE(cprj_ikn,(0,0)) 182 end if 183 184 ikptf = dtefield%i2fbz(ikpt) 185 ikgf = dtefield%fkgindex(ikptf) ! this is the shift for pwind 186 187 do idir = 1, 3 188 ! skip idir values for which efield_dot(idir)=0 189 if (abs(dtefield%efield_dot(idir)) < tol12) cycle 190 ! Implicitly, we use the gradient multiplied by the number of k points in the FBZ 191 fac = dtefield%efield_dot(idir)*dble(dtefield%fnkpt)/& 192 & (dble(dtefield%nstr(idir))*four_pi) 193 do ifor = 1, 2 194 ! Handle dtefield%i2fbz properly and ask whether t.r.s. is used 195 ikpt2f = dtefield%ikpt_dk(ikptf,ifor,idir) 196 if (dtefield%indkk_f2ibz(ikpt2f,6) == 1) then 197 itrs = 10 198 else 199 itrs = 0 200 end if 201 ikpt2 = dtefield%indkk_f2ibz(ikpt2f,1) 202 npw_k2 = npwarr(ikpt2) 203 ABI_ALLOCATE(cgq_k,(2,nbo*nspinor*npw_k2)) 204 pwind_k(1:npw) = pwind(ikgf+1:ikgf+npw,ifor,idir) 205 pwnsfac_k(1:2,1:npw) = pwnsfac(1:2,ikgf+1:ikgf+npw) 206 sflag_k(:) = dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir) 207 smat_k(:,:,:) = dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir) 208 if (mpi_enreg%nproc_cell > 1) then 209 icg1 = dtefield%cgqindex(2,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) 210 cgq_k(:,1:nbo*nspinor*npw_k2) = & 211 & cgq(:,icg1+1:icg1+nbo*nspinor*npw_k2) 212 idum1 = dtefield%cgqindex(3,ifor+2*(idir-1),ikpt+(isppol-1)*nkpt) 213 pwnsfac_k(3:4,1:npw_k2) = pwnsfacq(1:2,idum1+1:idum1+npw_k2) 214 else 215 icg1 = dtefield%cgindex(ikpt2,isppol) 216 cgq_k(:,1:nbo*nspinor*npw_k2) = & 217 & cg(:,icg1+1:icg1+nbo*nspinor*npw_k2) 218 idum1 = dtefield%fkgindex(ikpt2f) 219 pwnsfac_k(3:4,1:npw_k2) = pwnsfac(1:2,idum1+1:idum1+npw_k2) 220 end if 221 if (gs_hamk%usepaw == 1) then 222 icp2=nbo*(ikpt2-1)*nspinor 223 call pawcprj_get(gs_hamk%atindx1,cprj_kb,dtefield%cprj,natom,1,icp2,ikpt,0,isppol,& 224 & nbo,dtefield%fnkpt,natom,nbo,nbo,nspinor,nsppol,0,& 225 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 226 if (ikpt2 /= ikpt2f) then ! construct cprj_kb by symmetry 227 call pawcprj_copy(cprj_kb,cprj_ikn) 228 call pawcprj_symkn(cprj_fkn,cprj_ikn,dtefield%atom_indsym,dimlmn,-1,gs_hamk%indlmn,& 229 & dtefield%indkk_f2ibz(ikpt2f,2),dtefield%indkk_f2ibz(ikpt2f,6),& 230 & dtefield%fkptns(:,dtefield%i2fbz(ikpt2)),& 231 & dtefield%lmax,dtefield%lmnmax,mband,natom,nbo,nspinor,& 232 & dtefield%nsym,gs_hamk%ntypat,gs_hamk%typat,dtefield%zarot) 233 call pawcprj_copy(cprj_fkn,cprj_kb) 234 end if 235 call smatrix_k_paw(cprj_k,cprj_kb,dtefield,idir,ifor,mband,natom,smat_k_paw,gs_hamk%typat) 236 end if 237 238 icg1 = 0 ; ddkflag = 1 239 call smatrix(cg,cgq_k,cg1_k,ddkflag,dtm_k,icg,icg1,itrs,& 240 & job,iband,mcg,mcg_q,mcg1_k,iband,mpw,nbo,dtefield%nband_occ(isppol),& 241 & npw,npw_k2,nspinor,pwind_k,pwnsfac_k,sflag_k,& 242 & shiftbd,smat_inv,smat_k,smat_k_paw,gs_hamk%usepaw) 243 ABI_DEALLOCATE(cgq_k) 244 detovc(:,ifor,idir) = dtm_k(:) !store the determinant of the overlap 245 if (sqrt(dtm_k(1)*dtm_k(1) + dtm_k(2)*dtm_k(2)) < tol12) then 246 write(message,'(3a,i5,a,i3,a,a,a)') & 247 & ' (electric field)',ch10,& 248 & ' For k-point #',ikpt,' and band # ',iband,',',ch10,& 249 & ' the determinant of the overlap matrix is found to be 0. Fixing...' 250 ! REC try this: 251 write(std_out,*)message,dtm_k(1:2) 252 if(abs(dtm_k(1))<=1d-12)dtm_k(1)=1d-12 253 if(abs(dtm_k(2))<=1d-12)dtm_k(2)=1d-12 254 write(std_out,*)' Changing to:',dtm_k(1:2) 255 ! REC MSG_BUG(message) 256 end if 257 258 if (gs_hamk%usepaw == 1) then 259 ! this loop applies discretized derivative of projectors 260 ! note that qijb_kk is sorted by input atom order, but nonlop wants it sorted by type 261 do iatom = 1, natom 262 itypat = gs_hamk%typat(gs_hamk%atindx1(iatom)) 263 do klmn = 1, dtefield%lmn2_size(itypat) 264 ! note: D_ij-like terms have 4 spinor components: 11, 22, 12, and 21. Here the qijb is diagonal 265 ! in spin space so only the first two are nonzero and they are equal 266 do ispinor = 1, nspinor 267 qijbkk(2*klmn-1,iatom,ispinor) = dtefield%qijb_kk(1,klmn,gs_hamk%atindx1(iatom),idir) 268 qijbkk(2*klmn, iatom,ispinor) = dtefield%qijb_kk(2,klmn,gs_hamk%atindx1(iatom),idir) 269 if (ifor > 1) qijbkk(2*klmn,iatom,ispinor) = -qijbkk(2*klmn,iatom,ispinor) 270 end do 271 end do ! end loop over lmn2_size 272 end do ! end loop over natom 273 274 choice = 1 275 signs = 2 276 paw_opt = 1 277 cpopt = 2 ! use cprj_kb in memory 278 nspinortot=min(2,nspinor*(1+mpi_enreg%paral_spinor)) 279 do i_paw_band = 1, nbo 280 281 call pawcprj_get(gs_hamk%atindx,cprj_band_srt,cprj_kb,natom,i_paw_band,0,ikpt,1,& 282 & isppol,nbo,1,natom,1,nbo,nspinor,nsppol,0,& 283 & mpicomm=mpi_enreg%comm_kpt,proc_distrb=mpi_enreg%proc_distrb) 284 285 ! Pass dummy_enlout to avoid aliasing (enl, enlout) 286 call nonlop(choice,cpopt,cprj_band_srt,dummy_enlout,gs_hamk,idir,(/zero/),mpi_enreg,1,0,& 287 & paw_opt,signs,svectout_dum,0,direc,grad_berry_ev,enl=qijbkk) 288 289 ! Add i*fac*smat_inv(i_paw_band,iband)*grad_berry_ev to the gradient 290 do ipw = 1, npw*nspinor 291 292 grad_berry(1,ipw) = grad_berry(1,ipw) - & 293 & fac*(smat_inv(2,i_paw_band,iband)*grad_berry_ev(1,ipw) + & 294 & smat_inv(1,i_paw_band,iband)*grad_berry_ev(2,ipw)) 295 296 grad_berry(2,ipw) = grad_berry(2,ipw) + & 297 & fac*(smat_inv(1,i_paw_band,iband)*grad_berry_ev(1,ipw) - & 298 & smat_inv(2,i_paw_band,iband)*grad_berry_ev(2,ipw)) 299 300 end do 301 end do 302 end if ! end if PAW 303 304 ! Add i*fac*cg1_k to the gradient 305 do ipw = 1, npw*nspinor 306 grad_berry(1,ipw) = grad_berry(1,ipw) - fac*cg1_k(2,ipw) 307 grad_berry(2,ipw) = grad_berry(2,ipw) + fac*cg1_k(1,ipw) 308 end do 309 fac = -1._dp*fac 310 dtefield%sflag(:,ikpt+(isppol-1)*nkpt,ifor,idir) = sflag_k(:) 311 dtefield%sflag(iband,ikpt+(isppol-1)*nkpt,ifor,idir) = 0 312 dtefield%smat(:,:,:,ikpt+(isppol-1)*nkpt,ifor,idir) = smat_k(:,:,:) 313 end do ! ifor 314 315 ! if (gs_hamk%usepaw == 1) then 316 ! ! call nonlop to apply on-site dipole <EV> part to direc 317 ! ! note that rij is sorted by input atom order, but nonlop wants it sorted by type 318 ! do iatom = 1, natom 319 ! itypat = gs_hamk%typat(gs_hamk%atindx1(iatom)) 320 ! do klmn = 1, dtefield%lmn2_size(itypat) 321 ! ! note: D_ij-like terms have 4 spinor components: 11, 22, 12, and 21. Here the enl_rij is diagonal 322 ! ! in spin space so only the first two are nonzero and they are equal 323 ! do ispinor = 1, nspinor 324 ! if (nspinor == 1) then 325 ! enl_rij(klmn,iatom,ispinor) = dtefield%rij(klmn,itypat,idir) 326 ! else 327 ! enl_rij(2*klmn-1,iatom,ispinor) = dtefield%rij(klmn,itypat,idir) 328 ! end if 329 ! end do 330 ! end do ! end loop over lmn2_size 331 ! end do ! end loop over natom 332 ! cpopt = -1 ! compute cprj inside nonlop because we do not have them for direc 333 ! call nonlop(choice,cpopt,cprj_dum,dummy_enlout,gs_hamk,idir,zero,mpi_enreg,1,0,& 334 ! & paw_opt,signs,svectout_dum,0,direc,grad_berry_ev,enl=enl_rij) 335 ! grad_berry(:,:) = grad_berry(:,:) - dtefield%efield_dot(idir)*grad_berry_ev(:,:)/two_pi 336 ! end if 337 338 end do ! idir 339 340 !deallocations 341 if(gs_hamk%usepaw /= 0) then 342 call pawcprj_free(cprj_kb) 343 call pawcprj_free(cprj_band_srt) 344 if (nkpt /= dtefield%fnkpt) then 345 call pawcprj_free(cprj_fkn) 346 call pawcprj_free(cprj_ikn) 347 end if 348 end if 349 ABI_DEALLOCATE(grad_berry_ev) 350 ABI_DEALLOCATE(qijbkk) 351 ABI_DEALLOCATE(enl_rij) 352 ABI_DEALLOCATE(smat_k_paw) 353 ABI_DATATYPE_DEALLOCATE(cprj_kb) 354 ABI_DATATYPE_DEALLOCATE(cprj_band_srt) 355 ABI_DATATYPE_DEALLOCATE(cprj_fkn) 356 ABI_DATATYPE_DEALLOCATE(cprj_ikn) 357 358 !DBG_EXIT("COLL") 359 360 end subroutine make_grad_berry