Line data Source code
1 : ! Copyright (c) 2018 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
2 : ! This file is part of FLEUR and available as free software under the conditions
3 : ! of the MIT license as expressed in the LICENSE file in more detail.
4 : !--------------------------------------------------------------------------------
5 :
6 : MODULE m_magDiMom
7 :
8 : CONTAINS
9 :
10 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
11 : !
12 : ! This subroutine calculates intraatomic magnetic dipole moments.
13 : !
14 : ! GM'2018
15 : !
16 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
17 :
18 0 : SUBROUTINE magDiMom(sym,input,atoms,sphhar,noco,nococonv,l_fmpl2,rho,magDipoles,elecDipoles)
19 :
20 : USE m_constants
21 : USE m_types
22 : USE m_juDFT
23 : USE m_rotdenmat
24 : USE m_lattHarmsSphHarmsConv
25 : USE m_gaunt
26 : USE m_intgr
27 :
28 : IMPLICIT NONE
29 : TYPE(t_sym), INTENT(IN) :: sym
30 : TYPE(t_input), INTENT(IN) :: input
31 : TYPE(t_atoms), INTENT(IN) :: atoms
32 : TYPE(t_sphhar), INTENT(IN) :: sphhar
33 : TYPE(t_noco), INTENT(IN) :: noco
34 : TYPE(t_nococonv), INTENT(IN) :: nococonv
35 : REAL, INTENT(IN) :: rho(:,0:,:,:) ! if l_fmpl last dimension is 4, otherwise 2.
36 :
37 : LOGICAL, INTENT(IN) :: l_fmpl2
38 : REAL, INTENT(INOUT) :: magDipoles(:,:)
39 : REAL, INTENT(INOUT) :: elecDipoles(:,:)
40 :
41 :
42 0 : REAL, ALLOCATABLE :: inRho(:,:,:,:)
43 0 : COMPLEX, ALLOCATABLE :: rhoSphHarms(:,:,:,:), rhoTempSphHarms(:,:,:,:)
44 0 : COMPLEX, ALLOCATABLE :: rhoSphHarmsR(:,:,:)
45 :
46 : INTEGER :: iType, ilh, l, m, lm, lp, mp, lmp, i
47 : REAL :: theta, phi, cdn11, cdn22
48 : REAL :: magDipole(3), myCharge, elecDipole(3)
49 : REAL :: constA
50 : COMPLEX :: gauntA, gauntB, gauntC
51 : COMPLEX :: cdn21
52 :
53 0 : IF(input%jspins.EQ.1) RETURN
54 0 : IF(.NOT.noco%l_noco) RETURN
55 :
56 : !---> calculate the charge and magnetization density in the muffin tins
57 0 : ALLOCATE(inRho(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,4))
58 0 : DO iType = 1,atoms%ntype
59 0 : IF (.NOT.l_fmpl2) THEN
60 0 : theta = nococonv%beta(iType)
61 0 : phi = nococonv%alph(iType)
62 0 : inRho(:,:,iType,1) = rho(:,:,iType,1) + rho(:,:,iType,2)
63 0 : inRho(:,:,iType,2) = rho(:,:,iType,1) - rho(:,:,iType,2)
64 0 : inRho(:,:,iType,3) = inRho(:,:,iType,2) * SIN(phi)*SIN(theta)
65 0 : inRho(:,:,iType,4) = inRho(:,:,iType,2) * COS(theta)
66 0 : inRho(:,:,iType,2) = inRho(:,:,iType,2) * COS(phi)*SIN(theta)
67 : ELSE
68 0 : DO ilh = 0,sphhar%nlh(sym%ntypsy(atoms%firstAtom(iType)))
69 0 : DO i = 1,atoms%jri(iType)
70 :
71 0 : cdn11 = rho(i,ilh,iType,1)
72 0 : cdn22 = rho(i,ilh,iType,2)
73 0 : cdn21 = CMPLX(rho(i,ilh,iType,3),rho(i,ilh,iType,4))
74 0 : CALL rot_den_mat(nococonv%alph(iType),nococonv%beta(iType),cdn11,cdn22,cdn21)
75 0 : inRho(i,ilh,iType,1) = cdn11 + cdn22
76 0 : inRho(i,ilh,iType,2) = 2.0 * REAL(cdn21)
77 : ! Note: The minus sign in the following line is temporary to adjust for differences in the offdiagonal
78 : ! part of the density between this fleur version and ancient (v0.26) fleur.
79 0 : inRho(i,ilh,iType,3) = -2.0 * AIMAG(cdn21)
80 0 : inRho(i,ilh,iType,4) = cdn11 - cdn22
81 : END DO
82 : END DO
83 : END IF
84 : END DO
85 :
86 0 : ALLOCATE (rhoSphHarms(atoms%jmtd,(atoms%lmaxd+1)**2,atoms%ntype,4))
87 0 : ALLOCATE (rhoTempSphHarms(atoms%jmtd,(atoms%lmaxd+1)**2,atoms%ntype,4))
88 0 : ALLOCATE (rhoSphHarmsR(atoms%jmtd,(atoms%lmaxd+1)**2,atoms%ntype))
89 :
90 0 : rhoSphHarms = CMPLX(0.0,0.0)
91 0 : DO i = 1, 4
92 0 : DO iType = 1, atoms%ntype
93 0 : CALL lattHarmsRepToSphHarms(sym,atoms,sphhar,iType,inRho(:,0:,iType,i),rhoSphHarms(:,:,iType,i))
94 : END DO
95 : END DO
96 :
97 : ! electric dipole moment (start)
98 0 : DO iType = 1, atoms%ntype
99 0 : DO i = 1, atoms%jri(iType)
100 0 : rhoSphHarmsR(i,:,iType) = rhoSphHarms(i,:,iType,1) * atoms%rmsh(i,iType)
101 : END DO
102 : END DO
103 :
104 0 : constA = SQRT(2.0*pi_const/3.0)
105 0 : rhoTempSphHarms = CMPLX(0.0,0.0)
106 0 : DO iType = 1, atoms%ntype
107 0 : DO lp = 0, MIN(2,atoms%lmax(iType))
108 0 : DO mp = -lp, lp
109 0 : DO l = 0, MIN(2,atoms%lmax(iType))
110 0 : DO m = -l, l
111 : !note 1: For refinement maybe I could make use of selection rules.
112 : !note 2: ls for r^\hat is 1, ms is -1..1.
113 0 : gauntA = gaunt1(lp,1,l,mp,-1,m,atoms%lmaxd)
114 0 : gauntB = gaunt1(lp,1,l,mp,0,m,atoms%lmaxd)
115 0 : gauntC = gaunt1(lp,1,l,mp,1,m,atoms%lmaxd)
116 0 : lmp = lp*(lp+1) + mp + 1
117 0 : lm = l*(l+1) + m + 1
118 0 : DO i = 1, atoms%jri(iType)
119 : rhoTempSphHarms(i,lmp,iType,2) = rhoTempSphHarms(i,lmp,iType,2) +&
120 0 : constA*(gauntA-gauntC)*rhoSphHarmsR(i,lm,iType)
121 : rhoTempSphHarms(i,lmp,iType,3) = rhoTempSphHarms(i,lmp,iType,3) +&
122 0 : constA*CMPLX(0.0,1.0)*(gauntA+gauntC)*rhoSphHarmsR(i,lm,iType)
123 : rhoTempSphHarms(i,lmp,iType,4) = rhoTempSphHarms(i,lmp,iType,4) +&
124 0 : constA*SQRT(2.0)*gauntB*rhoSphHarmsR(i,lm,iType)
125 : END DO
126 : END DO
127 : END DO
128 : END DO
129 : END DO
130 : END DO
131 :
132 0 : elecDipole = 0.0
133 0 : DO iType = 1, atoms%ntype
134 0 : CALL intgr3(REAL(rhoTempSphHarms(:,1,iType,2)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),elecDipole(1))
135 0 : CALL intgr3(REAL(rhoTempSphHarms(:,1,iType,3)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),elecDipole(2))
136 0 : CALL intgr3(REAL(rhoTempSphHarms(:,1,iType,4)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),elecDipole(3))
137 0 : elecDipoles(:,iType) = elecDipole(:) * sfp_const
138 : END DO
139 :
140 : ! electric dipole moment (end)
141 :
142 : ! WRITE(7534,*) '===================================================='
143 0 : DO iType = 1, atoms%ntype
144 0 : DO l = 0, 2
145 0 : DO m = -l, l
146 0 : magDipole(:) = 0.0
147 : myCharge = 0.0
148 0 : lm = l*(l+1) + m + 1
149 0 : CALL intgr3(REAL(rhoSphHarms(:,lm,iType,1)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),myCharge)
150 0 : CALL intgr3(REAL(rhoSphHarms(:,lm,iType,2)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(1))
151 0 : CALL intgr3(REAL(rhoSphHarms(:,lm,iType,3)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(2))
152 0 : CALL intgr3(REAL(rhoSphHarms(:,lm,iType,4)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(3))
153 : ! WRITE(7534,'(3i7,4e24.8)') iType, l, m, myCharge, magDipole(:)
154 : END DO
155 : END DO
156 0 : rhoSphHarms(:,:,iType,1) = CMPLX(0.0,0.0)
157 : END DO
158 :
159 0 : constA = SQRT(2.0*pi_const/3.0)
160 0 : DO iType = 1, atoms%ntype
161 0 : DO lp = 0, MIN(2,atoms%lmax(iType))
162 0 : DO mp = -lp, lp
163 0 : DO l = 0, MIN(2,atoms%lmax(iType))
164 0 : DO m = -l, l
165 : !note 1: For refinement maybe I could make use of selection rules.
166 : !note 2: ls for r^\hat is 1, ms is -1..1.
167 0 : gauntA = gaunt1(lp,1,l,mp,-1,m,atoms%lmaxd)
168 0 : gauntB = gaunt1(lp,1,l,mp,0,m,atoms%lmaxd)
169 0 : gauntC = gaunt1(lp,1,l,mp,1,m,atoms%lmaxd)
170 0 : lmp = lp*(lp+1) + mp + 1
171 0 : lm = l*(l+1) + m + 1
172 0 : DO i = 1, atoms%jri(iType)
173 : rhoSphHarms(i,lmp,iType,1) = rhoSphHarms(i,lmp,iType,1) +&
174 0 : constA*(gauntA-gauntC)*rhoSphHarms(i,lm,iType,2)
175 : rhoSphHarms(i,lmp,iType,1) = rhoSphHarms(i,lmp,iType,1) +&
176 0 : constA*CMPLX(0.0,1.0)*(gauntA+gauntC)*rhoSphHarms(i,lm,iType,3)
177 : rhoSphHarms(i,lmp,iType,1) = rhoSphHarms(i,lmp,iType,1) +&
178 0 : constA*SQRT(2.0)*gauntB*rhoSphHarms(i,lm,iType,4)
179 : END DO
180 : END DO
181 : END DO
182 : END DO
183 : END DO
184 : END DO
185 :
186 0 : rhoTempSphHarms = CMPLX(0.0,0.0)
187 0 : DO iType = 1, atoms%ntype
188 0 : DO lp = 0, MIN(2,atoms%lmax(iType))
189 0 : DO mp = -lp, lp
190 0 : DO l = 0, MIN(2,atoms%lmax(iType))
191 0 : DO m = -l, l
192 : !note 1: For refinement maybe I could make use of selection rules.
193 : !note 2: ls for r^\hat is 1, ms is -1..1.
194 0 : gauntA = gaunt1(lp,1,l,mp,-1,m,atoms%lmaxd)
195 0 : gauntB = gaunt1(lp,1,l,mp,0,m,atoms%lmaxd)
196 0 : gauntC = gaunt1(lp,1,l,mp,1,m,atoms%lmaxd)
197 0 : lmp = lp*(lp+1) + mp + 1
198 0 : lm = l*(l+1) + m + 1
199 0 : DO i = 1, atoms%jri(iType)
200 : rhoTempSphHarms(i,lmp,iType,2) = rhoTempSphHarms(i,lmp,iType,2) +&
201 0 : constA*(gauntA-gauntC)*rhoSphHarms(i,lm,iType,1)
202 : rhoTempSphHarms(i,lmp,iType,3) = rhoTempSphHarms(i,lmp,iType,3) +&
203 0 : constA*CMPLX(0.0,1.0)*(gauntA+gauntC)*rhoSphHarms(i,lm,iType,1)
204 : rhoTempSphHarms(i,lmp,iType,4) = rhoTempSphHarms(i,lmp,iType,4) +&
205 0 : constA*SQRT(2.0)*gauntB*rhoSphHarms(i,lm,iType,1)
206 : END DO
207 : END DO
208 : END DO
209 : END DO
210 : END DO
211 : END DO
212 :
213 0 : DO iType = 1, atoms%ntype
214 0 : DO l = 0, atoms%lmax(iType)
215 0 : DO m = -l, l
216 0 : lm = l*(l+1) + m + 1
217 0 : DO i = 1, atoms%jri(iType)
218 0 : rhoSphHarms(i,lm,iType,2) = rhoSphHarms(i,lm,iType,2) - 3.0 * rhoTempSphHarms(i,lm,iType,2)
219 0 : rhoSphHarms(i,lm,iType,3) = rhoSphHarms(i,lm,iType,3) - 3.0 * rhoTempSphHarms(i,lm,iType,3)
220 0 : rhoSphHarms(i,lm,iType,4) = rhoSphHarms(i,lm,iType,4) - 3.0 * rhoTempSphHarms(i,lm,iType,4)
221 : END DO
222 : END DO
223 : END DO
224 : END DO
225 :
226 0 : DO iType = 1, atoms%ntype
227 0 : DO i = 1, atoms%jri(iType)
228 0 : IF (ANY(AIMAG(rhoSphHarms(i,1,iType,:)).GT.1.0e-11)) THEN
229 0 : WRITE(oUnit,*) 'imaginary part too large!'
230 : END IF
231 : END DO
232 : END DO
233 :
234 0 : magDipole = 0.0
235 0 : DO iType = 1, atoms%ntype
236 0 : CALL intgr3(REAL(rhoSphHarms(:,1,iType,2)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(1))
237 0 : CALL intgr3(REAL(rhoSphHarms(:,1,iType,3)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(2))
238 0 : CALL intgr3(REAL(rhoSphHarms(:,1,iType,4)),atoms%rmsh(:,iType),atoms%dx(iType),atoms%jri(iType),magDipole(3))
239 0 : magDipoles(:,iType) = magDipole(:) * sfp_const
240 : END DO
241 :
242 0 : END SUBROUTINE magDiMom
243 :
244 : END MODULE m_magDiMom
|
