Line data Source code
1 : MODULE m_tlmplm
2 : USE m_judft
3 :
4 : IMPLICIT NONE
5 :
6 : CONTAINS
7 2138 : SUBROUTINE tlmplm(n,sphhar,atoms,sym,enpara,nococonv,&
8 2138 : ilSpinPr,ilSpin,iSpinV,fmpi,v,vx,input,hub1inp,hub1data,td,ud,alpha_hybrid,one,l_dfpt,v1)
9 : ! Contruct the local potential matrices
10 : ! t_{L'L}^{\mu} = \sum_{lh} \int dV u_{l',order'}^{\mu}(r)Y_{l'}^{m'*}(\Omega)
11 : ! * V_{lh}(r)Y_{lh}(\Omega)
12 : ! * u_{l,order}^{\mu}(r)Y_{l}^{m}(\Omega)
13 : ! * i^{l-l'}
14 : ! of a real valued potential V(\bm{r}). The superindex L is defined as
15 : ! L := (l,m,order)
16 : ! with order = 0 refering to radial functions u and order = 1 denoting
17 : ! their energy derivatives (udot). This construction is not k-dependent
18 : ! and therefore executed only once each scf iteration.
19 :
20 : USE m_constants
21 : USE m_intgr, ONLY : intgr3
22 : USE m_genMTBasis
23 : USE m_tlo
24 : USE m_gaunt, ONLY: gaunt1
25 : USE m_types
26 :
27 : TYPE(t_input), INTENT(IN) :: input
28 : TYPE(t_sphhar), INTENT(IN) :: sphhar
29 : TYPE(t_atoms), INTENT(IN) :: atoms
30 : TYPE(t_sym), INTENT(IN) :: sym
31 : TYPE(t_enpara), INTENT(IN) :: enpara
32 : TYPE(t_nococonv), INTENT(IN) :: nococonv
33 : TYPE(t_mpi), INTENT(IN) :: fmpi
34 : TYPE(t_potden), INTENT(IN) :: v, vx
35 : TYPE(t_hub1inp), INTENT(IN) :: hub1inp
36 : TYPE(t_hub1data), INTENT(IN) :: hub1data
37 : TYPE(t_tlmplm), INTENT(INOUT) :: td
38 : TYPE(t_usdus), INTENT(INOUT) :: ud
39 :
40 : ! Indices of atom type, local spins, accessed spin of V and where to start summing lattice harmonics.
41 : INTEGER, INTENT(IN) :: n, ilSpinPr, ilSpin, iSpinV
42 : REAL, INTENT(IN) :: alpha_hybrid
43 :
44 : COMPLEX, INTENT(IN) :: one ! 1 for real part of a component, i if there is an imaginary one.
45 :
46 : REAL, OPTIONAL, INTENT(IN) :: v1(:,:) ! DFPT related. Get t of V1, not V.
47 :
48 : LOGICAL, OPTIONAL, INTENT(IN) :: l_dfpt
49 :
50 2138 : REAL, ALLOCATABLE :: uvu(:,:), uvd(:,:), dvu(:,:), dvd(:,:)
51 2138 : REAL, ALLOCATABLE :: f(:,:,:,:), g(:,:,:,:), flo(:,:,:,:)
52 2138 : REAL, ALLOCATABLE :: vr0(:,:), x(:)
53 :
54 : COMPLEX :: cil
55 : REAL :: temp
56 : INTEGER :: i,l,l2,lamda,lh,lm,lmin,lmin0,lmp,lmx,lp,lh0
57 : INTEGER :: lp1,lpl ,mem,mems,mp,mu,nh,na,m,nsym,s,lplmax
58 : LOGICAL :: l_remove
59 :
60 2138 : lplmax = atoms%lmaxd*(atoms%lmaxd+3)/2
61 :
62 3379306 : ALLOCATE(uvu(0:lplmax,0:sphhar%nlhd)); uvu = 0.0
63 3377168 : ALLOCATE(uvd(0:lplmax,0:sphhar%nlhd)); uvd = 0.0
64 3377168 : ALLOCATE(dvu(0:lplmax,0:sphhar%nlhd)); dvu = 0.0
65 3377168 : ALLOCATE(dvd(0:lplmax,0:sphhar%nlhd)); dvd = 0.0
66 :
67 19242 : ALLOCATE(f(atoms%jmtd,2,0:atoms%lmaxd,2),g(atoms%jmtd,2,0:atoms%lmaxd,2),x(atoms%jmtd))
68 10690 : ALLOCATE(flo(atoms%jmtd,2,atoms%nlod,2))
69 8552 : ALLOCATE( vr0(SIZE(v%mt,1),0:SIZE(v%mt,2)-1))
70 :
71 : !check if we need a l=0 contribution
72 2138 : lh0 = MERGE(1,0,iSpinV<3.and.alpha_hybrid==0)
73 2138 : if (atoms%l_nonpolbas(n)) lh0=0 !non-spin-pol basis
74 2138 : if (PRESENT(v1)) lh0=0 !DFPT code-path
75 :
76 : IF (.NOT.PRESENT(v1)) THEN
77 48719248 : vr0 = v%mt(:,:,n,iSpinV)
78 2138 : IF (iSpinV<3) THEN
79 1400528 : vr0(:,0)=0.0
80 2002 : IF (atoms%l_nonpolbas(n)) THEN
81 60640 : vr0(:,0)=v%mt(:,0,n,iSpinV)-(v%mt(:,0,n,1)+v%mt(:,0,n,2))/2.0
82 60640 : vr0(:,0)= vr0(:,0)/atoms%rmsh(:atoms%jri(n),n)*sfp_const
83 : ENDIF
84 808838 : IF (alpha_hybrid.NE.0) vr0=vr0-alpha_hybrid*vx%mt(:,:,n,iSpinV)
85 : ELSE
86 103008 : vr0(:,0)=vr0(:,0)-0.5*nococonv%b_con(iSpinV-2,n) !Add constraining field
87 : END IF
88 : ELSE
89 0 : vr0=v1
90 : END IF
91 :
92 :
93 4412 : DO i = MIN(ilSpinPr,ilSpin),MAX(ilSpinPr,ilSpin)
94 4412 : CALL genMTBasis(atoms,enpara,v,fmpi,n,i,ud,f(:,:,:,i),g(:,:,:,i),flo(:,:,:,i),hub1data=hub1data)
95 : END DO
96 :
97 2138 : na = atoms%firstAtom(n)
98 2138 : nsym = sym%ntypsy(na)
99 2138 : nh = sphhar%nlh(nsym)
100 :
101 : ! Generate the irreducible integrals
102 : ! <u_{l',order'}^{\mu}|V_{lh}|u_{l,order}^{\mu}>
103 : ! for l <= l' [lower triangle], but only those that will contribute!
104 :
105 21628 : DO lp = 0,atoms%lmax(n)
106 19490 : lp1 = (lp*(lp+1))/2
107 121826 : DO l = 0, lp
108 100198 : lpl = lp1 + l
109 : ! Remove non-spherical components for the orbitals treated with DFT+Hubbard-1
110 100198 : l_remove=.FALSE.
111 100198 : IF(l.EQ.lp.and.atoms%n_u+atoms%n_hia>0) THEN
112 2604 : if (atoms%n_hia>0.and.hub1inp%l_nonsphDC) then
113 0 : DO i = atoms%n_u+1, atoms%n_u+atoms%n_hia
114 0 : IF (atoms%lda_u(i)%atomType.EQ.n.AND.atoms%lda_u(i)%l.EQ.l) l_remove=.TRUE.
115 : END DO
116 : end if
117 2604 : if(atoms%n_u>0.and.input%ldauNonsphDC) then
118 0 : DO i = 1, atoms%n_u
119 0 : IF (atoms%lda_u(i)%atomType.EQ.n.AND.atoms%lda_u(i)%l.EQ.l) l_remove=.TRUE.
120 : END DO
121 : end if
122 : END IF
123 :
124 : ! Loop over the required components of the potential: must satisfy
125 : ! the triangular conditions and that l'+l+l_{lh} even (Gaunt
126 : ! coefficient selection rules)
127 3200658 : DO lh = lh0, nh
128 3080970 : lamda = sphhar%llh(lh,nsym)
129 3080970 : lmin = lp - l
130 3080970 : lmx = lp + l
131 3181168 : IF ((MOD(lamda+lmx,2).EQ.1) .OR. (lamda.LT.lmin) .OR. (lamda.GT.lmx) .OR. l_remove) THEN
132 2033882 : uvu(lpl,lh) = 0.0
133 2033882 : uvd(lpl,lh) = 0.0
134 2033882 : dvu(lpl,lh) = 0.0
135 2033882 : dvd(lpl,lh) = 0.0
136 : ELSE
137 764885396 : DO i = 1,atoms%jri(n)
138 764885396 : x(i) = (f(i,1,lp,ilSpinPr)*f(i,1,l,ilSpin)+f(i,2,lp,ilSpinPr)*f(i,2,l,ilSpin))* vr0(i,lh)
139 : END DO
140 1047088 : CALL intgr3(x,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),temp)
141 1047088 : uvu(lpl,lh) = temp
142 :
143 764885396 : DO i = 1,atoms%jri(n)
144 764885396 : x(i) = (f(i,1,lp,ilSpinPr)*g(i,1,l,ilSpin)+f(i,2,lp,ilSpinPr)*g(i,2,l,ilSpin))* vr0(i,lh)
145 : END DO
146 1047088 : CALL intgr3(x,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),temp)
147 1047088 : uvd(lpl,lh) = temp
148 :
149 764885396 : DO i = 1,atoms%jri(n)
150 764885396 : x(i) = (g(i,1,lp,ilSpinPr)*f(i,1,l,ilSpin)+g(i,2,lp,ilSpinPr)*f(i,2,l,ilSpin))* vr0(i,lh)
151 : END DO
152 1047088 : CALL intgr3(x,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),temp)
153 1047088 : dvu(lpl,lh) = temp
154 :
155 764885396 : DO i = 1,atoms%jri(n)
156 764885396 : x(i) = (g(i,1,lp,ilSpinPr)*g(i,1,l,ilSpin)+g(i,2,lp,ilSpinPr)*g(i,2,l,ilSpin))* vr0(i,lh)
157 : END DO
158 1047088 : CALL intgr3(x,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),temp)
159 1047088 : dvd(lpl,lh) = temp
160 : END IF
161 : END DO
162 : END DO
163 : END DO
164 :
165 : ! Generate the various t matrices for lm <= (lm)'
166 2138 : s = td%h_loc2(n) ! Offset for writing udot elements behind u
167 : ! (lm)' loop:
168 21628 : DO lp = 0,atoms%lmax(n)
169 19490 : lp1 = (lp*(lp+1))/2
170 202534 : DO mp = -lp,lp
171 180906 : lmp = lp* (lp+1) + mp
172 : ! lh loop:
173 5772070 : DO lh = lh0, nh
174 5571674 : lamda = sphhar%llh(lh,nsym)
175 5571674 : lmin0 = ABS(lp-lamda)
176 5571674 : IF (lmin0.GT.lp) CYCLE
177 : ! Ensure l+l'+lamda even
178 4898252 : lmx = lp - MOD(lamda,2)
179 4898252 : mems = sphhar%nmem(lh,nsym)
180 14699004 : DO mem = 1,mems
181 9619846 : mu = sphhar%mlh(mem,lh,nsym)
182 9619846 : m = mp - mu
183 9619846 : lmin = MAX(lmin0,ABS(m))
184 9619846 : l2 = ABS(lmx-lmin)
185 9619846 : lmin = lmin + MOD(l2,2)
186 29980350 : DO l = lmin,lmx,2
187 14788830 : lm = l* (l+1) + m
188 14788830 : IF (lm.GT.lmp) CYCLE
189 12871494 : lpl = lp1 + l
190 : cil = ImagUnit**(l-lp) * sphhar%clnu(mem,lh,nsym) &
191 12871494 : * gaunt1(lp,lamda,l,mp,mu,m,atoms%lmaxd)
192 :
193 12871494 : td%h_loc(lmp,lm,n,ilSpinPr,ilSpin) = td%h_loc(lmp,lm,n,ilSpinPr,ilSpin) + one*cil*uvu(lpl,lh)
194 12871494 : td%h_loc(lmp,lm+s,n,ilSpinPr,ilSpin) = td%h_loc(lmp,lm+s,n,ilSpinPr,ilSpin) + one*cil*uvd(lpl,lh)
195 12871494 : td%h_loc(lmp+s,lm,n,ilSpinPr,ilSpin) = td%h_loc(lmp+s,lm,n,ilSpinPr,ilSpin) + one*cil*dvu(lpl,lh)
196 12871494 : td%h_loc(lmp+s,lm+s,n,ilSpinPr,ilSpin) = td%h_loc(lmp+s,lm+s,n,ilSpinPr,ilSpin) + one*cil*dvd(lpl,lh)
197 : ! Use the fact that the t matrices are Hermitian by definition to construct the upper triangle as well.
198 22491340 : IF (lm.NE.lmp) THEN
199 12186880 : td%h_loc(lm,lmp,n,ilSpinPr,ilSpin) = td%h_loc(lm,lmp,n,ilSpinPr,ilSpin) + one*CONJG(cil*uvu(lpl,lh))
200 12186880 : td%h_loc(lm,lmp+s,n,ilSpinPr,ilSpin) = td%h_loc(lm,lmp+s,n,ilSpinPr,ilSpin) + one*CONJG(cil*dvu(lpl,lh))
201 12186880 : td%h_loc(lm+s,lmp,n,ilSpinPr,ilSpin) = td%h_loc(lm+s,lmp,n,ilSpinPr,ilSpin) + one*CONJG(cil*uvd(lpl,lh))
202 12186880 : td%h_loc(lm+s,lmp+s,n,ilSpinPr,ilSpin) = td%h_loc(lm+s,lmp+s,n,ilSpinPr,ilSpin) + one*CONJG(cil*dvd(lpl,lh))
203 : END IF
204 : END DO
205 : END DO
206 : END DO
207 : END DO
208 : END DO
209 :
210 : ! If necessary, set up the t-matrices for the local orbitals as well.
211 2138 : IF (atoms%nlo(n).GE.1) THEN
212 : CALL tlo(atoms,sym,sphhar,ilSpinPr,ilSpin,iSpinV,n,enpara,lh0,input,vr0,&
213 1166 : na,flo,f,g,ud, td, one, l_dfpt, PRESENT(v1))
214 : END IF
215 2138 : END SUBROUTINE tlmplm
216 : END MODULE m_tlmplm
|
