Line data Source code
1 : MODULE m_excSplitting
2 :
3 : use m_types
4 : use m_constants
5 : use m_trapz
6 : use m_xmlOutput
7 : use m_intgr
8 :
9 : IMPLICIT NONE
10 :
11 : CONTAINS
12 :
13 21 : SUBROUTINE excSplitting(gfinp,atoms,input,scalarGF,greensfImagPart,ef,vTot)
14 :
15 : TYPE(t_gfinp), INTENT(IN) :: gfinp
16 : TYPE(t_atoms), INTENT(IN) :: atoms
17 : TYPE(t_input), INTENT(IN) :: input
18 : TYPE(t_scalarGF), INTENT(IN) :: scalarGF(:)
19 : TYPE(t_greensfImagPart), INTENT(IN) :: greensfImagPart
20 : REAL, INTENT(IN) :: ef
21 : type(t_potden), intent(in) :: vTot
22 :
23 : INTEGER :: i_gf,i_elem,i_elemLO,ispin,m,l,atomType,nLO,iLO,iLOp
24 : LOGICAL :: l_sphavg
25 : REAL :: excSplit,del,delta
26 21 : REAL, ALLOCATABLE :: eMesh(:)
27 21 : COMPLEX, ALLOCATABLE :: imag(:,:,:)
28 21 : REAL, ALLOCATABLE :: intCOM(:,:), intNorm(:,:)
29 21 : REAL, ALLOCATABLE :: bxc(:)
30 : CHARACTER(LEN=20) :: attributes(4)
31 21 : CALL openXMLElementNoAttributes('onSiteExchangeSplitting')
32 :
33 : !Calculate the onsite exchange splitting from Bxc
34 17452 : allocate(bxc(atoms%jmtd), source=0.0)
35 47 : do atomType = 1, atoms%ntype
36 : !Get the Bxc part of the potential
37 : !L=0 of potential has an additional rescaling of r/sqrt(4pi)
38 : !The sqrt(4pi) is the part of the integral over the angular part
39 : bxc(:atoms%jri(atomType)) = (vTot%mt(:atoms%jri(atomType),0,atomType,2) - vTot%mt(:atoms%jri(atomType),0,atomType,1))/2.0 &
40 21348 : * atoms%rmsh(:atoms%jri(atomType),atomType)
41 :
42 26 : CALL intgr3(bxc,atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),delta)
43 :
44 130 : attributes = ''
45 26 : WRITE(attributes(1),'(i0)') atomType
46 26 : WRITE(attributes(2),'(f12.7)') delta * hartree_to_ev_const
47 26 : WRITE(attributes(3),'(a2)') 'eV'
48 : CALL writeXMLElementForm('bxcIntegral',['atomType','Delta ','units '],&
49 125 : attributes(:3),reshape([8,1,5,5,6,1,12,2],[3,2]))
50 :
51 : enddo
52 :
53 21 : IF(.NOT.gfinp%checkOnsite()) RETURN
54 :
55 21 : CALL gfinp%eMesh(ef,del=del,eMesh=eMesh)
56 :
57 84 : ALLOCATE(imag(SIZE(eMesh),-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const))
58 271326 : ALLOCATE(intNorm(SIZE(eMesh),input%jspins),source=0.0)
59 271305 : ALLOCATE(intCOM(SIZE(eMesh),input%jspins),source=0.0)
60 :
61 :
62 21 : WRITE(oUnit,9000)
63 : 9000 FORMAT(/,'Onsite Exchange Splitting (from imag. part of GF)')
64 21 : WRITE(oUnit,'(A)') '---------------------------------------------------'
65 :
66 530 : DO i_gf = 1, gfinp%n
67 :
68 509 : l = gfinp%elem(i_gf)%l
69 509 : atomType = gfinp%elem(i_gf)%atomType
70 509 : l_sphavg = gfinp%elem(i_gf)%l_sphavg
71 509 : nLO = gfinp%elem(i_gf)%countLOs(atoms)
72 : !Only onsite exchange splitting
73 509 : IF(gfinp%elem(i_gf)%isOffDiag()) CYCLE
74 33 : IF(gfinp%elem(i_gf)%l_kresolved_int) CYCLE
75 33 : IF(.NOT.gfinp%isUnique(i_gf, distinct_kresolved_int=.TRUE.)) CYCLE
76 :
77 31 : i_elem = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg)
78 31 : i_elemLO = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,lo=.TRUE.)
79 :
80 : !-------------------------------------------------
81 : ! Evaluate center of mass of the bands in question
82 : ! and take the difference between spin up/down
83 : ! We use the same cutoffs as for the kk-integration
84 : !-------------------------------------------------
85 : ! E_COM = 1/(int dE N_l(E)) * int dE E * N_l(E)
86 : !-------------------------------------------------
87 31 : excSplit = 0.0
88 93 : DO ispin = 1, input%jspins
89 813786 : intCOM = 0.0
90 813786 : intNorm = 0.0
91 62 : IF(l_sphavg) THEN
92 16760900 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,ispin,l_sphavg)
93 : ELSE
94 3175896 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,ispin,l_sphavg,imat=1)
95 3175896 : imag = imag + greensfImagPart%applyCutoff(i_elem,i_gf,ispin,l_sphavg,imat=2) * scalarGF(i_gf)%ddn(ispin,ispin)
96 12 : IF(nLO>0) THEN
97 18 : DO iLO = 1, nLO
98 2646580 : imag = imag + greensfImagPart%applyCutoff(i_elemLO,i_gf,ispin,l_sphavg,imat=1,iLO=iLO) * scalarGF(i_gf)%uulon(ilo,ispin,ispin)
99 2646580 : imag = imag + greensfImagPart%applyCutoff(i_elemLO,i_gf,ispin,l_sphavg,imat=2,iLO=iLO) * scalarGF(i_gf)%uloun(ilo,ispin,ispin)
100 2646580 : imag = imag + greensfImagPart%applyCutoff(i_elemLO,i_gf,ispin,l_sphavg,imat=3,iLO=iLO) * scalarGF(i_gf)%dulon(ilo,ispin,ispin)
101 2646580 : imag = imag + greensfImagPart%applyCutoff(i_elemLO,i_gf,ispin,l_sphavg,imat=4,iLO=iLO) * scalarGF(i_gf)%ulodn(ilo,ispin,ispin)
102 32 : DO iLOp = 1, nLO
103 3705222 : imag = imag + greensfImagPart%applyCutoff(i_elemLO,i_gf,ispin,l_sphavg,iLO=iLO,iLOp=iLOp) * scalarGF(i_gf)%uloulopn(ilo,ilop,ispin,ispin)
104 : ENDDO
105 : ENDDO
106 : ENDIF
107 : ENDIF
108 384 : DO m = -l, l
109 2077522 : intCOM(:,ispin) = intCOM(:,ispin) + eMesh*REAL(imag(:,m,m))
110 2077584 : intNorm(:,ispin) = intNorm(:,ispin) + REAL(imag(:,m,m))
111 : ENDDO
112 : excSplit = excSplit + (-1)**(ispin) * 1.0/(trapz(intNorm(:,ispin),del,SIZE(eMesh))) &
113 93 : * trapz(intCOM(:,ispin),del,SIZE(eMesh))
114 : ENDDO
115 31 : WRITE(oUnit,'(A,I4,A,I4,A,f10.4,A)') ' atom: ', atomType, ' l: ', l,&
116 62 : ' DeltaExc: ',excSplit * hartree_to_ev_const, ' eV'
117 :
118 155 : attributes = ''
119 31 : WRITE(attributes(1),'(i0)') atomType
120 31 : WRITE(attributes(2),'(i0)') l
121 31 : WRITE(attributes(3),'(f12.7)') excSplit * hartree_to_ev_const
122 31 : WRITE(attributes(4),'(a2)') 'eV'
123 : CALL writeXMLElementForm('excSplit',['atomType','l ','Delta ','units '],&
124 176 : attributes,reshape([8,1,5,5,6,1,12,2],[4,2]))
125 :
126 : ENDDO
127 21 : WRITE(oUnit,'(/)')
128 21 : CALL closeXMLElement('onSiteExchangeSplitting')
129 :
130 21 : END SUBROUTINE excSplitting
131 :
132 : END MODULE m_excSplitting
|
