Line data Source code
1 : MODULE m_crystalfield
2 :
3 : !------------------------------------------------------------------------------
4 : !
5 : ! MODULE: m_crystalfield
6 : !
7 : !> @author
8 : !> Henning Janßen
9 : !
10 : ! DESCRIPTION:
11 : !> -calculates the crystal-field-contrbution for the local hamiltonian
12 : !
13 : !------------------------------------------------------------------------------
14 : USE m_juDFT
15 : USE m_types
16 : USE m_constants
17 : USE m_trapz
18 : USE m_sgml
19 : USE m_rotMMPmat
20 :
21 : IMPLICIT NONE
22 :
23 : CONTAINS
24 :
25 0 : SUBROUTINE crystal_field(atoms,gfinp,input,noco,nococonv,greensfImagPart,v,ef,hub1data)
26 :
27 : !calculates the crystal-field matrix for the local hamiltonian
28 :
29 : TYPE(t_greensfImagPart), INTENT(IN) :: greensfImagPart
30 : TYPE(t_atoms), INTENT(IN) :: atoms
31 : TYPE(t_gfinp), INTENT(IN) :: gfinp
32 : TYPE(t_input), INTENT(IN) :: input
33 : TYPE(t_noco), INTENT(IN) :: noco
34 : TYPE(t_nococonv), INTENT(IN) :: nococonv
35 : TYPE(t_potden), INTENT(IN) :: v !LDA+U potential (should be removed from h_loc)
36 : REAL, INTENT(IN) :: ef
37 : TYPE(t_hub1data), INTENT(INOUT) :: hub1data
38 :
39 : !-Local Scalars
40 : INTEGER i_gf,l,nType,jspin,m,mp,ie,i_hia,i_u,isp,i_elem
41 : REAL tr,del,eb
42 : COMPLEX vso
43 : LOGICAL, PARAMETER :: l_correctMinus = .FALSE.
44 : REAL, PARAMETER :: excTolerance = 0.2/hartree_to_ev_const
45 : !-Local Arrays
46 0 : REAL :: h_loc(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,atoms%n_hia,input%jspins)
47 : REAL :: ex(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
48 : REAL :: shift(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
49 0 : REAL :: integrand(gfinp%ne), norm(gfinp%ne)
50 0 : COMPLEX, ALLOCATABLE :: imag(:,:,:)
51 0 : COMPLEX, ALLOCATABLE :: potmmpmat(:,:,:)
52 :
53 :
54 0 : ALLOCATE(potmmpmat(-lmaxU_const:lmaxU_const, -lmaxU_const:lmaxU_const, input%jspins))
55 0 : ALLOCATE(imag(gfinp%ne,-lmaxU_const:lmaxU_const, -lmaxU_const:lmaxU_const),source=cmplx_0)
56 :
57 0 : h_loc = 0.0
58 0 : DO i_hia = 1, atoms%n_hia
59 :
60 0 : l = atoms%lda_u(atoms%n_u+i_hia)%l
61 0 : nType = atoms%lda_u(atoms%n_u+i_hia)%atomType
62 :
63 0 : i_gf = gfinp%hiaElem(i_hia)
64 0 : i_elem = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=.TRUE., l_kresolved_int=.false.)
65 : !---------------------------------------------------------
66 : ! Perform the integration
67 : !---------------------------------------------------------
68 : ! \int_{E_b}^{E_c} dE E * N_LL'(E)
69 : !---------------------------------------------------------
70 : ! N_LL' is the l projected density of states and
71 : ! connected to the imaginary part of the greens function
72 : ! with a factor -1/pi
73 : !---------------------------------------------------------
74 0 : CALL gfinp%eMesh(ef,del,eb)
75 0 : DO jspin = 1, input%jspins
76 : !Use the same cutoffs as in the kramer kronigs integration
77 0 : norm = 0.0
78 0 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,.TRUE.)/(3.0-input%jspins)
79 0 : DO m = -l, l
80 0 : DO mp = -l, l
81 0 : integrand = 0.0
82 0 : DO ie = 1, gfinp%ne
83 0 : integrand(ie) = -1.0/pi_const * ((ie-1) * del+eb) * imag(ie,m,mp)
84 0 : IF(m.EQ.mp) norm(ie) = norm(ie) -1.0/pi_const * imag(ie,m,mp)
85 : ENDDO
86 0 : h_loc(m,mp,i_hia,jspin) = trapz(integrand,del,gfinp%ne)
87 : ENDDO
88 : ENDDO
89 : ENDDO
90 : #ifdef CPP_DEBUG
91 : WRITE(*,*) "UP"
92 : WRITE(*,"(7f7.3)") h_loc(-3:3,-3:3,i_hia,1)
93 : IF(input%jspins.EQ.2) THEN
94 : WRITE(*,*) "DOWN"
95 : WRITE(*,"(7f7.3)") h_loc(-3:3,-3:3,i_hia,2)
96 : ENDIF
97 : #endif
98 : !Remove LDA+U potential
99 0 : i_u = atoms%n_u+i_hia !position in the v%mmpmat array
100 : !Rotate the occupation matrix into the global frame in real-space
101 0 : IF(noco%l_noco) THEN
102 0 : potmmpmat = rotMMPmat(v%mmpmat(:,:,i_u,:input%jspins),0.0,-nococonv%beta(nType),-nococonv%alph(nType),l)
103 0 : ELSE IF(noco%l_soc) THEN
104 0 : potmmpmat = rotMMPmat(v%mmpmat(:,:,i_u,:input%jspins),0.0,-nococonv%theta,-nococonv%phi,l)
105 : ENDIF
106 0 : DO jspin = 1, input%jspins
107 0 : DO m = -l, l
108 0 : DO mp = -l, l
109 0 : IF(ABS(potmmpmat(m,mp,jspin)).GT.1e-4) THEN
110 0 : h_loc(m,mp,i_hia,jspin) = h_loc(m,mp,i_hia,jspin) - potmmpmat(m,mp,jspin)
111 : ENDIF
112 : ENDDO
113 : ENDDO
114 : ENDDO
115 0 : IF(ABS(hub1data%xi(i_hia)).GT.1e-12) THEN
116 : !Remove SOC potential (only spin-diagonal)
117 0 : DO jspin = 1, 2
118 0 : DO m = -l, l
119 0 : DO mp = -l, l
120 0 : isp = 3-2*jspin !1,-1
121 0 : vso = CMPLX(sgml(l,m,isp,l,mp,isp),0.0)
122 0 : h_loc(m,mp,i_hia,jspin) = h_loc(m,mp,i_hia,jspin) - REAL(vso)/2.0 * hub1data%xi(i_hia)/hartree_to_ev_const
123 : ENDDO
124 : ENDDO
125 : ENDDO
126 : ENDIF
127 : #ifdef CPP_DEBUG
128 : WRITE(*,*) "UP-REMOVED"
129 : WRITE(*,"(7f7.3)") h_loc(-3:3,-3:3,i_hia,1)
130 : IF(input%jspins.EQ.2) THEN
131 : WRITE(*,*) "DOWN-REMOVED"
132 : WRITE(*,"(7f7.3)") h_loc(-3:3,-3:3,i_hia,2)
133 : ENDIF
134 : #endif
135 0 : IF(input%jspins.EQ.2) THEN
136 : ex = 0.0
137 0 : DO m= -l, l
138 0 : DO mp = -l, l
139 0 : ex(m,mp) = h_loc(m,mp,i_hia,1)-h_loc(m,mp,i_hia,2)
140 : ENDDO
141 : ENDDO
142 : #ifdef CPP_DEBUG
143 : WRITE(*,*) "Exchange (eV)"
144 : WRITE(*,"(7f7.3)") ex(-3:3,-3:3)*hartree_to_ev_const*0.5
145 : #endif
146 : !------------------------------------------------------------------------------------
147 : ! If states move close to the cutoff we get some shift in the results on the diagonal
148 : ! The reason for this is a bit unclear we remove these results and replace them
149 : ! with either the -m or corresponding opposite spin result (only diagonal)
150 : !------------------------------------------------------------------------------------
151 : !Shift m=0 exchange to 0
152 : IF(l_correctMinus)THEN
153 : IF(ABS(ex(0,0)).LT.excTolerance) THEN
154 : !There is no big discrepancy between m=0 spin up/down => simply eleminate the exchange
155 : DO m = -l, l
156 : h_loc(m,m,i_hia,1) = h_loc(m,m,i_hia,1) - ex(0,0)/2.0
157 : h_loc(m,m,i_hia,2) = h_loc(m,m,i_hia,2) + ex(0,0)/2.0
158 : ENDDO
159 : ELSE
160 : !There is a big discrepancy due to numerical problems => Take the spin up part
161 : h_loc(0,0,i_hia,2) = h_loc(0,0,i_hia,1)
162 : ENDIF
163 :
164 : !Recalculate exchange
165 : ex = 0.0
166 : DO m= -l, l
167 : DO mp = -l, l
168 : ex(m,mp) = h_loc(m,mp,i_hia,1)-h_loc(m,mp,i_hia,2)
169 : ENDDO
170 : ENDDO
171 : #ifdef CPP_DEBUG
172 : WRITE(*,*) "Exchange shifted (eV)"
173 : WRITE(*,"(7f7.3)") ex(-3:3,-3:3)*hartree_to_ev_const*0.5
174 : #endif
175 : !Calculate shifts for numerically "troubled" elements
176 : shift=0.0
177 : DO m = -l, -1
178 : !Normal leftover from SOC+numerical problems
179 : IF(ABS(ex(m,m)).LT.excTolerance) CYCLE
180 : shift(m,m) = ex(m,m) + ex(-m,-m)
181 : ENDDO
182 : h_loc(:,:,i_hia,2) = h_loc(:,:,i_hia,2) + shift
183 :
184 : #ifdef CPP_DEBUG
185 : !Recalculate differences for verification
186 : ex = 0.0
187 : DO m= -l, l
188 : DO mp = -l, l
189 : ex(m,mp) = h_loc(m,mp,i_hia,1)-h_loc(m,mp,i_hia,2)
190 : ENDDO
191 : ENDDO
192 :
193 : WRITE(*,*) "Exchange after Correction (eV)"
194 : WRITE(*,"(7f7.3)") ex(-3:3,-3:3)*hartree_to_ev_const*0.5
195 : #endif
196 : ENDIF
197 : ENDIF
198 :
199 : !Average over spins
200 0 : hub1data%ccfmat(i_hia,:,:) = 0.0
201 0 : DO m = -l, l
202 0 : DO mp = -l, l
203 0 : hub1data%ccfmat(i_hia,m,mp) = SUM(h_loc(m,mp,i_hia,:))/2.0
204 : !For jspins.EQ.1 we need to take care of the fact that the spin-orbit coupling is opposite in spin 1/2
205 0 : IF(input%jspins.EQ.1) hub1data%ccfmat(i_hia,m,mp) = (h_loc(m,mp,i_hia,1)+h_loc(-m,-mp,i_hia,1))/2.0
206 : ENDDO
207 : ENDDO
208 : #ifdef CPP_DEBUG
209 : WRITE(*,*) "Average"
210 : WRITE(*,"(7f7.3)") hub1data%ccfmat(i_hia,-3:3,-3:3)
211 : #endif
212 :
213 : !-----------------------------------
214 : ! Symmetrize Matrix
215 : ! Delta_CF(m,mp) = Delta_CF(-m,-mp)
216 : !-----------------------------------
217 0 : IF(input%jspins.EQ.2) THEN
218 0 : DO m = -l, l
219 0 : DO mp = -l, l
220 0 : hub1data%ccfmat(i_hia,m,mp) = (hub1data%ccfmat(i_hia,m,mp)+hub1data%ccfmat(i_hia,-m,-mp))/2.0
221 0 : hub1data%ccfmat(i_hia,-m,-mp) = hub1data%ccfmat(i_hia,m,mp)
222 : ENDDO
223 : ENDDO
224 : ENDIF
225 : #ifdef CPP_DEBUG
226 : WRITE(*,*) "Average symmetrized"
227 : WRITE(*,"(7f7.3)") hub1data%ccfmat(i_hia,-3:3,-3:3)
228 : #endif
229 :
230 0 : tr = 0.0
231 : !calculate the trace
232 0 : DO m = -l, l
233 0 : tr = tr + hub1data%ccfmat(i_hia,m,m)
234 : ENDDO
235 : #ifdef CPP_DEBUG
236 : WRITE(*,*) "TRACE"
237 : WRITE(*,"(2f7.3)") tr, tr/(2*l+1)
238 : #endif
239 : !Remove trace
240 0 : DO m = -l, l
241 0 : hub1data%ccfmat(i_hia,m,m) = hub1data%ccfmat(i_hia,m,m) - tr/(2*l+1)
242 : ENDDO
243 :
244 : #ifdef CPP_DEBUG
245 : WRITE(*,*) "TRACELESS ccf (eV)"
246 : WRITE(*,"(7f7.3)") hub1data%ccfmat(i_hia,-3:3,-3:3)*hartree_to_ev_const
247 : #endif
248 : ENDDO
249 0 : END SUBROUTINE crystal_field
250 :
251 : END MODULE m_crystalfield
|
