Line data Source code
1 : !-------------------------------------------------------------------------------
2 : ! Copyright (c) 2022 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !-------------------------------------------------------------------------------
6 :
7 : MODULE m_dfpt_rhomt
8 : !! Module adapted from rhomt21, which in turn was adapted from rhomt, and
9 : !! from rhomtlo for the lo part.
10 : !!
11 : !! It contains subroutines for the non-LO and LO spherical contributions
12 : !! to the density coefficients used to construct \(\rho_{0}^{\alpha}(r)\).
13 : CONTAINS
14 8382 : SUBROUTINE dfpt_rhomt(atoms,we,we1,ne,ilSpinPr,ilSpin,qpoint,l_dfpt,eigVecCoeffs,eigVecCoeffs1,denCoeffs,eigVecCoeffs1m)
15 : !! Subroutine to construct all spherical MT density coefficients (for a
16 : !! density perturbation) without LOs in one routine. The spin input dictates,
17 : !! which element is gonna be built.
18 : !! The coefficients are of the form:
19 : !! \begin{aligned}
20 : !! d_{l,\lambda',\lambda}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha} &= \sum_{\nu\boldsymbol{k}}\sum_{m}\\
21 : !! &* A_{l,m,\lambda'}^{\sigma_{\alpha}',\nu\boldsymbol{k}*}\\
22 : !! &* (2\tilde{f}_{\nu\boldsymbol{k}}A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}\boldsymbol{q},j,\beta~(1)}+\tilde{f}_{\nu\boldsymbol{k}\boldsymbol{q}}^{j,\beta~(1)}A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}})
23 : !! \end{aligned}
24 : !! The k-point loop is performed outside this routine. In contrast to older
25 : !! routines, the arrays uu etc. and uu21 etc. are merged into one
26 : !! spin and u-order dependent array.
27 : !!
28 : !! \(\sigma_{\alpha}(')\): local spin indices \(\rightarrow\) ilSpinPr, ilSpin
29 : !!
30 : !! \(\lambda(')\): order of the radial basis function (0: u, 1: d)
31 : !!
32 : !! \(\nu\boldsymbol{k}\): State index (k-point and number of state)
33 : !!
34 : !! \(\boldsymbol{q},j,\beta\): Perturbation index; q-point, direction and atom
35 : !!
36 : !! \(\tilde{f}_{\nu\boldsymbol{k}}\): (Smeared) occupation number [perturbed for \(\tilde{f}^{(1)}\)]
37 : !!
38 : !! \(A\): Summed matching coefficients and eigenvectors [perturbed for \(A^{(1)}\)]
39 :
40 : USE m_types_atoms
41 : USE m_types_cdnval
42 :
43 : IMPLICIT NONE
44 :
45 : TYPE(t_atoms), INTENT(IN) :: atoms
46 :
47 : INTEGER, INTENT(IN) :: ne
48 : INTEGER, INTENT(IN) :: ilSpinPr !! \(\sigma_{\alpha}^{'}\)
49 : INTEGER, INTENT(IN) :: ilSpin !! \(\sigma_{\alpha}\)
50 : REAL, INTENT(IN) :: we(ne) !! \(\tilde{f}_{\nu\boldsymbol{k}}\)
51 : REAL, INTENT(IN) :: we1(ne) !! \(\tilde{f}_{\nu\boldsymbol{k}\boldsymbol{q}}^{j,\beta~(1)}\)
52 : REAL, INTENT(IN) :: qpoint(3) !! \(\boldsymbol{q}\)
53 : LOGICAL, INTENT(IN) :: l_dfpt
54 :
55 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs !! \(A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}}\)
56 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs1 !! \(A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}\boldsymbol{q},j,\beta~(1)}\)
57 :
58 : TYPE(t_denCoeffs), INTENT(INOUT) :: denCoeffs !! \(d_{l,\lambda',\lambda}^{\sigma_{\alpha}',\sigma_{\alpha},\alpha}\)
59 :
60 : TYPE(t_eigVecCoeffs), OPTIONAL, INTENT(IN) :: eigVecCoeffs1m !! \(A_{l,m,\lambda}^{\sigma_{\alpha},\nu\boldsymbol{k}\boldsymbol{-q},j,\beta~(1)}\)
61 :
62 : INTEGER i,l,lm,itype,na,natom,lo,lop,m
63 :
64 : COMPLEX :: temp
65 :
66 : LOGICAL :: l_minusq
67 :
68 8382 : l_minusq = PRESENT(eigVecCoeffs1m)
69 8382 : natom = 0
70 26088 : DO itype = 1,atoms%ntype
71 44088 : DO na = 1,atoms%neq(itype)
72 18000 : natom = natom + 1
73 199194 : DO l = 0,atoms%lmax(itype)
74 1706256 : DO m = -l,l
75 1524768 : lm = l* (l+1) + m
76 22907167 : DO i = 1,ne
77 : ! uu/du
78 21218911 : temp = we(i) * eigVecCoeffs1%abcof(i,lm,0,natom,ilSpin) ! If not DFPT, this is the base case for rhomt(21)
79 21218911 : IF (l_dfpt) THEN
80 0 : IF (.NOT.l_minusq) temp = temp * 2.0
81 0 : IF (norm2(qpoint)<=1e-8) THEN
82 0 : temp = temp + we1(i) * eigVecCoeffs%abcof(i,lm,0,natom,ilSpin)
83 : END IF
84 : END IF
85 : denCoeffs%mt_coeff(l,itype,0,0,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,0,0,ilSpinPr, ilSpin) &
86 21218911 : & + CONJG(eigVecCoeffs%abcof(i,lm,0,natom,ilSpinPr)) * temp
87 : denCoeffs%mt_coeff(l,itype,1,0,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,1,0,ilSpinPr, ilSpin) &
88 21218911 : & + CONJG(eigVecCoeffs%abcof(i,lm,1,natom,ilSpinPr)) * temp
89 21218911 : IF (l_minusq) THEN
90 : denCoeffs%mt_coeff(l,itype,0,0,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,0,0,ilSpinPr, ilSpin) &
91 : & + CONJG(eigVecCoeffs1m%abcof(i,lm,0,natom,ilSpinPr)) &
92 0 : & * eigVecCoeffs%abcof(i,lm,0,natom,ilSpin)
93 : denCoeffs%mt_coeff(l,itype,1,0,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,1,0,ilSpinPr, ilSpin) &
94 : & + CONJG(eigVecCoeffs1m%abcof(i,lm,1,natom,ilSpinPr)) &
95 0 : & * eigVecCoeffs%abcof(i,lm,0,natom,ilSpin)
96 : END IF
97 : ! ud/dd
98 21218911 : temp = we(i) * eigVecCoeffs1%abcof(i,lm,1,natom,ilSpin)
99 21218911 : IF (l_dfpt) THEN
100 0 : IF (.NOT.l_minusq) temp = temp * 2.0
101 0 : IF (norm2(qpoint)<=1e-8) THEN
102 0 : temp = temp + we1(i) * eigVecCoeffs%abcof(i,lm,1,natom,ilSpin)
103 : END IF
104 : END IF
105 : denCoeffs%mt_coeff(l,itype,0,1,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,0,1,ilSpinPr, ilSpin) &
106 21218911 : & + CONJG(eigVecCoeffs%abcof(i,lm,0,natom,ilSpinPr)) * temp
107 : denCoeffs%mt_coeff(l,itype,1,1,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,1,1,ilSpinPr, ilSpin) &
108 21218911 : & + CONJG(eigVecCoeffs%abcof(i,lm,1,natom,ilSpinPr)) * temp
109 22743679 : IF (l_minusq) THEN
110 : denCoeffs%mt_coeff(l,itype,0,1,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,0,1,ilSpinPr, ilSpin) &
111 : & + CONJG(eigVecCoeffs1m%abcof(i,lm,0,natom,ilSpinPr)) &
112 0 : & * eigVecCoeffs%abcof(i,lm,1,natom,ilSpin)
113 : denCoeffs%mt_coeff(l,itype,1,1,ilSpinPr, ilSpin) = denCoeffs%mt_coeff(l,itype,1,1,ilSpinPr, ilSpin) &
114 : & + CONJG(eigVecCoeffs1m%abcof(i,lm,1,natom,ilSpinPr)) &
115 0 : & * eigVecCoeffs%abcof(i,lm,1,natom,ilSpin)
116 : END IF
117 : END DO
118 : END DO
119 : END DO
120 : END DO
121 : END DO
122 8382 : END SUBROUTINE dfpt_rhomt
123 :
124 8382 : SUBROUTINE dfpt_rhomtlo(atoms,ne,we,we1,ilSpinPr,ilSpin,qpoint,l_dfpt,eigVecCoeffs,eigVecCoeffs1,denCoeffs,eigVecCoeffs1m)
125 : !! This is a complementary routine to the one above for \(\lambda(')\ge 2\),
126 : !! i.e. mixed or pure LO contributions.
127 :
128 : USE m_types_atoms
129 : USE m_types_cdnval
130 :
131 : IMPLICIT NONE
132 :
133 : TYPE(t_atoms), INTENT(IN) :: atoms
134 :
135 : INTEGER, INTENT(IN) :: ne, ilSpinPr, ilSpin
136 :
137 : REAL, INTENT(IN) :: we(:),we1(:)!(nobd)
138 : REAL, INTENT(IN) :: qpoint(3)
139 : LOGICAL, INTENT(IN) :: l_dfpt
140 :
141 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs, eigVecCoeffs1
142 : TYPE(t_denCoeffs), INTENT(INOUT) :: denCoeffs
143 :
144 : TYPE(t_eigVecCoeffs), OPTIONAL, INTENT(IN) :: eigVecCoeffs1m
145 :
146 : INTEGER i,l,lm,lo,lop ,natom,nn,ntyp,m
147 :
148 : COMPLEX :: temp
149 :
150 : LOGICAL :: l_minusq
151 :
152 8382 : l_minusq = PRESENT(eigVecCoeffs1m)
153 :
154 8382 : natom = 0
155 26088 : DO ntyp = 1,atoms%ntype
156 44088 : DO nn = 1,atoms%neq(ntyp)
157 18000 : natom = natom + 1
158 54186 : DO lo = 1,atoms%nlo(ntyp)
159 18480 : l = atoms%llo(lo,ntyp)
160 61700 : DO m = -l,l
161 43220 : lm = l* (l+1) + m
162 733703 : DO i = 1,ne
163 672003 : temp = we(i) * eigVecCoeffs1%ccof(m,i,lo,natom,ilSpin) ! If not DFPT, this is the base case for rhomt(21)
164 672003 : IF (l_dfpt) THEN
165 0 : temp = temp * 2.0
166 0 : IF (norm2(qpoint)<=1e-8) THEN
167 0 : temp = temp + we1(i) * eigVecCoeffs%ccof(m,i,lo,natom,ilSpin)
168 : END IF
169 : END IF
170 : denCoeffs%mt_ulo_coeff(lo,ntyp,0,ilSpinPr,ilSpin) = denCoeffs%mt_ulo_coeff(lo,ntyp,0,ilSpinPr,ilSpin) &
171 672003 : & + CONJG(eigVecCoeffs%abcof(i,lm,0,natom,ilSpinPr)) * temp
172 : denCoeffs%mt_ulo_coeff(lo,ntyp,1,ilSpinPr,ilSpin) = denCoeffs%mt_ulo_coeff(lo,ntyp,1,ilSpinPr,ilSpin) &
173 672003 : & + CONJG(eigVecCoeffs%abcof(i,lm,1,natom,ilSpinPr)) * temp
174 672003 : temp = we(i) * eigVecCoeffs1%abcof(i,lm,0,natom,ilSpin)
175 672003 : IF (l_dfpt) THEN
176 0 : temp = temp * 2.0
177 0 : IF (norm2(qpoint)<=1e-8) THEN
178 0 : temp = temp + we1(i) * eigVecCoeffs%abcof(i,lm,0,natom,ilSpin)
179 : END IF
180 : END IF
181 : denCoeffs%mt_lou_coeff(lo,ntyp,0,ilSpinPr,ilSpin) = denCoeffs%mt_lou_coeff(lo,ntyp,0,ilSpinPr,ilSpin) &
182 672003 : & + CONJG(eigVecCoeffs%ccof(m,i,lo,natom,ilSpinPr)) * temp
183 672003 : temp = we(i) * eigVecCoeffs1%abcof(i,lm,1,natom,ilSpin)
184 672003 : IF (l_dfpt) THEN
185 0 : temp = temp * 2.0
186 0 : IF (norm2(qpoint)<=1e-8) THEN
187 0 : temp = temp + we1(i) * eigVecCoeffs%abcof(i,lm,1,natom,ilSpin)
188 : END IF
189 : END IF
190 : denCoeffs%mt_lou_coeff(lo,ntyp,1,ilSpinPr,ilSpin) = denCoeffs%mt_lou_coeff(lo,ntyp,1,ilSpinPr,ilSpin) &
191 715223 : & + CONJG(eigVecCoeffs%ccof(m,i,lo,natom,ilSpinPr)) * temp
192 : END DO
193 : END DO
194 68460 : DO lop = 1,atoms%nlo(ntyp)
195 50460 : IF (atoms%llo(lop,ntyp).EQ.l) THEN
196 61860 : DO m = -l,l
197 739755 : DO i = 1,ne
198 677895 : temp = we(i) * eigVecCoeffs1%ccof(m,i,lo,natom,ilSpin) ! If not DFPT, this is the base case for rhomt(21)
199 677895 : IF (l_dfpt) THEN
200 0 : temp = temp * 2.0
201 0 : IF (norm2(qpoint)<=1e-8) THEN
202 0 : temp = temp + we1(i) * eigVecCoeffs%ccof(m,i,lo,natom,ilSpin)
203 : END IF
204 : END IF
205 : denCoeffs%mt_lolo_coeff(lop,lo,ntyp,ilSpinPr,ilSpin) = denCoeffs%mt_lolo_coeff(lop,lo,ntyp,ilSpinPr,ilSpin) &
206 721235 : & + CONJG(eigVecCoeffs%ccof(m,i,lop,natom,ilSpinPr)) * temp
207 : END DO
208 : END DO
209 : END IF
210 : END DO
211 : END DO
212 : END DO
213 : END DO
214 8382 : END SUBROUTINE dfpt_rhomtlo
215 : END MODULE m_dfpt_rhomt
|
