TABLE OF CONTENTS
- ABINIT/dterm_aij
- ABINIT/dterm_alloc
- ABINIT/dterm_BM
- ABINIT/dterm_free
- ABINIT/dterm_LR
- ABINIT/dterm_qij
- ABINIT/dterm_SOB1
- ABINIT/lamb_core
- ABINIT/local_fermie
- ABINIT/make_d
- ABINIT/make_pcg1
- ABINIT/orbmag
- ABINIT/orbmag_cc_k
- ABINIT/orbmag_nl1_k
- ABINIT/orbmag_nl_k
- ABINIT/orbmag_output
- ABINIT/orbmag_vv_k
- ABINIT/tt_me
- ABINIT/txt_me
ABINIT/dterm_aij [ Functions ]
NAME
dterm_aij
FUNCTION
transfer paw_ij to dterm%aij in more convenient format
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
PARENTS
m_orbmag
CHILDREN
SOURCE
2424 subroutine dterm_aij(atindx,dterm,dtset,paw_ij,pawtab) 2425 2426 !Arguments ------------------------------------ 2427 !scalars 2428 type(dterm_type),intent(inout) :: dterm 2429 type(dataset_type),intent(in) :: dtset 2430 2431 !arrays 2432 integer,intent(in) :: atindx(dtset%natom) 2433 type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom) 2434 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 2435 2436 !Local variables ------------------------- 2437 !scalars 2438 integer :: iat,iatom,idij,itypat,klmn 2439 !arrays 2440 !-------------------------------------------------------------------- 2441 2442 dterm%aij = czero 2443 2444 do iat=1,dtset%natom 2445 iatom=atindx(iat) 2446 itypat=dtset%typat(iat) 2447 do klmn=1,pawtab(itypat)%lmn2_size 2448 do idij = 1, dterm%ndij 2449 if (paw_ij(iatom)%cplex_dij .EQ. 2) then 2450 dterm%aij(iatom,klmn,idij) = & 2451 & CMPLX(paw_ij(iatom)%dij(2*klmn-1,idij),paw_ij(iatom)%dij(2*klmn,idij)) 2452 else 2453 dterm%aij(iatom,klmn,idij) = CMPLX(paw_ij(iatom)%dij(klmn,idij),zero) 2454 end if 2455 end do ! idij 2456 end do ! klmn 2457 end do ! iat 2458 2459 dterm%has_aij = 2 2460 2461 end subroutine dterm_aij
ABINIT/dterm_alloc [ Functions ]
NAME
dterm_alloc
FUNCTION
allocate space in dterm_type
COPYRIGHT
Copyright (C) 2003-2024 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
lmnmax=max value of lmn over all psps lmn2max=max value of lmn2 over all psps natom=number of atoms in cell ndij=spin channels in dij
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
PARENTS
m_orbmag
CHILDREN
SOURCE
2330 subroutine dterm_alloc(dterm,lmnmax,lmn2max,natom,ndij,pawspnorb) 2331 2332 !Arguments ------------------------------------ 2333 !scalars 2334 integer,intent(in) :: lmnmax,lmn2max,natom,ndij,pawspnorb 2335 type(dterm_type),intent(inout) :: dterm 2336 2337 !arrays 2338 2339 !Local variables ------------------------- 2340 !scalars 2341 2342 !arrays 2343 !-------------------------------------------------------------------- 2344 2345 dterm%lmnmax = lmnmax 2346 dterm%lmn2max = lmn2max 2347 dterm%natom = natom 2348 dterm%ndij = ndij 2349 2350 if(allocated(dterm%aij)) then 2351 ABI_FREE(dterm%aij) 2352 end if 2353 ABI_MALLOC(dterm%aij,(natom,lmn2max,ndij)) 2354 dterm%has_aij=1 2355 2356 if(allocated(dterm%qij)) then 2357 ABI_FREE(dterm%qij) 2358 end if 2359 ABI_MALLOC(dterm%qij,(natom,lmn2max,ndij)) 2360 dterm%has_qij=1 2361 2362 if(allocated(dterm%LR)) then 2363 ABI_FREE(dterm%LR) 2364 end if 2365 ABI_MALLOC(dterm%LR,(natom,lmn2max,ndij,3)) 2366 dterm%has_LR=1 2367 2368 if(allocated(dterm%BM)) then 2369 ABI_FREE(dterm%BM) 2370 end if 2371 ABI_MALLOC(dterm%BM,(natom,lmn2max,ndij,3)) 2372 dterm%has_BM=1 2373 2374 if(allocated(dterm%SOB1)) then 2375 ABI_FREE(dterm%SOB1) 2376 end if 2377 if (pawspnorb /= 0) then 2378 ABI_MALLOC(dterm%SOB1,(natom,lmn2max,ndij,3)) 2379 dterm%has_SOB1=1 2380 else 2381 dterm%has_SOB1=0 2382 end if 2383 2384 2385 end subroutine dterm_alloc
ABINIT/dterm_BM [ Functions ]
NAME
dterm_BM
FUNCTION
Compute onsite <A0.AN>
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset gntselect((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2)=nonzero gaunt integral indices gprimd(3,3)=reciprocal space lattice vectors my_lmax=augmented l_max over all psp pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data realgnt((2*my_lmax-1)**2*(my_lmax)**4)=nonzero gaunt integral values
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
ZTG23 Eq. 43
this term is A0.An = \frac{1}{2}(B x r).\alpha^2(m x r) which can be rewritten
as \frac{\alpha^2}{2} [(B.m)r^2 - B.rr.m]
the first term is bm1 below; the second term involves writing rr as a rank 2 cartesian
tensor, which leads to the rank 0 and rank 2 spherical components below. (that is, write
it as xx,xy,xz,... and write these terms in terms of their real spherical harmonic components)
SOURCE
1801 subroutine dterm_BM(atindx,dterm,dtset,gntselect,gprimd,my_lmax,pawrad,pawtab,realgnt) 1802 1803 !Arguments ------------------------------------ 1804 !scalars 1805 integer,intent(in) :: my_lmax 1806 type(dterm_type),intent(inout) :: dterm 1807 type(dataset_type),intent(in) :: dtset 1808 1809 !arrays 1810 integer,intent(in) :: atindx(dtset%natom) 1811 integer,intent(in) :: gntselect((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2) 1812 real(dp),intent(in) :: gprimd(3,3),realgnt((2*my_lmax-1)**2*(my_lmax)**4) 1813 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 1814 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1815 1816 !Local variables ------------------------- 1817 !scalars 1818 integer :: adir,gint,iat,iatom,ilm,ilmn,itypat,jlmn,jlm 1819 integer :: klmn,klm,kln,lpmp,mesh_size,pwave_size 1820 real(dp) :: a2,bm1,bm2,cfac,d00,d20,d22,dij,intg 1821 1822 !arrays 1823 complex(dpc) :: dij_cart(3),dij_red(3) 1824 real(dp),allocatable :: ff(:) 1825 1826 !-------------------------------------------------------------------- 1827 1828 dterm%BM = czero 1829 a2 = FineStructureConstant2 1830 d00 = sqrt(4.0*pi)/3.0 1831 dij = sqrt(4.0*pi/15.0) 1832 d20 = sqrt(16.0*pi/5.0)/6.0 1833 d22 = sqrt(16.0*pi/15.0)/2.0 1834 1835 do iat = 1, dtset%natom 1836 iatom = atindx(iat) 1837 itypat = dtset%typat(iat) 1838 1839 mesh_size=pawtab(itypat)%mesh_size 1840 pwave_size=size(pawtab(itypat)%phiphj(:,1)) 1841 ABI_MALLOC(ff,(mesh_size)) 1842 do klmn=1, pawtab(itypat)%lmn2_size 1843 klm = pawtab(itypat)%indklmn(1,klmn) 1844 kln = pawtab(itypat)%indklmn(2,klmn) 1845 ilmn = pawtab(itypat)%indklmn(7,klmn) 1846 ilm = pawtab(itypat)%indlmn(4,ilmn) 1847 jlmn = pawtab(itypat)%indklmn(8,klmn) 1848 jlm = pawtab(itypat)%indlmn(4,jlmn) 1849 1850 ff=0 1851 ff(2:pwave_size) = & 1852 & (pawtab(itypat)%phiphj(2:pwave_size,kln)-pawtab(itypat)%tphitphj(2:pwave_size,kln))/& 1853 & (pawrad(itypat)%rad(2:pwave_size)) 1854 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 1855 call simp_gen(intg,ff,pawrad(itypat)) 1856 1857 do adir = 1, 3 1858 if (ilm .EQ. jlm) then 1859 bm1 = half*a2*intg*dtset%nucdipmom(adir,iat) 1860 else 1861 bm1 = zero 1862 end if 1863 1864 bm2 = zero 1865 ! xx, yy, zz cases all have the same contribution from S00 1866 lpmp=1 1867 gint = gntselect(lpmp,klm) 1868 if (gint > 0) then 1869 bm2=bm2+half*a2*dtset%nucdipmom(adir,iat)*d00*realgnt(gint)*intg 1870 end if 1871 ! all other contributions involve Gaunt integrals of S_{2m} 1872 do lpmp = 5, 9 1873 gint = gntselect(lpmp,klm) 1874 if (gint > 0) then 1875 cfac = half*a2*realgnt(gint)*intg 1876 select case (lpmp) 1877 case (5) ! S_{2,-2} contributes to xy term 1878 select case (adir) 1879 case (1) 1880 bm2=bm2+cfac*dtset%nucdipmom(2,iat)*dij 1881 case (2) 1882 bm2=bm2+cfac*dtset%nucdipmom(1,iat)*dij 1883 end select 1884 case (6) ! S_{2,-1} contributes to yz term 1885 select case (adir) 1886 case (2) 1887 bm2=bm2+cfac*dtset%nucdipmom(3,iat)*dij 1888 case (3) 1889 bm2=bm2+cfac*dtset%nucdipmom(2,iat)*dij 1890 end select 1891 case (7) ! S_{2,0} contributes to xx, yy, and zz terms 1892 select case (adir) 1893 case (1) 1894 bm2=bm2-cfac*dtset%nucdipmom(1,iat)*d20 1895 case (2) 1896 bm2=bm2-cfac*dtset%nucdipmom(2,iat)*d20 1897 case (3) 1898 bm2=bm2+cfac*dtset%nucdipmom(3,iat)*2.0*d20 1899 end select 1900 case (8) ! S_{2,+1} contributes to xz term 1901 select case (adir) 1902 case (1) 1903 bm2=bm2+cfac*dtset%nucdipmom(3,iat)*dij 1904 case (3) 1905 bm2=bm2+cfac*dtset%nucdipmom(1,iat)*dij 1906 end select 1907 case (9) ! S_{2,2} contributes to xx, yy terms 1908 select case (adir) 1909 case (1) 1910 bm2=bm2+cfac*dtset%nucdipmom(1,iat)*d22 1911 case (2) 1912 bm2=bm2-cfac*dtset%nucdipmom(2,iat)*d22 1913 end select 1914 end select 1915 end if ! end check on nonzero gaunt integral 1916 end do ! end loop over lp,mp 1917 1918 dij_cart(adir) = -CMPLX(bm1 - bm2,zero) 1919 end do 1920 1921 dij_red = MATMUL(TRANSPOSE(gprimd),dij_cart) 1922 1923 dterm%BM(iatom,klmn,1,1:3) = dij_red(1:3) 1924 if (dterm%ndij > 1) then 1925 dterm%BM(iatom,klmn,2,1:3) = dij_red(1:3) 1926 end if 1927 1928 end do ! end loop over klmn 1929 ABI_FREE(ff) 1930 end do ! end loop over iatom 1931 1932 dterm%has_BM = 2 1933 1934 end subroutine dterm_BM
ABINIT/dterm_free [ Functions ]
NAME
dterm_free
FUNCTION
free space in dterm_type
COPYRIGHT
Copyright (C) 2003-2024 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
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
PARENTS
m_orbmag
CHILDREN
SOURCE
2252 subroutine dterm_free(dterm) 2253 2254 !Arguments ------------------------------------ 2255 !scalars 2256 type(dterm_type),intent(inout) :: dterm 2257 2258 !arrays 2259 2260 !Local variables ------------------------- 2261 !scalars 2262 2263 !arrays 2264 !-------------------------------------------------------------------- 2265 2266 if(allocated(dterm%aij)) then 2267 ABI_FREE(dterm%aij) 2268 end if 2269 dterm%has_aij=0 2270 2271 if(allocated(dterm%qij)) then 2272 ABI_FREE(dterm%qij) 2273 end if 2274 dterm%has_qij=0 2275 2276 if(allocated(dterm%LR)) then 2277 ABI_FREE(dterm%LR) 2278 end if 2279 dterm%has_LR=0 2280 2281 if(allocated(dterm%BM)) then 2282 ABI_FREE(dterm%BM) 2283 end if 2284 dterm%has_BM=0 2285 2286 if(allocated(dterm%SOB1)) then 2287 ABI_FREE(dterm%SOB1) 2288 end if 2289 dterm%has_SOB1=0 2290 2291 end subroutine dterm_free
ABINIT/dterm_LR [ Functions ]
NAME
dterm_LR
FUNCTION
Compute onsite <L_R/2>
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset gprimd(3,3)=reciprocal space lattice vectors pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
ZTG23 text after Eq 42, the on-site angular momentum
SOURCE
1971 subroutine dterm_LR(atindx,dterm,dtset,gprimd,pawrad,pawtab) 1972 1973 !Arguments ------------------------------------ 1974 !scalars 1975 type(dterm_type),intent(inout) :: dterm 1976 type(dataset_type),intent(in) :: dtset 1977 1978 !arrays 1979 integer,intent(in) :: atindx(dtset%natom) 1980 real(dp),intent(in) :: gprimd(3,3) 1981 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 1982 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1983 1984 !Local variables ------------------------- 1985 !scalars 1986 integer :: adir,iat,iatom,ilmn,il,im,itypat,jlmn,jl,jm 1987 integer :: klmn,kln,mesh_size,pwave_size 1988 real(dp) :: intg 1989 complex(dpc) :: orbl_me 1990 1991 !arrays 1992 complex(dpc) :: dij_cart(3),dij_red(3) 1993 real(dp),allocatable :: ff(:) 1994 1995 !-------------------------------------------------------------------- 1996 1997 dterm%LR = czero 1998 1999 do itypat=1,dtset%ntypat 2000 mesh_size=pawtab(itypat)%mesh_size 2001 pwave_size=size(pawtab(itypat)%phiphj(:,1)) 2002 ABI_MALLOC(ff,(mesh_size)) 2003 do klmn=1, pawtab(itypat)%lmn2_size 2004 2005 ilmn = pawtab(itypat)%indklmn(7,klmn) 2006 il=pawtab(itypat)%indlmn(1,ilmn) 2007 im=pawtab(itypat)%indlmn(2,ilmn) 2008 2009 jlmn = pawtab(itypat)%indklmn(8,klmn) 2010 jl=pawtab(itypat)%indlmn(1,jlmn) 2011 jm=pawtab(itypat)%indlmn(2,jlmn) 2012 2013 if ( il /= jl ) cycle ! <l'm'|L|lm> = 0 if l' /= l 2014 if ( il == 0 ) cycle ! <00|L|00> = 0 2015 2016 kln = pawtab(itypat)%indklmn(2,klmn) 2017 ff=0 2018 ff(2:pwave_size) = pawtab(itypat)%phiphj(2:pwave_size,kln)-pawtab(itypat)%tphitphj(2:pwave_size,kln) 2019 call pawrad_deducer0(ff,mesh_size,pawrad(itypat)) 2020 call simp_gen(intg,ff,pawrad(itypat)) 2021 2022 do adir = 1, 3 2023 ! compute <L_dir>/2 2024 call slxyzs(il,im,adir,jl,jm,orbl_me) 2025 dij_cart(adir) = -half*orbl_me*intg 2026 end do ! end loop over adir 2027 2028 ! convert to crystal frame 2029 dij_red = MATMUL(TRANSPOSE(gprimd),dij_cart) 2030 2031 do iat=1,dtset%natom 2032 iatom = atindx(iat) 2033 if(dtset%typat(iat) .EQ. itypat) then 2034 dterm%LR(iatom,klmn,1,1:3) = dij_red(1:3) 2035 if (dterm%ndij > 1) then 2036 dterm%LR(iatom,klmn,2,1:3) = dij_red(1:3) 2037 end if 2038 end if 2039 end do 2040 end do ! end loop over klmn 2041 ABI_FREE(ff) 2042 end do ! end loop over itypat 2043 2044 dterm%has_LR = 2 2045 2046 end subroutine dterm_LR
ABINIT/dterm_qij [ Functions ]
NAME
dterm_qij
FUNCTION
Transfer pawtab%sij to dterm, as complex, solely for convenience
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
Transfer pawtab%sij to dterm, as complex, solely for convenience
SOURCE
1723 subroutine dterm_qij(atindx,dterm,dtset,pawtab) 1724 1725 !Arguments ------------------------------------ 1726 !scalars 1727 type(dterm_type),intent(inout) :: dterm 1728 type(dataset_type),intent(in) :: dtset 1729 1730 !arrays 1731 integer,intent(in) :: atindx(dtset%natom) 1732 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1733 1734 !Local variables ------------------------- 1735 !scalars 1736 integer :: iat,iatom,itypat,lmn2_size 1737 1738 !arrays 1739 1740 !-------------------------------------------------------------------- 1741 1742 dterm%qij = czero 1743 do iat = 1, dtset%natom 1744 iatom = atindx(iat) 1745 itypat = dtset%typat(iat) 1746 lmn2_size = pawtab(itypat)%lmn2_size 1747 dterm%qij(iatom,1:lmn2_size,1) = & 1748 & CMPLX(pawtab(itypat)%sij(1:lmn2_size),zero) 1749 if (dterm%ndij > 1) then 1750 dterm%qij(iatom,1:lmn2_size,2) = & 1751 & CMPLX(pawtab(itypat)%sij(1:lmn2_size),zero) 1752 end if 1753 end do ! iat 1754 1755 dterm%has_qij=2 1756 1757 end subroutine dterm_qij
ABINIT/dterm_SOB1 [ Functions ]
NAME
dterm_SOB1
FUNCTION
compute DIJ SO B1 if called for
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
PARENTS
m_orbmag
CHILDREN
SOURCE
2501 subroutine dterm_SOB1(atindx,cplex_dij,dterm,dtset,paw_an,pawang,pawrad,pawtab,qphase) 2502 2503 !Arguments ------------------------------------ 2504 !scalars 2505 integer,intent(in) :: cplex_dij,qphase 2506 type(dterm_type),intent(inout) :: dterm 2507 type(dataset_type),intent(in) :: dtset 2508 2509 !arrays 2510 integer,intent(in) :: atindx(dtset%natom) 2511 type(paw_an_type),intent(inout) :: paw_an(dtset%natom) 2512 type(pawang_type),intent(in) :: pawang 2513 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 2514 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 2515 2516 !Local variables ------------------------- 2517 !scalars 2518 integer :: angl_size,gint,iat,iatom,idij,idir,igf,ij_size,ipts,itypat 2519 integer :: klmn,klmn1,klm,kln,lmn2_size,mesh_size,ndij 2520 real(dp), parameter :: c1=sqrt(four_pi/five) 2521 real(dp), parameter :: c2=one/sqrt(three) 2522 real(dp), parameter :: c3=sqrt(four_pi) 2523 real(dp), parameter :: HalfFineStruct2=half*FineStructureConstant2 2524 !arrays 2525 character(len=500) :: msg 2526 real(dp) :: rgnt(9),svgt(6) 2527 real(dp),allocatable :: dijsob1(:,:,:),dijsob1_rad(:),dv1dr(:),ff(:) 2528 !-------------------------------------------------------------------- 2529 2530 dterm%SOB1 = czero 2531 !Check data consistency 2532 if (qphase/=1) then 2533 msg='qphase=2 not yet available in Dij SO B1 ' 2534 ABI_BUG(msg) 2535 end if 2536 if (cplex_dij/=2) then 2537 msg='cplex_dij must be 2 for Dij SO B1 ' 2538 ABI_BUG(msg) 2539 end if 2540 if (dterm%ndij/=4) then 2541 msg='ndij must be 4 for Dij SO B1 ' 2542 ABI_BUG(msg) 2543 end if 2544 ! if (size(vxc1,1)/=qphase*mesh_size.or.size(vxc1,3)/=nspden.or.& 2545 ! & (size(vxc1,2)/=angl_size.and.pawxcdev==0).or.& 2546 ! & (size(vxc1,2)/=lm_size.and.pawxcdev/=0)) then 2547 ! msg='invalid sizes for vxc1!' 2548 ! ABI_BUG(msg) 2549 ! end if 2550 2551 do iat=1,dtset%natom 2552 iatom=atindx(iat) 2553 itypat=dtset%typat(iat) 2554 lmn2_size=pawtab(itypat)%lmn2_size 2555 ndij=dterm%ndij 2556 ij_size=pawtab(itypat)%ij_size 2557 mesh_size=pawtab(itypat)%mesh_size 2558 angl_size=pawang%angl_size 2559 2560 if (size(paw_an(iatom)%vh1,1)/=qphase*mesh_size .or. & 2561 & size(paw_an(iatom)%vh1,2)<1 .or. & 2562 & size(paw_an(iatom)%vh1,3)<1) then 2563 msg='invalid sizes for vh1!' 2564 ABI_BUG(msg) 2565 end if 2566 2567 ! compute <Phi_i|r.dV/dr|Phi_j>*alpha2/4*Y_00 (for spin-orbit B1) 2568 ABI_MALLOC(dv1dr,(mesh_size)) 2569 ABI_MALLOC(dijsob1_rad,(ij_size)) 2570 ABI_MALLOC(ff,(mesh_size)) 2571 !fact=one/sqrt(four_pi) ! Y_00 = 1/c3 2572 ff(1:mesh_size)=zero 2573 if (dtset%nspden==1) then 2574 do ipts=1,angl_size 2575 ff(1:mesh_size)=ff(1:mesh_size)+& 2576 & paw_an(iatom)%vxc1(1:mesh_size,ipts,1)*pawang%angwgth(ipts) 2577 end do 2578 else 2579 do ipts=1,angl_size 2580 ff(1:mesh_size)=ff(1:mesh_size)+half*& 2581 & (paw_an(iatom)%vxc1(1:mesh_size,ipts,1)+paw_an(iatom)%vxc1(1:mesh_size,ipts,2))*& 2582 & pawang%angwgth(ipts) 2583 end do 2584 end if 2585 ! ff(1:mesh_size)=c3*ff(1:mesh_size) 2586 ff(1:mesh_size)=(c3*ff(1:mesh_size)+paw_an(iatom)%vh1(1:mesh_size,1,1))/c3 2587 call nderiv_gen(dv1dr,ff,pawrad(itypat)) 2588 dv1dr(2:mesh_size)=half*HalfFineStruct2*& 2589 & (one/(one-ff(2:mesh_size)*half/InvFineStruct**2)**2) & 2590 & *dv1dr(2:mesh_size)*pawrad(itypat)%rad(2:mesh_size) 2591 call pawrad_deducer0(dv1dr,mesh_size,pawrad(itypat)) 2592 do kln=1,ij_size 2593 ff(1:mesh_size)= dv1dr(1:mesh_size)*pawtab(itypat)%phiphj(1:mesh_size,kln) 2594 call simp_gen(dijsob1_rad(kln),ff,pawrad(itypat)) 2595 end do 2596 ABI_FREE(dv1dr) 2597 ABI_FREE(ff) 2598 dijsob1_rad(:)=dtset%spnorbscl*dijsob1_rad(:) 2599 2600 ABI_MALLOC(dijsob1,(3,cplex_dij*qphase*lmn2_size,ndij)) 2601 dijsob1=zero 2602 klmn1=1 2603 do klmn=1,lmn2_size 2604 klm=pawtab(itypat)%indklmn(1,klmn);kln=pawtab(itypat)%indklmn(2,klmn) 2605 ! s00=zero;s2m2=zero;s2m1=zero;s20=zero;s21=zero;s22=zero 2606 rgnt=zero 2607 do igf=1,9 2608 gint=pawang%gntselect(igf,klm) 2609 if (gint .NE. 0) rgnt(igf)=pawang%realgnt(gint) 2610 end do 2611 !rgnt(1)=s00 factor, rgnt(5)=s2m2 etc 2612 ! svgt(1)=xx,svgt(2)=yy,svgt(3)=zz,svgt(4)=yz,svgt(5)=xz,svgt(6)=xy 2613 ! (1-x2/r2) 2614 svgt(1)=c3*rgnt(1)-(c1*c2*rgnt(9)+(c3*rgnt(1)-c1*rgnt(7))/three) 2615 ! (1-y2/r2) 2616 svgt(2)=c3*rgnt(1)-((c3*rgnt(1)-c1*rgnt(7))/three - c1*c2*rgnt(9)) 2617 ! (1-z2/r2) 2618 svgt(3)=c3*rgnt(1)-(two*c1*rgnt(7)+c3*rgnt(1))/three 2619 ! zy/r2 2620 svgt(4)=c1*c2*rgnt(6) 2621 ! zx/r2 2622 svgt(5)=c1*c2*rgnt(8) 2623 ! yx/r2 2624 svgt(6)=c1*c2*rgnt(5) 2625 2626 svgt(:) = svgt(:)*half*dijsob1_rad(kln) 2627 2628 ! bx: sxx - syx - szx 2629 dijsob1(1,klmn1,3)=dijsob1(1,klmn1,3)+svgt(1) 2630 dijsob1(1,klmn1,4)=dijsob1(1,klmn1,4)+svgt(1) 2631 dijsob1(1,klmn1+1,3)=dijsob1(1,klmn1+1,3)+svgt(6) 2632 dijsob1(1,klmn1+1,4)=dijsob1(1,klmn1+1,4)-svgt(6) 2633 dijsob1(1,klmn1,1)=dijsob1(1,klmn1,1)-svgt(5) 2634 dijsob1(1,klmn1,2)=dijsob1(1,klmn1,2)+svgt(5) 2635 ! by: syy - syx - szy 2636 dijsob1(2,klmn1,3)=dijsob1(2,klmn1,3)-svgt(6) 2637 dijsob1(2,klmn1,4)=dijsob1(2,klmn1,4)-svgt(6) 2638 dijsob1(2,klmn1+1,3)=dijsob1(2,klmn1+1,3)-svgt(2) 2639 dijsob1(2,klmn1+1,4)=dijsob1(2,klmn1+1,4)+svgt(2) 2640 dijsob1(2,klmn1,1)=dijsob1(2,klmn1,1)-svgt(4) 2641 dijsob1(2,klmn1,2)=dijsob1(2,klmn1,2)+svgt(4) 2642 ! bz: szz - szy - szx 2643 dijsob1(3,klmn1,3)=dijsob1(3,klmn1,3)-svgt(5) 2644 dijsob1(3,klmn1,4)=dijsob1(3,klmn1,4)-svgt(5) 2645 dijsob1(3,klmn1+1,3)=dijsob1(3,klmn1+1,3)+svgt(4) 2646 dijsob1(3,klmn1+1,4)=dijsob1(3,klmn1+1,4)-svgt(4) 2647 dijsob1(3,klmn1,1)=dijsob1(3,klmn1,1)+svgt(3) 2648 dijsob1(3,klmn1,2)=dijsob1(3,klmn1,2)-svgt(3) 2649 klmn1=klmn1+cplex_dij 2650 2651 do idij = 1, ndij 2652 do idir = 1, 3 2653 dterm%SOB1(iatom,klmn,idij,idir) = & 2654 & CMPLX(dijsob1(idir,2*klmn-1,idij),& 2655 & dijsob1(idir,2*klmn,idij)) 2656 end do ! idir 2657 end do ! idij 2658 2659 end do ! klmn 2660 ABI_FREE(dijsob1_rad) 2661 ABI_FREE(dijsob1) 2662 end do ! iat 2663 2664 dterm%has_SOB1 = 2 2665 2666 end subroutine dterm_SOB1
ABINIT/lamb_core [ Functions ]
NAME
lamb_core
FUNCTION
add core electron contribution to the orbital magnetic moment
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(dtset%natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data
OUTPUT
omlamb(2,3)=contribution of Lamb shielding to magnetic moment
SIDE EFFECTS
TODO
NOTES
lamb shielding of core electrons contributes -m.lambsig to orbital magnetic moment
SOURCE
1452 subroutine lamb_core(atindx,dtset,omlamb,pawtab) 1453 1454 !Arguments ------------------------------------ 1455 !scalars 1456 type(dataset_type),intent(in) :: dtset 1457 1458 !arrays 1459 integer,intent(in) :: atindx(dtset%natom) 1460 real(dp),intent(out) :: omlamb(2,3) 1461 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1462 1463 !Local variables ------------------------- 1464 !scalars 1465 integer :: adir,iat,iatom,itypat 1466 real(dp) :: lambsig 1467 1468 !-------------------------------------------------------------------- 1469 1470 omlamb = zero 1471 do iat=1,dtset%natom 1472 iatom = atindx(iat) 1473 itypat = dtset%typat(iat) 1474 ! if user input lambsig specifically in the input file, use it 1475 if (abs(dtset%lambsig(itypat)).GT.tol8) then 1476 lambsig=dtset%lambsig(itypat) 1477 ! else use the value read in to pawtab structure (which might well be zero) 1478 else 1479 lambsig=pawtab(itypat)%lamb_shielding 1480 end if 1481 do adir = 1, 3 1482 omlamb(1,adir) = omlamb(1,adir) - lambsig*dtset%nucdipmom(adir,iat) 1483 end do ! end loop over adir 1484 end do ! end loop over iat 1485 1486 end subroutine lamb_core
ABINIT/local_fermie [ Functions ]
NAME
local_fermie
FUNCTION
estimate Fermi energy as max value of all occupied bands/kpts
COPYRIGHT
Copyright (C) 2003-2024 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
dtset <type(dataset_type)>=all input variables for this dataset eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)=ground state eigenvalues at each band and kpt mpi_enreg<type(MPI_type)>=information about MPI parallelization occ(dtset%mband*dtset%nkpt*dtset%nsppol)=occup number for each band (often 2) at each k point
OUTPUT
fermie=maximum energy (real(dp)) found over all occupied input bands
SIDE EFFECTS
TODO
NOTES
this routine is parallelized over kpts only, not bands
PARENTS
m_orbmag
CHILDREN
SOURCE
2803 subroutine local_fermie(dtset,eigen0,fermie,mpi_enreg,occ) 2804 2805 !Arguments ------------------------------------ 2806 !scalars 2807 real(dp),intent(out) :: fermie 2808 type(dataset_type),intent(in) :: dtset 2809 type(MPI_type), intent(inout) :: mpi_enreg 2810 2811 !arrays 2812 real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol) 2813 real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 2814 2815 !Local variables ------------------------- 2816 !scalars 2817 integer :: bdtot_index,ierr,ikpt,isppol,me 2818 integer :: nband_k,nn,nproc,spaceComm 2819 real(dp) :: fermie_proc 2820 2821 !arrays 2822 real(dp),allocatable :: eig_k(:),occ_k(:) 2823 2824 !-------------------------------------------------------------------- 2825 2826 spaceComm=mpi_enreg%comm_cell 2827 nproc=xmpi_comm_size(spaceComm) 2828 me = mpi_enreg%me_kpt 2829 2830 bdtot_index=0 2831 fermie_proc = -1.0D99 2832 do isppol = 1, dtset%nsppol 2833 do ikpt = 1, dtset%nkpt 2834 2835 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 2836 2837 ! if the current kpt is not on the current processor, cycle 2838 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 2839 bdtot_index=bdtot_index+nband_k 2840 cycle 2841 end if 2842 2843 ABI_MALLOC(occ_k,(nband_k)) 2844 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 2845 2846 ABI_MALLOC(eig_k,(nband_k)) 2847 eig_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 2848 2849 do nn = 1, nband_k 2850 if ( (abs(occ_k(nn)).GT.tol8) .AND. (eig_k(nn).GT.fermie_proc) ) then 2851 fermie_proc = eig_k(nn) 2852 end if 2853 end do ! nn 2854 2855 ABI_FREE(occ_k) 2856 ABI_FREE(eig_k) 2857 2858 bdtot_index=bdtot_index+nband_k 2859 2860 end do ! end loop over kpts 2861 end do ! end loop over isppol 2862 2863 call xmpi_max(fermie_proc,fermie,spaceComm,ierr) 2864 2865 end subroutine local_fermie
ABINIT/make_d [ Functions ]
NAME
make_d
FUNCTION
this is a driver to compute different onsite terms, in convenient (complex) format
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) dtset <type(dataset_type)>=all input variables for this dataset gprimd(3,3)=reciprocal space lattice vectors paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials
OUTPUT
SIDE EFFECTS
dterm <type(dterm_type)> data related to onsite interactions
TODO
NOTES
PARENTS
m_orbmag
CHILDREN
SOURCE
2708 subroutine make_d(atindx,dterm,dtset,gprimd,paw_an,paw_ij,pawang,pawrad,pawtab,psps) 2709 2710 !Arguments ------------------------------------ 2711 !scalars 2712 type(dterm_type),intent(inout) :: dterm 2713 type(dataset_type),intent(in) :: dtset 2714 type(pseudopotential_type), intent(in) :: psps 2715 2716 !arrays 2717 integer,intent(in) :: atindx(dtset%natom) 2718 real(dp),intent(in) :: gprimd(3,3) 2719 type(paw_an_type),intent(inout) :: paw_an(dtset%natom) 2720 type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom) 2721 type(pawang_type),intent(in) :: pawang 2722 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat) 2723 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 2724 2725 !Local variables ------------------------- 2726 !scalars 2727 integer :: my_lmax,ngnt 2728 2729 !arrays 2730 integer,allocatable :: gntselect(:,:) 2731 real(dp),allocatable :: realgnt(:) 2732 !-------------------------------------------------------------------- 2733 2734 ! make Gaunt integrals 2735 my_lmax = psps%mpsang + 1 2736 ABI_MALLOC(realgnt,((2*my_lmax-1)**2*(my_lmax)**4)) 2737 ABI_MALLOC(gntselect,((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2)) 2738 call realgaunt(my_lmax,ngnt,gntselect,realgnt) 2739 2740 ! generate d terms 2741 2742 ! normal PAW sij overlap, in complex form because it's convenient 2743 call dterm_qij(atindx,dterm,dtset,pawtab) 2744 2745 ! onsite angular momentum expectation values 2746 call dterm_LR(atindx,dterm,dtset,gprimd,pawrad,pawtab) 2747 2748 ! onsite <A_0.A_n> interaction between magnetic field and nuclear dipole 2749 call dterm_BM(atindx,dterm,dtset,gntselect,gprimd,my_lmax,pawrad,pawtab,realgnt) 2750 2751 ! transfers paw_ij to dterm%aij because it's convenient 2752 call dterm_aij(atindx,dterm,dtset,paw_ij,pawtab) 2753 2754 ! generate dterm%SOB1 if necessary 2755 if (dtset%pawspnorb /= 0) then 2756 call dterm_SOB1(atindx,paw_ij(1)%cplex_dij,dterm,dtset,paw_an,pawang,& 2757 & pawrad,pawtab,paw_ij(1)%qphase) 2758 end if 2759 2760 ABI_FREE(realgnt) 2761 ABI_FREE(gntselect) 2762 2763 end subroutine make_d
ABINIT/make_pcg1 [ Functions ]
NAME
make_pcg1
FUNCTION
compute Pc|cg1> from |cg1> and |cg>
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) cg_k(2,mcgk)=ground state wavefunctions at this k point cg1_k(2,mcgk,3)=DDK wavefunctions at this k point, all 3 directions cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k dimlmn(dtset%natom)=cprj lmn dimensions dtset <type(dataset_type)>=all input variables for this dataset gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k ikpt=current k pt isppol=current spin polarization mcgk=dimension of cg_k mcprjk=dimension of cprj_k mkmem_rbz=kpts in memory mpi_enreg<type(MPI_type)>=information about MPI parallelization nband_k=bands at this kpt npw_k=number of planewaves at this kpt occ_k=band occupations at this kpt
OUTPUT
pcg1_k(2,mcgk,3)=cg1_k projected on conduction space
SIDE EFFECTS
TODO
NOTES
see Audouze et al PRB 78, 035105 (2008) Eq. 40
SOURCE
1329 subroutine make_pcg1(atindx,cg_k,cg1_k,cprj_k,dimlmn,dtset,gs_hamk,& 1330 & ikpt,isppol,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,occ_k,pcg1_k) 1331 1332 !Arguments ------------------------------------ 1333 !scalars 1334 integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k 1335 type(dataset_type),intent(in) :: dtset 1336 type(gs_hamiltonian_type),intent(inout) :: gs_hamk 1337 type(MPI_type), intent(inout) :: mpi_enreg 1338 1339 !arrays 1340 integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom) 1341 real(dp),intent(in) :: cg_k(2,mcgk),cg1_k(2,mcgk,3),occ_k(nband_k) 1342 real(dp),intent(out) :: pcg1_k(2,mcgk,3) 1343 type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk) 1344 1345 !Local variables ------------------------- 1346 !scalars 1347 integer :: adir,choice,cpopt,iband,jband 1348 integer :: ndat,nnlout,npwsp,paw_opt,signs,tim_nonlop 1349 real(dp) :: doti,dotr 1350 !arrays 1351 real(dp) :: lambda(1) 1352 real(dp),allocatable :: cwavef(:,:),enlout(:),svectout(:,:) 1353 real(dp),allocatable :: vcg1(:,:),vectout(:,:) 1354 type(pawcprj_type),allocatable :: cwaveprj(:,:) 1355 1356 !-------------------------------------------------------------------- 1357 1358 choice = 5 ! dS/dk in nonlop 1359 cpopt = 4 ! cprj and derivatives already in memory 1360 paw_opt = 3 1361 signs = 2 1362 tim_nonlop = 0 1363 lambda = zero 1364 nnlout = 0 1365 ndat = 1 1366 1367 npwsp = npw_k*dtset%nspinor 1368 1369 ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor)) 1370 call pawcprj_alloc(cwaveprj,3,dimlmn) 1371 ABI_MALLOC(cwavef,(2,npwsp)) 1372 ABI_MALLOC(vectout,(2,npwsp)) 1373 ABI_MALLOC(svectout,(2,npwsp)) 1374 ABI_MALLOC(vcg1,(2,npwsp)) 1375 1376 pcg1_k = zero 1377 1378 do adir = 1, 3 1379 1380 do iband = 1, nband_k 1381 1382 cwavef(1:2,1:npwsp)=cg_k(1:2,(iband-1)*npwsp+1:iband*npwsp) 1383 1384 call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,iband,0,ikpt,0,isppol,dtset%mband,& 1385 & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0) 1386 1387 ! compute S^1|u_i^0> where S^1 = \partial S/\partial k_adir, the k derivative of S in 1388 ! direction adir 1389 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,adir,lambda,mpi_enreg,ndat,& 1390 & nnlout,paw_opt,signs,svectout,tim_nonlop,cwavef,vectout) 1391 1392 !! form vcg1 = -1/2 \sum |u_j^0><u_j^0|S^1|u_i^0>, the valence band part of cg1 1393 vcg1 = zero 1394 do jband = 1, nband_k 1395 if(abs(occ_k(jband)).LT.tol8) cycle 1396 cwavef(1:2,1:npwsp)=cg_k(1:2,(jband-1)*npwsp+1:jband*npwsp) 1397 dotr = DOT_PRODUCT(cwavef(1,:),svectout(1,:))+DOT_PRODUCT(cwavef(2,:),svectout(2,:)) 1398 doti = DOT_PRODUCT(cwavef(1,:),svectout(2,:))-DOT_PRODUCT(cwavef(2,:),svectout(1,:)) 1399 vcg1(1,:) = vcg1(1,:) - half*( dotr*cwavef(1,:) - doti*cwavef(2,:)) 1400 vcg1(2,:) = vcg1(2,:) - half*( dotr*cwavef(2,:) + doti*cwavef(1,:)) 1401 end do 1402 1403 ! subtract vcg1 from cg1_k to obtain pcg1, the conduction band part of cg1 1404 pcg1_k(1:2,(iband-1)*npwsp+1:iband*npwsp,adir) =cg1_k(1:2,(iband-1)*npwsp+1:iband*npwsp,adir)-& 1405 & vcg1(1:2,1:npwsp) 1406 1407 end do 1408 end do 1409 1410 ABI_FREE(cwavef) 1411 ABI_FREE(vectout) 1412 ABI_FREE(vcg1) 1413 ABI_FREE(svectout) 1414 call pawcprj_free(cwaveprj) 1415 ABI_FREE(cwaveprj) 1416 1417 end subroutine make_pcg1
ABINIT/orbmag [ Functions ]
NAME
orbmag
FUNCTION
This routine computes the orbital magnetization and Berry curvature based on input wavefunctions and DDK wavefuntions.
COPYRIGHT
Copyright (C) 2003-2024 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)=all ground state wavefunctions cg1(2,mcg1,3)=all DDK wavefunctions in all 3 directions cprj(dtset%natom,mcprj)<type(pawcprj_type)>=all ground state cprj dtset <type(dataset_type)>=all input variables for this dataset eigen0(dtset%mband*dtset%nkpt*dtset%nsppol)=all ground state eigenvalues gsqcut=large sphere cut-off kg(3,mpw*mkmem_rbz)=basis sphere of planewaves at k mcg=dimension of cg mcg1=dimension of cg1 mcprj=dimension of cprj mkmem_rbz=kpts in memory mpi_enreg<type(MPI_type)>=information about MPI parallelization mpw=max number of planewaves at k nfftf=(effective) number of FFT grid points (for this proc) for the "fine" grid (see NOTES in respfn.F90) ngfftf(18)=FFT grid size information (from pawfgr%ngfft) npwarr(dtset%nkpt)=npw_k at each kpt occ(dtset%mband*dtset%nkpt*dtset%nsppol)=occup number for each band (often 2) at each k point paw_an(my_natom) <type(paw_an_type)>=paw arrays given on angular mesh for the GS paw_ij(dtset%natom) <type(paw_ij_type)>=paw arrays given on (i,j) channels for the GS pawang <type(pawang_type)>=paw angular mesh and related data pawfgr <type(pawfgr_type)>=fine grid parameters and related data pawrad(dtset%ntypat) <type(pawrad_type)>=paw radial mesh and related data pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials rprimd(3,3)=real space translation vectors vtrial(nfftf,dtset%nspden)=GS potential (Hartree) vxctau(nfftf,nspden,4*usekden)=derivative of e_xc with respect to kinetic energy density, for mGGA xred(3,dtset%natom)=reduced dimensionless atomic coordinates ylm(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm)=all ylm's ylmgr(mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm)=gradients of ylm's
OUTPUT
only printing in call to orbmag_output
SIDE EFFECTS
TODO
NOTES
See Zwanziger, Torrent, and Gonze Phys Rev B 107, 165157 (2023), "ZTG23" DDK wavefunctions are used for the derivatives.
SOURCE
223 subroutine orbmag(cg,cg1,cprj,dtset,eigen0,gsqcut,kg,mcg,mcg1,mcprj,mkmem_rbz,mpi_enreg,mpw,& 224 & nfftf,ngfftf,npwarr,occ,paw_an,paw_ij,pawang,pawfgr,pawrad,& 225 & pawtab,psps,rprimd,vtrial,vxctau,xred,ylm,ylmgr) 226 227 !Arguments ------------------------------------ 228 !scalars 229 integer,intent(in) :: mcprj,mcg,mcg1,mkmem_rbz,mpw,nfftf 230 real(dp),intent(in) :: gsqcut 231 type(dataset_type),intent(in) :: dtset 232 type(MPI_type), intent(inout) :: mpi_enreg 233 type(pawfgr_type),intent(in) :: pawfgr 234 type(pseudopotential_type), intent(in) :: psps 235 236 !arrays 237 integer,intent(in) :: kg(3,mpw*mkmem_rbz),ngfftf(18),npwarr(dtset%nkpt) 238 real(dp),intent(in) :: cg(2,mcg),cg1(2,mcg1,3) 239 real(dp),intent(in) :: eigen0(dtset%mband*dtset%nkpt*dtset%nsppol) 240 real(dp), intent(in) :: occ(dtset%mband*dtset%nkpt*dtset%nsppol) 241 real(dp),intent(in) :: rprimd(3,3),xred(3,dtset%natom) 242 real(dp),intent(inout) :: vtrial(nfftf,dtset%nspden) 243 real(dp),intent(inout) :: vxctau(nfftf,dtset%nspden,4*dtset%usekden) 244 real(dp),intent(in) :: ylm(mpw*mkmem_rbz,psps%mpsang*psps%mpsang*psps%useylm) 245 real(dp),intent(in) :: ylmgr(mpw*mkmem_rbz,3,psps%mpsang*psps%mpsang*psps%useylm) 246 type(pawcprj_type),intent(in) :: cprj(dtset%natom,mcprj) 247 type(paw_an_type),intent(inout) :: paw_an(dtset%natom) 248 type(paw_ij_type),intent(inout) :: paw_ij(dtset%natom*psps%usepaw) 249 type(pawang_type),intent(in) :: pawang 250 type(pawrad_type),intent(in) :: pawrad(dtset%ntypat*psps%usepaw) 251 type(pawtab_type),intent(in) :: pawtab(psps%ntypat*psps%usepaw) 252 253 !Local 254 !scalars 255 integer :: adir,bdtot_index,buff_size,choice,cpopt,dimffnl,exchn2n3d 256 integer :: iat,iatom,icg,icmplx,icprj,ider,idir,ierr 257 integer :: ikg,ikg1,ikpt,ilm,indx,isppol,istwf_k,iterm,itypat,lmn2max 258 integer :: me,mcgk,mcprjk,my_lmax,my_nspinor,nband_k,nband_me,ngfft1,ngfft2,ngfft3,ngfft4 259 integer :: ngfft5,ngfft6,ngnt,nl1_option,nn,nkpg,npw_k,npwsp,nproc,spaceComm 260 real(dp) :: arg,ecut_eff,fermie,ucvol 261 logical :: has_nucdip,has_vxctau 262 type(dterm_type) :: dterm 263 type(gs_hamiltonian_type) :: gs_hamk 264 265 !arrays 266 integer,allocatable :: atindx(:),atindx1(:),dimlmn(:),gntselect(:,:),kg_k(:,:),nattyp(:) 267 real(dp) :: gmet(3,3),gprimd(3,3),kpoint(3),omlamb(2,3),rmet(3,3) 268 real(dp),allocatable :: buffer1(:),buffer2(:) 269 real(dp),allocatable :: b1_k(:,:,:),b2_k(:,:,:) 270 real(dp),allocatable :: cg_k(:,:),cg1_k(:,:,:),cwavef(:,:) 271 real(dp),allocatable :: eig_k(:),ffnl_k(:,:,:,:),kinpw(:),kpg_k(:,:) 272 real(dp),allocatable :: m1_k(:,:,:),m1_mu_k(:,:,:),m2_k(:,:,:),m2_mu_k(:,:,:) 273 real(dp),allocatable :: occ_k(:),orbmag_terms(:,:,:,:,:),orbmag_trace(:,:,:) 274 real(dp),allocatable :: pcg1_k(:,:,:),ph1d(:,:),ph3d(:,:,:),phkxred(:,:),realgnt(:) 275 real(dp),allocatable :: vectornd(:,:,:),vectornd_pac(:,:,:,:,:),vlocal(:,:,:,:) 276 real(dp),allocatable :: vxctaulocal(:,:,:,:,:) 277 real(dp),allocatable :: ylm_k(:,:),ylmgr_k(:,:,:) 278 type(pawcprj_type),allocatable :: cprj_k(:,:),cprj1_k(:,:,:),cwaveprj(:,:) 279 280 !---------------------------------------------- 281 282 ! set up basic FFT parameters 283 my_nspinor=max(1,dtset%nspinor/mpi_enreg%nproc_spinor) 284 istwf_k = 1 285 spaceComm=mpi_enreg%comm_cell 286 nproc=xmpi_comm_size(spaceComm) 287 me = mpi_enreg%me_kpt 288 ngfft1=dtset%ngfft(1) ; ngfft2=dtset%ngfft(2) ; ngfft3=dtset%ngfft(3) 289 ngfft4=dtset%ngfft(4) ; ngfft5=dtset%ngfft(5) ; ngfft6=dtset%ngfft(6) 290 ecut_eff = dtset%ecut*(dtset%dilatmx)**2 291 exchn2n3d = 0; ikg1 = 0 292 293 ! Fermi energy 294 call local_fermie(dtset,eigen0,fermie,mpi_enreg,occ) 295 !fermie = dtset%userra 296 !write(std_out,'(a,es16.8)')'JWZ debug using fermi e = ',fermie 297 298 !Definition of atindx array 299 !Generate an index table of atoms, in order for them to be used type after type. 300 ABI_MALLOC(atindx,(dtset%natom)) 301 ABI_MALLOC(atindx1,(dtset%natom)) 302 ABI_MALLOC(nattyp,(psps%ntypat)) 303 indx=1 304 do itypat=1,psps%ntypat 305 nattyp(itypat)=0 306 do iatom=1,dtset%natom 307 if(dtset%typat(iatom)==itypat)then 308 atindx(iatom)=indx 309 atindx1(indx)=iatom 310 indx=indx+1 311 nattyp(itypat)=nattyp(itypat)+1 312 end if 313 end do 314 end do 315 ABI_MALLOC(ph1d,(2,dtset%natom*(2*(ngfft1+ngfft2+ngfft3)+3))) 316 call getph(atindx,dtset%natom,ngfft1,ngfft2,ngfft3,ph1d,xred) 317 318 ABI_MALLOC(kg_k,(3,mpw)) 319 ABI_MALLOC(kinpw,(mpw)) 320 321 ABI_MALLOC(dimlmn,(dtset%natom)) 322 call pawcprj_getdim(dimlmn,dtset%natom,nattyp,dtset%ntypat,dtset%typat,pawtab,'O') 323 324 ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor)) 325 call pawcprj_alloc(cwaveprj,0,dimlmn) 326 327 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 328 329 my_lmax = psps%mpsang + 1 330 ABI_MALLOC(realgnt,((2*my_lmax-1)**2*(my_lmax)**4)) 331 ABI_MALLOC(gntselect,((2*my_lmax-1)**2,my_lmax**2*(my_lmax**2+1)/2)) 332 call realgaunt(my_lmax,ngnt,gntselect,realgnt) 333 334 lmn2max = psps%lmnmax*(psps%lmnmax+1)/2 335 ! note: in make_d, terms will be filled as iatom using atindx 336 call dterm_alloc(dterm,psps%lmnmax,lmn2max,dtset%natom,paw_ij(1)%ndij,dtset%pawspnorb) 337 call make_d(atindx,dterm,dtset,gprimd,paw_an,paw_ij,pawang,pawrad,pawtab,psps) 338 339 ABI_MALLOC(orbmag_terms,(2,dtset%mband,dtset%nsppol,3,nterms)) 340 orbmag_terms = zero 341 342 !==== Initialize most of the Hamiltonian ==== 343 !Allocate all arrays and initialize quantities that do not depend on k and spin. 344 !gs_hamk is the normal hamiltonian at k 345 call init_hamiltonian(gs_hamk,psps,pawtab,dtset%nspinor,dtset%nsppol,dtset%nspden,dtset%natom,& 346 & dtset%typat,xred,dtset%nfft,dtset%mgfft,dtset%ngfft,rprimd,dtset%nloalg,nucdipmom=dtset%nucdipmom,& 347 & paw_ij=paw_ij) 348 349 ! iterate over spin channels 350 bdtot_index=0 351 icg = 0 352 icprj = 0 353 do isppol = 1, dtset%nsppol 354 355 !========= construct local potential ================== 356 ABI_MALLOC(vlocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc)) 357 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, & 358 & dtset%nspden, gs_hamk%nvloc, 1, pawfgr, mpi_enreg, vtrial, vlocal) 359 call gs_hamk%load_spin(isppol,vlocal=vlocal,with_nonlocal=.true.) 360 361 !======== compute nuclear dipole vector potential (may be zero) ========== 362 has_nucdip = ANY( ABS(dtset%nucdipmom) .GT. tol8 ) 363 if(has_nucdip) then 364 ABI_MALLOC(vectornd,(nfftf,dtset%nspden,3)) 365 vectornd = zero 366 call make_vectornd(1,gsqcut,psps%usepaw,mpi_enreg,dtset%natom,nfftf,ngfftf,& 367 & dtset%nspden,dtset%nucdipmom,rprimd,vectornd,xred) 368 ABI_MALLOC(vectornd_pac,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,3)) 369 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, & 370 & dtset%nspden, gs_hamk%nvloc, 3, pawfgr, mpi_enreg, vectornd,vectornd_pac) 371 ABI_FREE(vectornd) 372 call gs_hamk%load_spin(isppol,vectornd=vectornd_pac) 373 end if 374 375 !======== compute vxctaulocal if vxctau present ===================== 376 377 has_vxctau = ( size(vxctau) > 0 ) 378 if(has_vxctau) then 379 ABI_MALLOC(vxctaulocal,(ngfft4,ngfft5,ngfft6,gs_hamk%nvloc,4)) 380 call gspot_transgrid_and_pack(isppol, psps%usepaw, dtset%paral_kgb, dtset%nfft, dtset%ngfft, nfftf, & 381 & dtset%nspden, gs_hamk%nvloc, 4, pawfgr, mpi_enreg, vxctau, vxctaulocal) 382 call gs_hamk%load_spin(isppol, vxctaulocal=vxctaulocal) 383 end if 384 385 ikg = 0 386 !============= BIG FAT KPT LOOP :) =========================== 387 do ikpt = 1, dtset%nkpt 388 389 nband_k=dtset%nband(ikpt+(isppol-1)*dtset%nkpt) 390 nband_me = proc_distrb_nband(mpi_enreg%proc_distrb,ikpt,nband_k,isppol,me) 391 392 ! if the current kpt is not on the current processor, cycle 393 if(proc_distrb_cycle(mpi_enreg%proc_distrb,ikpt,1,nband_k,isppol,me)) then 394 bdtot_index=bdtot_index+nband_k 395 cycle 396 end if 397 398 kpoint(:)=dtset%kptns(:,ikpt) 399 npw_k = npwarr(ikpt) 400 npwsp = npw_k*dtset%nspinor 401 402 ! retrieve kg_k at this k point 403 kg_k(1:3,1:npw_k) = kg(1:3,ikg+1:ikg+npw_k) 404 405 ! retrieve ylm at this k point 406 ABI_MALLOC(ylm_k,(npw_k,psps%mpsang*psps%mpsang)) 407 ABI_MALLOC(ylmgr_k,(npw_k,3,psps%mpsang*psps%mpsang*psps%useylm)) 408 do ilm=1,psps%mpsang*psps%mpsang 409 ylm_k(1:npw_k,ilm)=ylm(1+ikg:npw_k+ikg,ilm) 410 ylmgr_k(1:npw_k,1:3,ilm)=ylmgr(1+ikg:npw_k+ikg,1:3,ilm) 411 end do 412 413 ! retrieve occupation numbers at this k point 414 ABI_MALLOC(occ_k,(nband_k)) 415 occ_k(:)=occ(1+bdtot_index:nband_k+bdtot_index) 416 417 ! Compute kinetic energy at kpt 418 kinpw(:) = zero 419 call mkkin(dtset%ecut,dtset%ecutsm,dtset%effmass_free,gmet,kg_k,kinpw,kpoint,npw_k,0,0) 420 421 ! Compute k+G at this k point 422 nkpg = 3 423 ABI_MALLOC(kpg_k,(npw_k,nkpg)) 424 call mkkpg(kg_k,kpg_k,kpoint,nkpg,npw_k) 425 426 ! Make 3d phase factors 427 ABI_MALLOC(phkxred,(2,dtset%natom)) 428 do iat = 1, dtset%natom 429 iatom = atindx(iat) 430 arg=two_pi*DOT_PRODUCT(kpoint,xred(:,iat)) 431 phkxred(1,iatom)=DCOS(arg);phkxred(2,iatom)=DSIN(arg) 432 end do 433 ABI_MALLOC(ph3d,(2,npw_k,dtset%natom)) 434 call ph1d3d(1,dtset%natom,kg_k,dtset%natom,dtset%natom,& 435 & npw_k,ngfft1,ngfft2,ngfft3,phkxred,ph1d,ph3d) 436 437 ! Compute nonlocal form factors ffnl at all (k+G): 438 ider=1 ! ffnl and 1st derivatives 439 idir=4 ! ignored when ider = 0; idir=0 means d ffnl/ dk in reduced units referenced 440 ! to reciprocal translations 441 ! idir=4 meand d ffnl / dk in reduced units referenced to real space 442 ! translations. rfddk = 1 wavefunctions are computed using this convention. 443 dimffnl=4 ! 1 + number of derivatives 444 ABI_MALLOC(ffnl_k,(npw_k,dimffnl,psps%lmnmax,dtset%ntypat)) 445 call mkffnl(psps%dimekb,dimffnl,psps%ekb,ffnl_k,psps%ffspl,& 446 & gmet,gprimd,ider,idir,psps%indlmn,kg_k,kpg_k,kpoint,psps%lmnmax,& 447 & psps%lnmax,psps%mpsang,psps%mqgrid_ff,nkpg,& 448 & npw_k,dtset%ntypat,psps%pspso,psps%qgrid_ff,rmet,& 449 & psps%usepaw,psps%useylm,ylm_k,ylmgr_k) 450 ! - Load k-dependent quantities in the Hamiltonian 451 call gs_hamk%load_k(kpt_k=kpoint(:),istwf_k=istwf_k,npw_k=npw_k,& 452 & kinpw_k=kinpw,kg_k=kg_k,kpg_k=kpg_k,ffnl_k=ffnl_k,ph3d_k=ph3d,& 453 & compute_gbound=.TRUE.) 454 455 ! retrieve ground state wavefunctions at this k point and isppol 456 mcgk = npw_k*nband_k*dtset%nspinor 457 ABI_MALLOC(cg_k,(2,mcgk)) 458 cg_k = cg(1:2,icg+1:icg+mcgk) 459 460 ! retrieve first order wavefunctions at this k point and isppol 461 ABI_MALLOC(cg1_k,(2,mcgk,3)) 462 cg1_k = cg1(1:2,icg+1:icg+mcgk,1:3) 463 464 ! retrieve zeroth order eigenvalues at this k point and isppol 465 ABI_MALLOC(eig_k,(nband_k)) 466 eig_k(:)=eigen0(1+bdtot_index:nband_k+bdtot_index) 467 468 ! retrieve cprj_k at this k point and isppol 469 mcprjk = nband_k*dtset%nspinor 470 ABI_MALLOC(cprj_k,(dtset%natom,mcprjk)) 471 call pawcprj_alloc(cprj_k,cprj(1,1)%ncpgr,dimlmn) 472 call pawcprj_get(atindx,cprj_k,cprj,dtset%natom,1,icprj,ikpt,0,isppol,dtset%mband,& 473 & mkmem_rbz,dtset%natom,nband_k,nband_k,dtset%nspinor,dtset%nsppol,0) 474 475 ! compute P_c|cg1> 476 ABI_MALLOC(pcg1_k,(2,mcgk,3)) 477 if (dtset%berryopt /= -2) then 478 ! DDK wavefunctions computed in DFPT, must be projected onto conduction space 479 call make_pcg1(atindx,cg_k,cg1_k,cprj_k,dimlmn,dtset,gs_hamk,ikpt,isppol,& 480 & mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,occ_k,pcg1_k) 481 else 482 ! we are using berryopt PEAD DDK functions which are already projected by construction 483 pcg1_k(1:2,1:mcgk,1:3) = cg1_k(1:2,1:mcgk,1:3) 484 end if 485 486 ! compute <p|Pc cg1> cprjs 487 ABI_MALLOC(cprj1_k,(dtset%natom,mcprjk,3)) 488 do adir = 1, 3 489 call pawcprj_alloc(cprj1_k(:,:,adir),0,dimlmn) 490 end do 491 choice = 1 492 cpopt = 0 493 idir = 0 494 ABI_MALLOC(cwavef,(2,npwsp)) 495 do nn = 1, nband_k 496 do adir = 1, 3 497 cwavef(1:2,1:npwsp) = pcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,adir) 498 call getcprj(choice,cpopt,cwavef,cwaveprj,ffnl_k,idir,psps%indlmn,istwf_k,& 499 & kg_k,kpg_k,kpoint,psps%lmnmax,dtset%mgfft,mpi_enreg,dtset%natom,nattyp,dtset%ngfft,& 500 & dtset%nloalg,npw_k,dtset%nspinor,dtset%ntypat,phkxred,ph1d,ph3d,ucvol,psps%useylm) 501 call pawcprj_put(atindx,cwaveprj,cprj1_k(:,:,adir),dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,& 502 & mkmem_rbz,dtset%natom,1,nband_k,dimlmn,dtset%nspinor,dtset%nsppol,0) 503 end do 504 end do 505 ABI_FREE(cwavef) 506 507 !-------------------------------------------------------------------------------- 508 ! Finally ready to compute contributions to orbital magnetism and Berry curvature 509 !-------------------------------------------------------------------------------- 510 511 !! check aij against paw_ij 512 !call check_eig_k(atindx,cg_k,cprj_k,dimlmn,dterm,dtset,eig_k,& 513 ! & gs_hamk,ikpt,isppol,mcgk,mpi_enreg,my_nspinor,nband_k,npw_k,pawtab) 514 515 ABI_MALLOC(b1_k,(2,nband_k,3)) 516 ABI_MALLOC(b2_k,(2,nband_k,3)) 517 ABI_MALLOC(m1_k,(2,nband_k,3)) 518 ABI_MALLOC(m1_mu_k,(2,nband_k,3)) 519 ABI_MALLOC(m2_k,(2,nband_k,3)) 520 ABI_MALLOC(m2_mu_k,(2,nband_k,3)) 521 522 ! ZTG23 Eq. 36 term 2 and Eq. 46 term 1 523 call orbmag_cc_k(atindx,b1_k,cprj1_k,dimlmn,dtset,eig_k,fermie,gs_hamk,ikpt,isppol,& 524 & m1_k,m1_mu_k,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,occ_k,pcg1_k,ucvol) 525 526 ! ZTG Eq. 46 term 1 527 orbmag_terms(:,:,isppol,:,ibcc) = orbmag_terms(:,:,isppol,:,ibcc) + b1_k(:,:,:) 528 529 orbmag_terms(:,:,isppol,:,imcc) = orbmag_terms(:,:,isppol,:,imcc) + m1_k(:,:,:) 530 orbmag_terms(:,:,isppol,:,imcc) = orbmag_terms(:,:,isppol,:,imcc) + m1_mu_k(:,:,:) 531 532 ! ZTG23 Eq. 36 terms 3 and 4 and Eq. 46 term 2 533 call orbmag_vv_k(atindx,b1_k,b2_k,cg_k,cprj_k,dimlmn,dtset,eig_k,fermie,gs_hamk,& 534 & ikpt,isppol,m1_k,m1_mu_k,m2_k,m2_mu_k,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,& 535 & npw_k,occ_k,pcg1_k,ucvol) 536 537 ! ZTG23 Eq. 46 term 2 538 orbmag_terms(:,:,isppol,:,ibvv1) = orbmag_terms(:,:,isppol,:,ibvv1) + b1_k(:,:,:) 539 orbmag_terms(:,:,isppol,:,ibvv2) = orbmag_terms(:,:,isppol,:,ibvv2) + b2_k(:,:,:) 540 541 ! ZTG23 Eq. 36 term 3 542 orbmag_terms(:,:,isppol,:,imvv1) = orbmag_terms(:,:,isppol,:,imvv1) + m1_k(:,:,:) 543 orbmag_terms(:,:,isppol,:,imvv1) = orbmag_terms(:,:,isppol,:,imvv1) + m1_mu_k(:,:,:) 544 545 ! ZTG23 Eq. 36 term 4 546 orbmag_terms(:,:,isppol,:,imvv2) = orbmag_terms(:,:,isppol,:,imvv2) + m2_k(:,:,:) 547 orbmag_terms(:,:,isppol,:,imvv2) = orbmag_terms(:,:,isppol,:,imvv2) + m2_mu_k(:,:,:) 548 549 ! ZTG23 Eq. 36 term 1 550 call orbmag_nl_k(atindx,cprj_k,dimlmn,dterm,dtset,eig_k,ikpt,isppol,& 551 & m1_k,mcprjk,mkmem_rbz,nband_k,occ_k,pawtab,ucvol) 552 orbmag_terms(:,:,isppol,:,imnl) = orbmag_terms(:,:,isppol,:,imnl) + m1_k(:,:,:) 553 554 ! ZTG23 text after Eq. 42 555 nl1_option = 1 ! LR 556 call orbmag_nl1_k(atindx,cprj_k,dimlmn,dterm,dtset,ikpt,isppol,m1_k,mcprjk,mkmem_rbz,& 557 & nband_k,nl1_option,occ_k,pawtab,ucvol) 558 orbmag_terms(:,:,isppol,:,imlr) = orbmag_terms(:,:,isppol,:,imlr) + m1_k 559 560 ! ZTG23 Eq. 43 561 nl1_option = 2 ! BM 562 call orbmag_nl1_k(atindx,cprj_k,dimlmn,dterm,dtset,ikpt,isppol,m1_k,mcprjk,mkmem_rbz,& 563 & nband_k,nl1_option,occ_k,pawtab,ucvol) 564 orbmag_terms(:,:,isppol,:,imbm) = orbmag_terms(:,:,isppol,:,imbm) + m1_k 565 566 ! SOB1 567 if (dterm%has_SOB1 == 2) then 568 nl1_option = 3 ! SOB1 569 call orbmag_nl1_k(atindx,cprj_k,dimlmn,dterm,dtset,ikpt,isppol,m1_k,mcprjk,mkmem_rbz,& 570 & nband_k,nl1_option,occ_k,pawtab,ucvol) 571 orbmag_terms(:,:,isppol,:,imsob1) = orbmag_terms(:,:,isppol,:,imsob1) + m1_k 572 end if 573 574 ABI_FREE(b1_k) 575 ABI_FREE(b2_k) 576 ABI_FREE(m1_k) 577 ABI_FREE(m1_mu_k) 578 ABI_FREE(m2_k) 579 ABI_FREE(m2_mu_k) 580 581 icg = icg + mcgk 582 icprj = icprj + mcprjk 583 ikg = ikg + npw_k 584 bdtot_index=bdtot_index+nband_k 585 586 ABI_FREE(cg_k) 587 ABI_FREE(cg1_k) 588 ABI_FREE(pcg1_k) 589 ABI_FREE(eig_k) 590 ABI_FREE(occ_k) 591 ABI_FREE(ylm_k) 592 ABI_FREE(ylmgr_k) 593 ABI_FREE(kpg_k) 594 ABI_FREE(ffnl_k) 595 ABI_FREE(ph3d) 596 ABI_FREE(phkxred) 597 call pawcprj_free(cprj_k) 598 ABI_FREE(cprj_k) 599 do adir = 1, 3 600 call pawcprj_free(cprj1_k(:,:,adir)) 601 end do 602 ABI_FREE(cprj1_k) 603 604 end do ! end loop over kpts 605 606 ABI_FREE(vlocal) 607 if(allocated(vectornd_pac)) then 608 ABI_FREE(vectornd_pac) 609 end if 610 if(allocated(vxctaulocal)) then 611 ABI_FREE(vxctaulocal) 612 end if 613 614 end do ! end loop over isppol 615 616 !! collect orbmag_terms if distributed over different processes 617 if (nproc > 1) then 618 buff_size=size(orbmag_terms) 619 ABI_MALLOC(buffer1,(buff_size)) 620 ABI_MALLOC(buffer2,(buff_size)) 621 buffer1=zero;buffer2=zero 622 buffer1(1:buff_size) = reshape(orbmag_terms,(/2*nband_k*dtset%nsppol*3*nterms/)) 623 call xmpi_sum(buffer1,buffer2,buff_size,spaceComm,ierr) 624 orbmag_terms(1:2,1:nband_k,1:dtset%nsppol,1:3,1:nterms)=reshape(buffer2,(/2,nband_k,dtset%nsppol,3,nterms/)) 625 ABI_FREE(buffer1) 626 ABI_FREE(buffer2) 627 end if 628 629 !! convert to cartesian frame from reduced triclinic 630 ! general results: [rprimd]*r_red = r_cart 631 ! [gprimd]*k_red = k_cart 632 ! [gprimd]*grad_(r_red) = grad_(r_cart) 633 ! [rprimd]*grad_(k_red) = grad_(k_cart) 634 ! most terms in orbital magnetism look like 635 ! \grad_(k_cart) x \grad_(k_cart) but are calculated in reduced k 636 ! so [rprimd]*grad_(k_red) x [rprimd]*grad_(k_red) = det[rprimd]*[rprimd]^{-1,T} grad_(k_red) x grad_(k_red) 637 ! 638 do iterm = 1, nterms 639 if ((iterm.EQ.iomlmb)) cycle 640 if ((iterm.EQ.imsob1)) cycle 641 do isppol = 1, dtset%nsppol 642 do nn = 1, nband_k 643 do icmplx = 1, 2 644 if((iterm.EQ.imlr).OR.(iterm.EQ.imbm)) then 645 orbmag_terms(icmplx,nn,isppol,1:3,iterm) = & 646 & MATMUL(rprimd,orbmag_terms(icmplx,nn,isppol,1:3,iterm)) 647 else 648 orbmag_terms(icmplx,nn,isppol,1:3,iterm) = & 649 & ucvol*MATMUL(gprimd,orbmag_terms(icmplx,nn,isppol,1:3,iterm)) 650 end if 651 end do !icmplx 652 end do ! nn 653 end do !isppol 654 end do 655 656 !! convert orbmag magnetization to orbital moment 657 !! Berry curvature terms are ignored 658 !! Lamb term ignored 659 do iterm = 1, nterms 660 if (iterm .EQ. ibcc) cycle 661 if (iterm .EQ. ibvv1) cycle 662 if (iterm .EQ. ibvv2) cycle 663 if (iterm .EQ. iomlmb) cycle 664 orbmag_terms(:,:,:,:,iterm) = ucvol*orbmag_terms(:,:,:,:,iterm) 665 end do 666 667 ! compute trace over filled states of each term 668 ABI_MALLOC(orbmag_trace,(2,3,nterms)) 669 orbmag_trace = zero 670 do isppol = 1, dtset%nsppol 671 do nn = 1, nband_k 672 orbmag_trace = orbmag_trace + orbmag_terms(:,nn,isppol,:,:) 673 end do ! nn 674 end do ! isppol 675 676 ! get the Lamb term 677 call lamb_core(atindx,dtset,omlamb,pawtab) 678 orbmag_trace(:,:,iomlmb) = omlamb 679 680 call orbmag_output(dtset,orbmag_terms,orbmag_trace) 681 682 !--------------------------------------------------- 683 ! deallocate memory 684 !--------------------------------------------------- 685 686 call gs_hamk%free() 687 688 ABI_FREE(kg_k) 689 ABI_FREE(kinpw) 690 ABI_FREE(ph1d) 691 692 ABI_FREE(realgnt) 693 ABI_FREE(gntselect) 694 695 ABI_FREE(atindx) 696 ABI_FREE(atindx1) 697 ABI_FREE(nattyp) 698 699 ABI_FREE(dimlmn) 700 call pawcprj_free(cwaveprj) 701 ABI_FREE(cwaveprj) 702 703 call dterm_free(dterm) 704 705 ABI_FREE(orbmag_terms) 706 ABI_FREE(orbmag_trace) 707 708 end subroutine orbmag
ABINIT/orbmag_cc_k [ Functions ]
NAME
orbmag_cc_k
FUNCTION
computes <P_c du/dk|H + E*S|P_c du/dk> term in orbital magnetism
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) cprj1_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for Pc cg1_k dimlmn(dtset%natom)=cprj lmn dimensions dtset <type(dataset_type)>=all input variables for this dataset eig_k(nband_k)=gs eigenvalues at this kpt fermie=offset energy to use gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k ikpt=current k pt isppol=current spin polarization mcgk=dimension of cg_k mcprjk=dimension of cprj_k mkmem_rbz=kpts in memory mpi_enreg<type(MPI_type)>=information about MPI parallelization nband_k=bands at this kpt npw_k=number of planewaves at this kpt occ_k=band occupations at this kpt pcg1_k(2,mcgk,3)=cg1_k projected on conduction space ucvol=unit cell volume
OUTPUT
b1_k(2,nband_k,3)=Berry curvature contribution from <P_cu|Pc_u> terms m1_k(2,nband_k,3)=Orb mag contribution from <P_cu|Pc_u> terms m1_mu_k(2,nband_k,3)=Orb mag contribution from <P_cu|Pc_u> terms with offset
SIDE EFFECTS
TODO
NOTES
ZTG23 Eq. 36 term 2 and Eq. 46 term 1
SOURCE
978 subroutine orbmag_cc_k(atindx,b1_k,cprj1_k,dimlmn,dtset,eig_k,fermie,gs_hamk,ikpt,isppol,& 979 & m1_k,m1_mu_k,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,occ_k,pcg1_k,ucvol) 980 981 !Arguments ------------------------------------ 982 !scalars 983 integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k 984 real(dp),intent(in) :: fermie,ucvol 985 type(dataset_type),intent(in) :: dtset 986 type(gs_hamiltonian_type),intent(inout) :: gs_hamk 987 type(MPI_type), intent(inout) :: mpi_enreg 988 989 !arrays 990 integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom) 991 real(dp),intent(in) :: eig_k(nband_k),occ_k(nband_k),pcg1_k(2,mcgk,3) 992 real(dp),intent(out) :: b1_k(2,nband_k,3),m1_k(2,nband_k,3),m1_mu_k(2,nband_k,3) 993 type(pawcprj_type),intent(in) :: cprj1_k(dtset%natom,mcprjk,3) 994 995 !Local variables ------------------------- 996 !scalars 997 integer :: adir,bdir,cpopt,gdir,ndat,nn,npwsp,sij_opt,tim_getghc,type_calc 998 real(dp) :: doti,dotr,epsabg,lams,trnrm 999 complex(dpc) :: b1,m1,m1_mu,prefac_b,prefac_m 1000 !arrays 1001 real(dp),allocatable :: bra(:,:),ghc(:,:),gsc(:,:),gvnlxc(:,:),ket(:,:) 1002 type(pawcprj_type),allocatable :: cwaveprj1(:,:) 1003 1004 !-------------------------------------------------------------------- 1005 1006 b1_k = zero 1007 m1_k = zero 1008 m1_mu_k = zero 1009 npwsp = npw_k*dtset%nspinor 1010 1011 ABI_MALLOC(bra,(2,npwsp)) 1012 ABI_MALLOC(ket,(2,npwsp)) 1013 ABI_MALLOC(ghc,(2,npwsp)) 1014 ABI_MALLOC(gsc,(2,npwsp)) 1015 ABI_MALLOC(gvnlxc,(2,npwsp)) 1016 ABI_MALLOC(cwaveprj1,(dtset%natom,dtset%nspinor)) 1017 call pawcprj_alloc(cwaveprj1,0,dimlmn) 1018 1019 tim_getghc = 0 1020 lams = zero 1021 ndat = 1 1022 sij_opt = 1 1023 type_calc = 0 1024 1025 do adir = 1, 3 1026 do nn = 1, nband_k 1027 trnrm = occ_k(nn)*dtset%wtk(ikpt)/ucvol 1028 if (abs(trnrm).LT.tol8) cycle 1029 1030 m1 = czero 1031 m1_mu = czero 1032 b1 = czero 1033 1034 do bdir = 1, 3 1035 do gdir = 1, 3 1036 epsabg = eijk(adir,bdir,gdir) 1037 if (ABS(epsabg) .LT. half) cycle 1038 prefac_b = cbc*c2*epsabg 1039 prefac_m = com*c2*epsabg 1040 1041 cpopt = 2 1042 ket(1:2,1:npwsp) = pcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,gdir) 1043 1044 call pawcprj_get(atindx,cwaveprj1,cprj1_k(:,:,gdir),dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,& 1045 & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0) 1046 1047 ! compute H|Pc du> and S|Pc du> 1048 call getghc(cpopt,ket,cwaveprj1,ghc,gsc,gs_hamk,gvnlxc,lams,mpi_enreg,& 1049 & ndat,dtset%prtvol,sij_opt,tim_getghc,type_calc) 1050 1051 bra(1:2,1:npwsp) = pcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,bdir) 1052 1053 dotr = DOT_PRODUCT(bra(1,:),ghc(1,:))+DOT_PRODUCT(bra(2,:),ghc(2,:)) 1054 doti = DOT_PRODUCT(bra(1,:),ghc(2,:))-DOT_PRODUCT(bra(2,:),ghc(1,:)) 1055 m1 = m1 + prefac_m*CMPLX(dotr,doti) 1056 1057 dotr = DOT_PRODUCT(bra(1,:),gsc(1,:))+DOT_PRODUCT(bra(2,:),gsc(2,:)) 1058 doti = DOT_PRODUCT(bra(1,:),gsc(2,:))-DOT_PRODUCT(bra(2,:),gsc(1,:)) 1059 b1 = b1 -two*prefac_b*CMPLX(dotr,doti) 1060 m1 = m1 + prefac_m*CMPLX(dotr,doti)*eig_k(nn) 1061 m1_mu = m1_mu - two*prefac_m*CMPLX(dotr,doti)*fermie 1062 1063 end do !gdir 1064 end do !bdir 1065 1066 b1_k(1,nn,adir) = trnrm*real(b1); b1_k(2,nn,adir) = trnrm*aimag(b1) 1067 m1_k(1,nn,adir) = trnrm*real(m1); m1_k(2,nn,adir) = trnrm*aimag(m1) 1068 m1_mu_k(1,nn,adir) = trnrm*real(m1_mu); m1_mu_k(2,nn,adir) = trnrm*aimag(m1_mu) 1069 1070 end do !nn 1071 end do !adir 1072 1073 ABI_FREE(bra) 1074 ABI_FREE(ket) 1075 ABI_FREE(ghc) 1076 ABI_FREE(gsc) 1077 ABI_FREE(gvnlxc) 1078 call pawcprj_free(cwaveprj1) 1079 ABI_FREE(cwaveprj1) 1080 1081 end subroutine orbmag_cc_k
ABINIT/orbmag_nl1_k [ Functions ]
NAME
orbmag_nl1_k
FUNCTION
make NL(1) term at k
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k dimlmn(dtset%natom)=cprj lmn dimensions dterm <type(dterm_type)> data related to onsite interactions dtset <type(dataset_type)>=all input variables for this dataset ikpt=current k pt isppol=current spin polarization mcprjk=dimension of cprj_k mkmem_rbz=kpts in memory nband_k=bands at this kpt nl1_option=chooses which onsite term to apply occ_k=band occupations at this kpt pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data ucvol=unit cell volume
OUTPUT
if nl1_option = 1, orbmag contribution of <L_R> is returned if nl1_option = 2, orbmag contribution of <A0.An> is returned m1_k(2,nband_k,3)=Orb mag contribution
SIDE EFFECTS
TODO
NOTES
returns \sum_{Rij}<u|p_i>a_ij<p_j|u> for various a_ij inputs
SOURCE
755 subroutine orbmag_nl1_k(atindx,cprj_k,dimlmn,dterm,dtset,ikpt,isppol,m1_k,mcprjk,& 756 & mkmem_rbz,nband_k,nl1_option,occ_k,pawtab,ucvol) 757 758 !Arguments ------------------------------------ 759 !scalars 760 integer,intent(in) :: ikpt,isppol,mcprjk,mkmem_rbz,nband_k,nl1_option 761 real(dp),intent(in) :: ucvol 762 type(dterm_type),intent(in) :: dterm 763 type(dataset_type),intent(in) :: dtset 764 765 !arrays 766 integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom) 767 real(dp),intent(in) :: occ_k(nband_k) 768 real(dp),intent(out) :: m1_k(2,nband_k,3) 769 type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk) 770 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 771 772 !Local variables ------------------------- 773 !scalars 774 integer :: adir,nn 775 real(dp) :: trnrm 776 complex(dpc) :: tt 777 !arrays 778 type(pawcprj_type),allocatable :: cwaveprj(:,:) 779 780 !-------------------------------------------------------------------- 781 782 m1_k = zero 783 784 ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor)) 785 call pawcprj_alloc(cwaveprj,cprj_k(1,1)%ncpgr,dimlmn) 786 787 do adir = 1, 3 788 do nn = 1, nband_k 789 trnrm = occ_k(nn)*dtset%wtk(ikpt)/ucvol 790 if(abs(trnrm).LT.tol8) cycle 791 792 call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,& 793 & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0) 794 795 select case (nl1_option) 796 case(1) 797 call tt_me(dterm%LR(:,:,:,adir),atindx,cwaveprj,dtset,cwaveprj,dterm%lmn2max,dterm%ndij,pawtab,tt) 798 case(2) 799 call tt_me(dterm%BM(:,:,:,adir),atindx,cwaveprj,dtset,cwaveprj,dterm%lmn2max,dterm%ndij,pawtab,tt) 800 case(3) 801 call tt_me(dterm%SOB1(:,:,:,adir),atindx,cwaveprj,dtset,cwaveprj,dterm%lmn2max,dterm%ndij,pawtab,tt) 802 case default 803 tt = czero 804 end select 805 806 m1_k(1,nn,adir) = trnrm*real(tt); m1_k(2,nn,adir) = trnrm*aimag(tt) 807 808 end do !nn 809 end do !adir 810 811 call pawcprj_free(cwaveprj) 812 ABI_FREE(cwaveprj) 813 814 end subroutine orbmag_nl1_k
ABINIT/orbmag_nl_k [ Functions ]
NAME
orbmag_nl_k
FUNCTION
make NL term at k
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k dimlmn(dtset%natom)=cprj lmn dimensions dterm <type(dterm_type)> data related to onsite interactions dtset <type(dataset_type)>=all input variables for this dataset eig_k(nband_k)=gs eigenvalues at this kpt ikpt=current k pt isppol=current spin polarization mcprjk=dimension of cprj_k mkmem_rbz=kpts in memory nband_k=bands at this kpt occ_k=band occupations at this kpt pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data ucvol=unit cell volume
OUTPUT
m1_k(2,nband_k,3)=Orb mag contribution
SIDE EFFECTS
TODO
NOTES
computes -\frac{i}{2}\sum_{Rij}<u|d_b p_i>D^0_{ij} - E^0s^0_{ij}<d_g p_j|u>
This is ZTG23 Eq. 36 term 1
SOURCE
862 subroutine orbmag_nl_k(atindx,cprj_k,dimlmn,dterm,dtset,eig_k,ikpt,isppol,& 863 & m1_k,mcprjk,mkmem_rbz,nband_k,occ_k,pawtab,ucvol) 864 865 !Arguments ------------------------------------ 866 !scalars 867 integer,intent(in) :: ikpt,isppol,mcprjk,mkmem_rbz,nband_k 868 real(dp),intent(in) :: ucvol 869 type(dterm_type),intent(in) :: dterm 870 type(dataset_type),intent(in) :: dtset 871 872 !arrays 873 integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom) 874 real(dp),intent(in) :: eig_k(nband_k),occ_k(nband_k) 875 real(dp),intent(out) :: m1_k(2,nband_k,3) 876 type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk) 877 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 878 879 !Local variables ------------------------- 880 !scalars 881 integer :: adir,bdir,gdir,nn 882 real(dp) :: epsabg,trnrm 883 complex(dpc) :: m1,prefac_m,txt_d,txt_q 884 !arrays 885 type(pawcprj_type),allocatable :: cwaveprj(:,:) 886 887 !-------------------------------------------------------------------- 888 889 m1_k = zero 890 ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor)) 891 call pawcprj_alloc(cwaveprj,cprj_k(1,1)%ncpgr,dimlmn) 892 893 do adir = 1, 3 894 do nn = 1, nband_k 895 trnrm = occ_k(nn)*dtset%wtk(ikpt)/ucvol 896 if(abs(trnrm).LT.tol8) cycle 897 call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,& 898 & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0) 899 900 m1 = czero 901 do bdir = 1, 3 902 do gdir = 1, 3 903 epsabg = eijk(adir,bdir,gdir) 904 if (ABS(epsabg) .LT. half) cycle 905 prefac_m = com*c2*epsabg 906 907 call txt_me(dterm%aij,atindx,cwaveprj,bdir,dtset,gdir,cwaveprj,& 908 & dterm%lmn2max,dterm%ndij,pawtab,txt_d) 909 910 call txt_me(dterm%qij,atindx,cwaveprj,bdir,dtset,gdir,cwaveprj,& 911 & dterm%lmn2max,dterm%ndij,pawtab,txt_q) 912 913 ! note rho^0 H^1 term has opposite sign of rho^1 H^0 914 m1 = m1 - prefac_m*(txt_d - eig_k(nn)*txt_q) 915 916 end do !gdir 917 end do !bdir 918 919 m1_k(1,nn,adir) = trnrm*real(m1); m1_k(2,nn,adir) = trnrm*aimag(m1) 920 921 end do !nn 922 end do !adir 923 924 call pawcprj_free(cwaveprj) 925 ABI_FREE(cwaveprj) 926 927 end subroutine orbmag_nl_k
ABINIT/orbmag_output [ Functions ]
NAME
orbmag_output
FUNCTION
Only printing. This routine outputs orbmag terms to the normal abi out file
COPYRIGHT
Copyright (C) 2003-2024 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
dtset <type(dataset_type)>=all input variables for this dataset orbmag_terms(2,dtset%mband,dtset%nsppol,3,nterms)=all computed terms per band of orb mag orbmag_trace(2,3,nterms)=traces of orbmag_terms
OUTPUT
SIDE EFFECTS
TODO
NOTES
SOURCE
2078 subroutine orbmag_output(dtset,orbmag_terms,orbmag_trace) 2079 2080 2081 !Arguments ------------------------------------ 2082 !scalars 2083 type(dataset_type),intent(in) :: dtset 2084 2085 !arrays 2086 real(dp),intent(in) :: orbmag_terms(2,dtset%mband,dtset%nsppol,3,nterms),orbmag_trace(2,3,nterms) 2087 2088 !Local variables ------------------------- 2089 !scalars 2090 integer :: adir,iband,isppol,iterms 2091 character(len=500) :: message 2092 2093 !arrays 2094 real(dp) :: berry_bb(2,dtset%mband,3),berry_total(2,3),orbmag_bb(2,dtset%mband,3),orbmag_total(2,3) 2095 2096 ! *********************************************************************** 2097 2098 orbmag_bb=zero;orbmag_total=zero 2099 do iterms = ibvv2+1, nterms 2100 orbmag_total(1:2,1:3)=orbmag_total(1:2,1:3) + orbmag_trace(1:2,1:3,iterms) 2101 do isppol = 1, dtset%nsppol 2102 do iband=1, dtset%mband 2103 orbmag_bb(1:2,iband,1:3) = orbmag_bb(1:2,iband,1:3) + orbmag_terms(1:2,iband,isppol,1:3,iterms) 2104 end do ! iband 2105 end do ! isppol 2106 end do 2107 berry_bb=zero;berry_total=zero 2108 do iterms = ibcc,ibvv2 2109 berry_total(1:2,1:3)=berry_total(1:2,1:3) + orbmag_trace(1:2,1:3,iterms) 2110 do isppol = 1, dtset%nsppol 2111 do iband=1, dtset%mband 2112 berry_bb(1:2,iband,1:3) = berry_bb(1:2,iband,1:3) + orbmag_terms(1:2,iband,isppol,1:3,iterms) 2113 end do ! iband 2114 end do ! isppol 2115 end do 2116 2117 write(message,'(a,a,a)')ch10,'====================================================',ch10 2118 call wrtout(ab_out,message,'COLL') 2119 2120 write(message,'(a,a)')' Orbital magnetic moment computed with DFPT derivative wavefunctions ',ch10 2121 call wrtout(ab_out,message,'COLL') 2122 2123 write(message,'(a)')' Orbital magnetic moment, Cartesian directions : ' 2124 call wrtout(ab_out,message,'COLL') 2125 write(message,'(3es16.8)') (orbmag_total(1,adir),adir=1,3) 2126 call wrtout(ab_out,message,'COLL') 2127 write(message,'(a)')ch10 2128 call wrtout(ab_out,message,'COLL') 2129 write(message,'(a)')' Chern vector, Cartesian directions : ' 2130 call wrtout(ab_out,message,'COLL') 2131 write(message,'(3es16.8)') (berry_total(1,adir),adir=1,3) 2132 call wrtout(ab_out,message,'COLL') 2133 2134 if(dtset%orbmag .EQ. 2) then 2135 write(message,'(a)')ch10 2136 call wrtout(ab_out,message,'COLL') 2137 write(message,'(a)')' Orbital magnetic moment, term-by-term breakdown : ' 2138 call wrtout(ab_out,message,'COLL') 2139 write(message,'(a,3es16.8)') ' rho(1) CC : ',(orbmag_trace(1,adir,imcc),adir=1,3) 2140 call wrtout(ab_out,message,'COLL') 2141 write(message,'(a,3es16.8)') ' rho(1) VV1 : ',(orbmag_trace(1,adir,imvv1),adir=1,3) 2142 call wrtout(ab_out,message,'COLL') 2143 write(message,'(a,3es16.8)') ' rho(1) VV2 : ',(orbmag_trace(1,adir,imvv2),adir=1,3) 2144 call wrtout(ab_out,message,'COLL') 2145 write(message,'(a,3es16.8)') ' rho(0) NL : ',(orbmag_trace(1,adir,imnl),adir=1,3) 2146 call wrtout(ab_out,message,'COLL') 2147 write(message,'(a,3es16.8)') ' <L_R> terms : ',(orbmag_trace(1,adir,imlr),adir=1,3) 2148 call wrtout(ab_out,message,'COLL') 2149 write(message,'(a,3es16.8)') ' <A0.An> terms : ',(orbmag_trace(1,adir,imbm),adir=1,3) 2150 call wrtout(ab_out,message,'COLL') 2151 if (dtset%pawspnorb /= 0) then 2152 write(message,'(a,3es16.8)') ' <SO B1> terms : ',(orbmag_trace(1,adir,imsob1),adir=1,3) 2153 call wrtout(ab_out,message,'COLL') 2154 end if 2155 write(message,'(a,3es16.8)') ' Lamb terms : ',(orbmag_trace(1,adir,iomlmb),adir=1,3) 2156 call wrtout(ab_out,message,'COLL') 2157 write(message,'(a)')' Chern vector, term-by-term breakdown : ' 2158 call wrtout(ab_out,message,'COLL') 2159 write(message,'(a,3es16.8)') ' Ch CC : ',(orbmag_trace(1,adir,ibcc),adir=1,3) 2160 call wrtout(ab_out,message,'COLL') 2161 write(message,'(a,3es16.8)') ' Ch vv1 : ',(orbmag_trace(1,adir,ibvv1),adir=1,3) 2162 call wrtout(ab_out,message,'COLL') 2163 write(message,'(a,3es16.8)') ' Ch vv2 : ',(orbmag_trace(1,adir,ibvv2),adir=1,3) 2164 call wrtout(ab_out,message,'COLL') 2165 end if 2166 2167 if(dtset%orbmag .EQ. 2) then 2168 write(message,'(a)')ch10 2169 call wrtout(ab_out,message,'COLL') 2170 write(message,'(a)')' Term-by-term breakdowns for each band : ' 2171 call wrtout(ab_out,message,'COLL') 2172 do iband = 1, dtset%mband 2173 do isppol = 1, dtset%nsppol 2174 write(message,'(a)')ch10 2175 call wrtout(ab_out,message,'COLL') 2176 if (dtset%nsppol .EQ. 2) then 2177 write(message,'(a,i3,a,i3,a,i2,a,i2)') ' band ',iband,' of ',dtset%mband,'; spin polarization ',isppol,' of ',dtset%nsppol 2178 else 2179 write(message,'(a,i3,a,i3)') ' band ',iband,' of ',dtset%mband 2180 end if 2181 call wrtout(ab_out,message,'COLL') 2182 write(message,'(a,3es16.8)') ' Orbital magnetic moment : ',(orbmag_bb(1,iband,adir),adir=1,3) 2183 call wrtout(ab_out,message,'COLL') 2184 write(message,'(a,3es16.8)') ' rho(1) CC : ',(orbmag_terms(1,iband,isppol,adir,imcc),adir=1,3) 2185 call wrtout(ab_out,message,'COLL') 2186 write(message,'(a,3es16.8)') ' rho(1) VV1 : ',(orbmag_terms(1,iband,isppol,adir,imvv1),adir=1,3) 2187 call wrtout(ab_out,message,'COLL') 2188 write(message,'(a,3es16.8)') ' rho(1) VV2 : ',(orbmag_terms(1,iband,isppol,adir,imvv2),adir=1,3) 2189 call wrtout(ab_out,message,'COLL') 2190 write(message,'(a,3es16.8)') ' rho(0) NL : ',(orbmag_terms(1,iband,isppol,adir,imnl),adir=1,3) 2191 call wrtout(ab_out,message,'COLL') 2192 write(message,'(a,3es16.8)') ' <L_R> terms : ',(orbmag_terms(1,iband,isppol,adir,imlr),adir=1,3) 2193 call wrtout(ab_out,message,'COLL') 2194 write(message,'(a,3es16.8)') ' <A0.An> terms : ',(orbmag_terms(1,iband,isppol,adir,imbm),adir=1,3) 2195 call wrtout(ab_out,message,'COLL') 2196 if (dtset%pawspnorb /= 0) then 2197 write(message,'(a,3es16.8)') ' <SO B1> terms : ',(orbmag_terms(1,iband,isppol,adir,imsob1),adir=1,3) 2198 call wrtout(ab_out,message,'COLL') 2199 end if 2200 write(message,'(a)')ch10 2201 call wrtout(ab_out,message,'COLL') 2202 write(message,'(a,3es16.8)') ' Chern vector : ',(berry_bb(1,iband,adir),adir=1,3) 2203 call wrtout(ab_out,message,'COLL') 2204 write(message,'(a,3es16.8)') ' Ch CC : ',(orbmag_terms(1,iband,isppol,adir,ibcc),adir=1,3) 2205 call wrtout(ab_out,message,'COLL') 2206 write(message,'(a,3es16.8)') ' Ch VV1 : ',(orbmag_terms(1,iband,isppol,adir,ibvv1),adir=1,3) 2207 call wrtout(ab_out,message,'COLL') 2208 write(message,'(a,3es16.8)') ' Ch VV2 : ',(orbmag_terms(1,iband,isppol,adir,ibvv2),adir=1,3) 2209 call wrtout(ab_out,message,'COLL') 2210 end do !isppol 2211 end do ! iband 2212 end if 2213 2214 write(message,'(a,a,a)')ch10,'====================================================',ch10 2215 call wrtout(ab_out,message,'COLL') 2216 2217 end subroutine orbmag_output
ABINIT/orbmag_vv_k [ Functions ]
NAME
orbmag_vv_k
FUNCTION
orbmag_vv_k
COPYRIGHT
Copyright (C) 2003-2024 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
atindx(natom)=index table for atoms (see gstate.f) cg_k(2,mcgk)=ground state wavefunctions at this k point cprj_k(dtset%natom,mcprjk)<type(pawcprj_type)>=cprj for cg_k dimlmn(dtset%natom)=cprj lmn dimensions dtset <type(dataset_type)>=all input variables for this dataset eig_k(nband_k)=gs eigenvalues at this kpt fermie=offset energy to use gs_hamk<type(gs_hamiltonian_type)>=ground state Hamiltonian at this k ikpt=current k pt isppol=current spin polarization mcgk=dimension of cg_k mcprjk=dimension of cprj_k mkmem_rbz=kpts in memory mpi_enreg<type(MPI_type)>=information about MPI parallelization nband_k=bands at this kpt npw_k=number of planewaves at this kpt occ_k=band occupations at this kpt pcg1_k(2,mcgk,3)=cg1_k projected on conduction space ucvol=unit cell volume
OUTPUT
b1_k(2,nband_k,3)=Berry curvature contribution from <u|Pc_u> terms b2_k(2,nband_k,3)=Berry curvature contribution from <u|d S|u> terms m1_k(2,nband_k,3)=Orb mag contribution from <u|Pc_u> terms m1_mu_k(2,nband_k,3)=Orb mag contribution from <u|Pc_u> terms with offset m2_k(2,nband_k,3)=Orb mag contribution from <u|dS|u> terms m2_mu_k(2,nband_k,3)=Orb mag contribution from <u|dS|u> terms with offset
SIDE EFFECTS
TODO
NOTES
contributions (1) <Pc d_b u|E d_gS|u> + <u|E d_bS|Pc d_g u> and (2) \sum_n' <u |d_b ES|u_n'><u_n'|d_g ES|u> to orbital magnetization these are ZTG23 Eq 36 terms 3 and 4, and Eq. 46 term 2
SOURCE
1138 subroutine orbmag_vv_k(atindx,b1_k,b2_k,cg_k,cprj_k,dimlmn,dtset,eig_k,fermie,gs_hamk,& 1139 & ikpt,isppol,m1_k,m1_mu_k,m2_k,m2_mu_k,mcgk,mcprjk,mkmem_rbz,mpi_enreg,nband_k,npw_k,& 1140 & occ_k,pcg1_k,ucvol) 1141 1142 !Arguments ------------------------------------ 1143 !scalars 1144 integer,intent(in) :: ikpt,isppol,mcgk,mcprjk,mkmem_rbz,nband_k,npw_k 1145 real(dp),intent(in) :: fermie,ucvol 1146 type(dataset_type),intent(in) :: dtset 1147 type(gs_hamiltonian_type),intent(inout) :: gs_hamk 1148 type(MPI_type), intent(inout) :: mpi_enreg 1149 1150 !arrays 1151 integer,intent(in) :: atindx(dtset%natom),dimlmn(dtset%natom) 1152 real(dp),intent(in) :: cg_k(2,mcgk),eig_k(nband_k),occ_k(nband_k),pcg1_k(2,mcgk,3) 1153 real(dp),intent(out) :: b1_k(2,nband_k,3),b2_k(2,nband_k,3) 1154 real(dp),intent(out) :: m1_k(2,nband_k,3),m1_mu_k(2,nband_k,3) 1155 real(dp),intent(out) :: m2_k(2,nband_k,3),m2_mu_k(2,nband_k,3) 1156 type(pawcprj_type),intent(in) :: cprj_k(dtset%natom,mcprjk) 1157 1158 !Local variables ------------------------- 1159 !scalars 1160 integer :: adir,bdir,choice,cpopt,gdir,ndat,nn,nnlout,np,npwsp 1161 integer :: paw_opt,signs,tim_getghc 1162 real(dp) :: doti,dotr,epsabg,trnrm 1163 complex(dpc) :: b1,bv2b,m1,m1_mu,mb,mg,mv2b,mv2b_mu,prefac_b,prefac_m 1164 !arrays 1165 real(dp) :: enlout(1),lamv(1) 1166 real(dp),allocatable :: bra(:,:),ket(:,:),svectoutb(:,:),svectoutg(:,:),vectout(:,:) 1167 type(pawcprj_type),allocatable :: cwaveprj(:,:) 1168 1169 !-------------------------------------------------------------------- 1170 1171 m1_k = zero; m2_k = zero 1172 m1_mu_k = zero; m2_mu_k = zero 1173 b1_k = zero; b2_k = zero 1174 npwsp = npw_k*dtset%nspinor 1175 1176 ABI_MALLOC(bra,(2,npwsp)) 1177 ABI_MALLOC(ket,(2,npwsp)) 1178 ABI_MALLOC(svectoutb,(2,npwsp)) 1179 ABI_MALLOC(svectoutg,(2,npwsp)) 1180 ABI_MALLOC(vectout,(2,npwsp)) 1181 ABI_MALLOC(cwaveprj,(dtset%natom,dtset%nspinor)) 1182 call pawcprj_alloc(cwaveprj,cprj_k(1,1)%ncpgr,dimlmn) 1183 1184 tim_getghc = 0 1185 lamv = zero 1186 ndat = 1 1187 cpopt = 4 ! cprj and derivs in memory 1188 choice = 5 ! apply dS/dk 1189 paw_opt = 3 ! retain dS/dk|u> 1190 signs = 2 1191 nnlout = 1 1192 1193 do adir = 1, 3 1194 do nn = 1, nband_k 1195 1196 m1 = czero; mv2b = czero 1197 m1_mu = czero; mv2b_mu = czero 1198 b1 = czero; bv2b = czero 1199 trnrm = occ_k(nn)*dtset%wtk(ikpt)/ucvol 1200 if(abs(trnrm).LT.tol8) cycle 1201 1202 do bdir = 1, 3 1203 do gdir = 1, 3 1204 epsabg = eijk(adir,bdir,gdir) 1205 if (ABS(epsabg) .LT. half) cycle 1206 prefac_b = cbc*c2*epsabg 1207 prefac_m = com*c2*epsabg 1208 1209 ! extract |u_nk> 1210 ket(1:2,1:npwsp) = cg_k(1:2,(nn-1)*npwsp+1:nn*npwsp) 1211 call pawcprj_get(atindx,cwaveprj,cprj_k,dtset%natom,nn,0,ikpt,0,isppol,dtset%mband,& 1212 & mkmem_rbz,dtset%natom,1,nband_k,dtset%nspinor,dtset%nsppol,0) 1213 1214 ! compute dS/dk_g |u_nk> 1215 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,gdir,lamv,mpi_enreg,ndat,nnlout,& 1216 & paw_opt,signs,svectoutg,tim_getghc,ket,vectout) 1217 ! compute dS/dk_b|u_nk> 1218 call nonlop(choice,cpopt,cwaveprj,enlout,gs_hamk,bdir,lamv,mpi_enreg,ndat,nnlout,& 1219 & paw_opt,signs,svectoutb,tim_getghc,ket,vectout) 1220 1221 ! extract |Pc du/dk_b> 1222 bra(1:2,1:npwsp) = pcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,bdir) 1223 !overlap 1224 dotr = DOT_PRODUCT(bra(1,:),svectoutg(1,:))+DOT_PRODUCT(bra(2,:),svectoutg(2,:)) 1225 doti = DOT_PRODUCT(bra(1,:),svectoutg(2,:))-DOT_PRODUCT(bra(2,:),svectoutg(1,:)) 1226 ! add <Pc du/dk_b|dS/dk_g|u_nk>*E_nk 1227 b1 = b1 - prefac_b*CMPLX(dotr,doti) 1228 m1 = m1 + prefac_m*CMPLX(dotr,doti)*eig_k(nn) 1229 m1_mu = m1_mu - prefac_m*CMPLX(dotr,doti)*fermie 1230 1231 ! extract |Pc du/dk_g> 1232 bra(1:2,1:npwsp) = pcg1_k(1:2,(nn-1)*npwsp+1:nn*npwsp,gdir) 1233 !overlap 1234 dotr = DOT_PRODUCT(bra(1,:),svectoutb(1,:))+DOT_PRODUCT(bra(2,:),svectoutb(2,:)) 1235 doti = DOT_PRODUCT(bra(1,:),svectoutb(2,:))-DOT_PRODUCT(bra(2,:),svectoutb(1,:)) 1236 1237 ! add CONJG(<Pc du/dk_b|dS/dk_g|u_nk>)*E_nk 1238 b1 = b1 - prefac_b*CMPLX(dotr,-doti) 1239 m1 = m1 + prefac_m*CMPLX(dotr,-doti)*eig_k(nn) 1240 m1_mu = m1_mu - prefac_m*CMPLX(dotr,-doti)*fermie 1241 1242 do np = 1, nband_k 1243 if(abs(occ_k(np)).LT.tol8) cycle 1244 1245 bra(1:2,1:npwsp) = cg_k(1:2,(np-1)*npwsp+1:np*npwsp) 1246 1247 dotr = DOT_PRODUCT(bra(1,:),svectoutg(1,:))+DOT_PRODUCT(bra(2,:),svectoutg(2,:)) 1248 doti = DOT_PRODUCT(bra(1,:),svectoutg(2,:))-DOT_PRODUCT(bra(2,:),svectoutg(1,:)) 1249 mg = CMPLX(dotr,doti) 1250 1251 dotr = DOT_PRODUCT(bra(1,:),svectoutb(1,:))+DOT_PRODUCT(bra(2,:),svectoutb(2,:)) 1252 doti = DOT_PRODUCT(bra(1,:),svectoutb(2,:))-DOT_PRODUCT(bra(2,:),svectoutb(1,:)) 1253 mb = CMPLX(dotr,doti) 1254 1255 ! terms in <u|dS|u'><u'|dS|u> 1256 bv2b = bv2b + prefac_b*CONJG(mb)*mg 1257 mv2b = mv2b - prefac_m*CONJG(mb)*mg*eig_k(nn) 1258 mv2b_mu = mv2b_mu + prefac_m*CONJG(mb)*mg*fermie 1259 end do ! np 1260 1261 end do !gdir 1262 end do !bdir 1263 1264 b1_k(1,nn,adir) = trnrm*real(b1); b1_k(2,nn,adir) = trnrm*aimag(b1) 1265 b2_k(1,nn,adir) = trnrm*real(bv2b); b2_k(2,nn,adir) = trnrm*aimag(bv2b) 1266 m1_k(1,nn,adir) = trnrm*real(m1); m1_k(2,nn,adir) = trnrm*aimag(m1) 1267 m1_mu_k(1,nn,adir) = trnrm*real(m1_mu); m1_mu_k(2,nn,adir) = trnrm*aimag(m1_mu) 1268 m2_k(1,nn,adir) = trnrm*real(mv2b); m2_k(2,nn,adir) = trnrm*aimag(mv2b) 1269 m2_mu_k(1,nn,adir) = trnrm*real(mv2b_mu); m2_mu_k(2,nn,adir) = trnrm*aimag(mv2b_mu) 1270 1271 end do !nn 1272 end do !adir 1273 1274 ABI_FREE(bra) 1275 ABI_FREE(ket) 1276 ABI_FREE(svectoutb) 1277 ABI_FREE(svectoutg) 1278 ABI_FREE(vectout) 1279 call pawcprj_free(cwaveprj) 1280 ABI_FREE(cwaveprj) 1281 1282 end subroutine orbmag_vv_k
ABINIT/tt_me [ Functions ]
NAME
tt_me
FUNCTION
Onsite part of matrix element <u_n|a_ij|u_m>
COPYRIGHT
Copyright (C) 2003-2024 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
aij(dtset%natom,lmn2max,ndij)=(complex)scalar ij couplings atindx(natom)=index table for atoms (see gstate.f) bcp(dtset%natom,dtset%nspinor)<type(pawcprj_type)>=bra side cprj <p|bra> dtset <type(dataset_type)>=all input variables for this dataset kcp(dtset%natom,dtset%nspinor)<type(pawcprj_type)>=ket side cprj <p|ket> lmn2max=max value of lmn2 over all psps pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data ndij=spin channels in dij
OUTPUT
tt=(complex) computed matrix element
SIDE EFFECTS
TODO
NOTES
computes on-site \sum_{Rij}<u|p_i>aij<p_j|u> for generic aij input
SOURCE
1629 subroutine tt_me(aij,atindx,bcp,dtset,kcp,lmn2max,ndij,pawtab,tt) 1630 1631 !Arguments ------------------------------------ 1632 !scalars 1633 integer,intent(in) :: lmn2max,ndij 1634 complex(dpc),intent(out) :: tt 1635 type(dataset_type),intent(in) :: dtset 1636 1637 !arrays 1638 integer,intent(in) :: atindx(dtset%natom) 1639 complex(dpc),intent(in) :: aij(dtset%natom,lmn2max,ndij) 1640 type(pawcprj_type),intent(in) :: bcp(dtset%natom,dtset%nspinor),kcp(dtset%natom,dtset%nspinor) 1641 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1642 1643 !Local variables ------------------------- 1644 !scalars 1645 integer :: iat,iatom,isp,itypat,ilmn,jlmn,klmn 1646 complex(dpc) :: cpi,cpj,dij 1647 1648 !arrays 1649 1650 !-------------------------------------------------------------------- 1651 1652 tt = czero 1653 do iat = 1, dtset%natom 1654 iatom = atindx(iat) 1655 itypat=dtset%typat(iat) 1656 do isp = 1, dtset%nspinor 1657 do ilmn = 1, pawtab(itypat)%lmn_size 1658 do jlmn = 1, pawtab(itypat)%lmn_size 1659 klmn=MATPACK(ilmn,jlmn) 1660 ! in ndij = 4 case, isp 1 delivers up-up, isp 2 delivers down-down 1661 dij = aij(iatom,klmn,isp) 1662 cpi = CMPLX(bcp(iatom,isp)%cp(1,ilmn),bcp(iatom,isp)%cp(2,ilmn)) 1663 cpj = CMPLX(kcp(iatom,isp)%cp(1,jlmn),kcp(iatom,isp)%cp(2,jlmn)) 1664 ! see note at top of file near definition of MATPACK macro 1665 if (ilmn .GT. jlmn) dij = CONJG(dij) 1666 ! note use of CONJG(cpi), because cpi is from the bra side cprj 1667 tt = tt + CONJG(cpi)*dij*cpj 1668 if (ndij == 4) then 1669 if (isp == 1) then 1670 dij = aij(iatom,klmn,3) ! up-down 1671 ! D^ss'_ij=D^s's_ji^* 1672 if (ilmn .GT. jlmn) dij = CONJG(aij(iatom,klmn,4)) 1673 cpi = CMPLX(bcp(iatom,1)%cp(1,ilmn),bcp(iatom,1)%cp(2,ilmn)) 1674 cpj = CMPLX(kcp(iatom,2)%cp(1,jlmn),kcp(iatom,2)%cp(2,jlmn)) 1675 else 1676 dij = aij(iatom,klmn,4) ! down-up 1677 ! D^ss'_ij=D^s's_ji^* 1678 if (ilmn .GT. jlmn) dij = CONJG(aij(iatom,klmn,3)) 1679 cpi = CMPLX(bcp(iatom,2)%cp(1,ilmn),bcp(iatom,2)%cp(2,ilmn)) 1680 cpj = CMPLX(kcp(iatom,1)%cp(1,jlmn),kcp(iatom,1)%cp(2,jlmn)) 1681 end if 1682 tt = tt + CONJG(cpi)*cpj*dij 1683 end if 1684 end do !jlmn 1685 end do !ilmn 1686 end do ! isp 1687 end do !iat 1688 1689 end subroutine tt_me
ABINIT/txt_me [ Functions ]
NAME
txt_me
FUNCTION
Onsite part of matrix element <u_n|dp>a_ij<dp|u_m>
COPYRIGHT
Copyright (C) 2003-2024 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
aij(dtset%natom,lmn2max,ndij)=(complex)scalar ij couplings atindx(natom)=index table for atoms (see gstate.f) bcp(dtset%natom,dtset%nspinor)<type(pawcprj_type)>=bra side cprj <p|bra> bdir=direction of bcp derivative to use dtset <type(dataset_type)>=all input variables for this dataset gdir=direction of kcp derivative to use kcp(dtset%natom,dtset%nspinor)<type(pawcprj_type)>=ket side cprj <p|ket> lmn2max=max value of lmn2 over all psps pawtab(dtset%ntypat) <type(pawtab_type)>=paw tabulated starting data ndij=spin channels in dij
OUTPUT
txt=(complex) computed matrix element
OUTPUT
SIDE EFFECTS
TODO
NOTES
computes on-site \sum_{Rij}<u|d_bdir p_i>aij<d_gdir p_j|u> for generic aij input
SOURCE
1529 subroutine txt_me(aij,atindx,bcp,bdir,dtset,gdir,kcp,lmn2max,ndij,pawtab,txt) 1530 1531 !Arguments ------------------------------------ 1532 !scalars 1533 integer,intent(in) :: bdir,gdir,lmn2max,ndij 1534 complex(dpc),intent(out) :: txt 1535 type(dataset_type),intent(in) :: dtset 1536 1537 !arrays 1538 integer,intent(in) :: atindx(dtset%natom) 1539 complex(dpc),intent(in) :: aij(dtset%natom,lmn2max,ndij) 1540 type(pawcprj_type),intent(in) :: bcp(dtset%natom,dtset%nspinor) 1541 type(pawcprj_type),intent(in) :: kcp(dtset%natom,dtset%nspinor) 1542 type(pawtab_type),intent(in) :: pawtab(dtset%ntypat) 1543 1544 !Local variables ------------------------- 1545 !scalars 1546 integer :: iat,iatom,itypat,ilmn,isp,jlmn,klmn 1547 complex(dpc) :: dcpi,dcpj,dij 1548 1549 !arrays 1550 1551 !-------------------------------------------------------------------- 1552 1553 txt = czero 1554 do iat = 1, dtset%natom 1555 itypat=dtset%typat(iat) 1556 iatom = atindx(iat) 1557 do isp = 1, dtset%nspinor 1558 do ilmn = 1, pawtab(itypat)%lmn_size 1559 do jlmn = 1, pawtab(itypat)%lmn_size 1560 klmn = MATPACK(ilmn,jlmn) 1561 ! in ndij = 4 case, isp 1 delivers up-up, isp 2 delivers down-down 1562 dij = aij(iatom,klmn,isp) 1563 ! see note at top of file near definition of MATPACK macro 1564 if (ilmn .GT. jlmn) dij = CONJG(dij) 1565 dcpi = CMPLX(bcp(iatom,isp)%dcp(1,bdir,ilmn),bcp(iatom,isp)%dcp(2,bdir,ilmn)) 1566 dcpj = CMPLX(kcp(iatom,isp)%dcp(1,gdir,jlmn),kcp(iatom,isp)%dcp(2,gdir,jlmn)) 1567 txt = txt + CONJG(dcpi)*dcpj*dij 1568 if (ndij == 4) then 1569 if (isp == 1) then 1570 dij = aij(iatom,klmn,3) ! up-down 1571 ! D^ss'_ij=D^s's_ji^* 1572 if (ilmn .GT. jlmn) dij = CONJG(aij(iatom,klmn,4)) 1573 dcpi = CMPLX(bcp(iatom,1)%dcp(1,bdir,ilmn),bcp(iatom,1)%dcp(2,bdir,ilmn)) 1574 dcpj = CMPLX(kcp(iatom,2)%dcp(1,gdir,jlmn),kcp(iatom,2)%dcp(2,gdir,jlmn)) 1575 else 1576 dij = aij(iatom,klmn,4) ! down-up 1577 ! D^ss'_ij=D^s's_ji^* 1578 if (ilmn .GT. jlmn) dij = CONJG(aij(iatom,klmn,3)) 1579 dcpi = CMPLX(bcp(iatom,2)%dcp(1,bdir,ilmn),bcp(iatom,2)%dcp(2,bdir,ilmn)) 1580 dcpj = CMPLX(kcp(iatom,1)%dcp(1,gdir,jlmn),kcp(iatom,1)%dcp(2,gdir,jlmn)) 1581 end if 1582 txt = txt + CONJG(dcpi)*dcpj*dij 1583 end if 1584 end do !jlmn 1585 end do !ilmn 1586 end do ! isp 1587 end do !iat 1588 1589 end subroutine txt_me