Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2020 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 : MODULE m_forcea21
7 : CONTAINS
8 62 : SUBROUTINE force_a21(input,atoms,sym ,cell,we,jsp,epar,ne,eig,usdus,tlmplm,&
9 62 : vtot,eigVecCoeffs,aveccof,bveccof,cveccof,f_a21,f_b4,results)
10 : !--------------------------------------------------------------------------
11 : ! Pulay 2nd and 3rd term force contributions à la Rici et al.
12 : !
13 : ! Equation A17 and A20 combined, Phys. Rev. B 43, 6411
14 : !
15 : ! Note1: We do NOT include the i**l factors in the alm, blm coming from
16 : ! to_pulay anymore. Therefore, we can use matrix elements from file 28, 38
17 : ! DIRECTLY.
18 : !
19 : ! Note2: The present version only yields forces for the highest energy window
20 : ! (=valence states). If semicore forces are wanted as well the tmas and tmat
21 : ! files have to be saved, indexed and properly used here in force_a21.
22 : !
23 : ! 22/june/97: Probably found symmetrization error replacing S^-1 by S
24 : ! (IS instead of isinv)
25 : !
26 : ! Force contribution B4 added following
27 : ! Madsen, Blaha, Schwarz, Sjostedt, Nordstrom
28 : ! GMadsen FZJ 20/3-01
29 : !--------------------------------------------------------------------------
30 :
31 : USE m_forcea21lo
32 : USE m_forcea21U
33 : USE m_types_setup
34 : USE m_types_misc
35 : USE m_types_usdus
36 : USE m_types_tlmplm
37 : USE m_types_cdnval
38 : USE m_types_potden
39 : USE m_constants
40 : USE m_juDFT
41 :
42 : IMPLICIT NONE
43 :
44 : TYPE(t_input), INTENT(IN) :: input
45 : TYPE(t_atoms), INTENT(IN) :: atoms
46 : TYPE(t_sym), INTENT(IN) :: sym
47 :
48 : TYPE(t_cell), INTENT(IN) :: cell
49 : TYPE(t_usdus), INTENT(IN) :: usdus
50 : TYPE(t_tlmplm), INTENT(IN) :: tlmplm
51 : TYPE(t_potden), INTENT(IN) :: vtot
52 : TYPE(t_eigVecCoeffs), INTENT(IN) :: eigVecCoeffs
53 : TYPE(t_results), INTENT(INOUT) :: results
54 :
55 : INTEGER, INTENT(IN) :: jsp, ne
56 :
57 : REAL, INTENT(IN) :: we(ne), epar(0:atoms%lmaxd,atoms%ntype)
58 : REAL, INTENT(IN) :: eig(input%neig)
59 : COMPLEX, INTENT(IN) :: aveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
60 : COMPLEX, INTENT(IN) :: bveccof(3,ne,0:atoms%lmaxd*(atoms%lmaxd+2),atoms%nat)
61 : COMPLEX, INTENT(IN) :: cveccof(3,-atoms%llod:atoms%llod,ne,atoms%nlod,atoms%nat)
62 : COMPLEX, INTENT(INOUT) :: f_a21(3,atoms%ntype)
63 : COMPLEX, INTENT(INOUT) :: f_b4(3,atoms%ntype)
64 :
65 : REAL, PARAMETER :: zero=0.0
66 : COMPLEX, PARAMETER :: czero=CMPLX(0.,0.)
67 : COMPLEX dtd, dtu, utd, utu
68 : INTEGER lo
69 : INTEGER i, ie, im, l1, l2, ll1, ll2, lm1, lm2, m1, m2, n, natom, m
70 : INTEGER natrun, is, isinv, j, irinv, it, lmplmd
71 :
72 62 : REAL, ALLOCATABLE :: a21(:,:), b4(:,:)
73 : COMPLEX forc_a21(3), forc_b4(3)
74 : REAL starsum(3), starsum2(3), gvint(3), gvint2(3)
75 : REAL vec(3), vec2(3), vecsum(3), vecsum2(3)
76 :
77 62 : CALL timestart("force_a21")
78 :
79 62 : lmplmd = (atoms%lmaxd*(atoms%lmaxd+2)* (atoms%lmaxd*(atoms%lmaxd+2)+3))/2
80 :
81 248 : ALLOCATE(a21(3,atoms%nat),b4(3,atoms%nat) )
82 :
83 190 : DO n = 1,atoms%ntype
84 128 : natom = atoms%firstAtom(n)
85 190 : IF (atoms%l_geo(n)) THEN
86 128 : forc_a21(:) = czero
87 128 : forc_b4(:) = czero
88 :
89 344 : DO natrun = natom,natom + atoms%neq(n) - 1
90 864 : a21(:,natrun) = zero
91 992 : b4(:,natrun) = zero
92 : END DO
93 :
94 1116 : DO ie = 1,ne
95 9592 : DO l1 = 0,atoms%lmax(n)
96 8476 : ll1 = l1* (l1+1)
97 82836 : DO m1 = -l1,l1
98 73372 : lm1 = ll1 + m1
99 713336 : DO l2 = 0,atoms%lmax(n)
100 639964 : ll2 = l2* (l2+1)
101 6330324 : DO m2 = -l2,l2
102 5616988 : lm2 = ll2 + m2
103 24169640 : DO natrun = natom,natom + atoms%neq(n) - 1
104 17912688 : utu = CONJG(tlmplm%h_loc(lm2,lm1,n,jsp,jsp))
105 17912688 : dtd = CONJG(tlmplm%h_loc(lm2+tlmplm%h_loc2(n),lm1+tlmplm%h_loc2(n),n,jsp,jsp))
106 17912688 : utd = CONJG(tlmplm%h_loc(lm2+tlmplm%h_loc2(n),lm1,n,jsp,jsp))
107 17912688 : dtu = CONJG(tlmplm%h_loc(lm2,lm1+tlmplm%h_loc2(n),n,jsp,jsp))
108 77267740 : DO i = 1,3
109 : a21(i,natrun) = a21(i,natrun) + 2.0*&
110 : AIMAG( CONJG(eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)) *utu*aveccof(i,ie,lm2,natrun)&
111 : +CONJG(eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)) *utd*bveccof(i,ie,lm2,natrun)&
112 : +CONJG(eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)) *dtu*aveccof(i,ie,lm2,natrun)&
113 71650752 : +CONJG(eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)) *dtd*bveccof(i,ie,lm2,natrun))*we(ie)/atoms%neq(n)
114 : END DO ! i (spatial directions)
115 : END DO ! natrun
116 : END DO ! m2
117 : END DO ! l2
118 :
119 : ! Correct spherical part
120 73372 : utu = -eig(ie)
121 73372 : utd = 0.0
122 73372 : dtu = 0.0
123 73372 : dtd = utu*usdus%ddn(l1,n,jsp)
124 301964 : DO i = 1,3
125 975040 : DO natrun = natom,natom + atoms%neq(n) - 1
126 : a21(i,natrun) = a21(i,natrun) + 2.0*AIMAG(&
127 : CONJG(eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp))*utu*aveccof(i,ie,lm1,natrun)&
128 : +CONJG(eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp))*utd*bveccof(i,ie,lm1,natrun)&
129 : +CONJG(eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp))*dtu*aveccof(i,ie,lm1,natrun)&
130 : +CONJG(eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp))*dtd*bveccof(i,ie,lm1,natrun)&
131 901668 : )*we(ie) /atoms%neq(n)
132 : END DO
133 : END DO
134 : END DO ! m1
135 : END DO ! l1
136 : END DO ! ie
137 :
138 : ! Add the local orbital and U contribution to a21:
139 :
140 128 : CALL force_a21_lo(atoms,jsp,n,we,eig,ne,eigVecCoeffs,aveccof,bveccof,cveccof,tlmplm,usdus,a21)
141 :
142 128 : IF (atoms%n_u+atoms%n_hia>0) THEN
143 12 : CALL force_a21_U(atoms,n,jsp,we,ne,usdus,vTot%mmpMat(:,:,:,jsp),eigVecCoeffs,aveccof,bveccof,cveccof,a21)
144 : END IF
145 :
146 128 : IF (input%l_useapw) THEN
147 : ! B4 force
148 0 : DO ie = 1,ne
149 0 : DO l1 = 0,atoms%lmax(n)
150 0 : ll1 = l1* (l1+1)
151 0 : DO m1 = -l1,l1
152 0 : lm1 = ll1 + m1
153 0 : DO i = 1,3
154 0 : DO natrun = natom,natom + atoms%neq(n) - 1
155 : b4(i,natrun) = b4(i,natrun) + 0.5 *&
156 : we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
157 : CONJG(eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)*usdus%us(l1,n,jsp)&
158 : +eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)*usdus%uds(l1,n,jsp))*&
159 : (aveccof(i,ie,lm1,natrun)*usdus%dus(l1,n,jsp)&
160 : +bveccof(i,ie,lm1,natrun)*usdus%duds(l1,n,jsp) )&
161 : -CONJG(aveccof(i,ie,lm1,natrun)*usdus%us(l1,n,jsp)&
162 : +bveccof(i,ie,lm1,natrun)*usdus%uds(l1,n,jsp) )*&
163 : (eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)*usdus%dus(l1,n,jsp)&
164 0 : +eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)*usdus%duds(l1,n,jsp)) )
165 : END DO
166 : END DO
167 : END DO
168 : END DO
169 :
170 0 : DO lo = 1,atoms%nlo(n)
171 0 : l1 = atoms%llo(lo,n)
172 0 : DO m = -l1,l1
173 0 : lm1 = l1* (l1+1) + m
174 0 : DO i=1,3
175 0 : DO natrun = natom,natom + atoms%neq(n) - 1
176 : b4(i,natrun) = b4(i,natrun) + 0.5 *&
177 : we(ie)/atoms%neq(n)*atoms%rmt(n)**2*AIMAG(&
178 : CONJG( eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)* usdus%us(l1,n,jsp)&
179 : + eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)* usdus%uds(l1,n,jsp) ) *&
180 : cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp)&
181 : + CONJG(eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)*usdus%ulos(lo,n,jsp)) *&
182 : ( aveccof(i,ie,lm1,natrun)* usdus%dus(l1,n,jsp)&
183 : + bveccof(i,ie,lm1,natrun)* usdus%duds(l1,n,jsp)&
184 : + cveccof(i,m,ie,lo,natrun)*usdus%dulos(lo,n,jsp) ) &
185 : - (CONJG( aveccof(i,ie,lm1,natrun) *usdus%us(l1,n,jsp)&
186 : + bveccof(i,ie,lm1,natrun) *usdus%uds(l1,n,jsp) ) *&
187 : eigVecCoeffs%ccof(m,ie,lo,natrun,jsp) *usdus%dulos(lo,n,jsp)&
188 : + CONJG(cveccof(i,m,ie,lo,natrun)*usdus%ulos(lo,n,jsp)) *&
189 : ( eigVecCoeffs%abcof(ie,lm1,0,natrun,jsp)*usdus%dus(l1,n,jsp)&
190 : + eigVecCoeffs%abcof(ie,lm1,1,natrun,jsp)*usdus%duds(l1,n,jsp)&
191 0 : + eigVecCoeffs%ccof(m,ie,lo,natrun,jsp)*usdus%dulos(lo,n,jsp) ) ) )
192 : END DO
193 : END DO
194 : END DO
195 : END DO
196 : END DO
197 : END IF
198 :
199 344 : DO natrun = natom,natom + atoms%neq(n) - 1
200 : ! to complete summation over stars of k now sum
201 : ! over all operations which leave (k+G)*R(natrun)*taual(natrun)
202 : ! invariant. We sum over ALL these operations and not only
203 : ! the ones needed for the actual star of k. Should be
204 : ! ok if we divide properly by the number of operations
205 : ! First, we find operation S where RS=T. T -like R- leaves
206 : ! the above scalar product invariant (if S=1 then R=T).
207 : ! R is the operation which generates position of equivalent atom
208 : ! out of position of representative
209 : ! S=R^(-1) T
210 : ! number of ops which leave (k+G)*op*taual invariant: invarind
211 : ! index of inverse operation of R: irinv
212 : ! index of operation T: invarop
213 : ! now, we calculate index of operation S: is
214 : !
215 : ! note, that vector in expression A17,A20 + A21 is a
216 : ! reciprocal lattice vector! other transformation rules
217 : !
218 : ! transform recip vector g-g' into internal coordinates
219 :
220 864 : vec(:) = a21(:,natrun)
221 864 : vec2(:) = b4(:,natrun)
222 :
223 3672 : gvint=MATMUL(cell%bmat,vec)/tpi_const
224 3672 : gvint2=MATMUL(cell%bmat,vec2)/tpi_const
225 216 : vecsum(:) = zero
226 216 : vecsum2(:) = zero
227 :
228 : !-gb2002
229 : ! irinv = invtab(ngopr(natrun))
230 : ! DO it = 1,invarind(natrun)
231 : ! is = multab(irinv,invarop(natrun,it))
232 : !c note, actually we need the inverse of S but -in principle
233 : !c because {S} is a group and we sum over all S- S should also
234 : !c work; to be lucid we take the inverse:
235 : ! isinv = invtab(is)
236 : !! isinv = is
237 : ! Rotation is alreadt done in to_pulay, here we work only in the
238 : ! coordinate system of the representative atom (natom):
239 :
240 752 : DO it = 1,sym%invarind(natom)
241 536 : is =sym%invarop(natom,it)
242 536 : isinv = sym%invtab(is)
243 : !-gb 2002
244 : ! now we have the wanted index of operation with which we have
245 : ! to rotate gv. Note gv is given in cart. coordinates but
246 : ! mrot acts on internal ones
247 2144 : DO i = 1,3
248 1608 : vec(i) = zero
249 1608 : vec2(i) = zero
250 6968 : DO j = 1,3
251 :
252 4824 : vec(i) = vec(i) + sym%mrot(i,j,isinv)*gvint(j)
253 6432 : vec2(i) = vec2(i) + sym%mrot(i,j,isinv)*gvint2(j)
254 :
255 : END DO
256 : END DO
257 2360 : DO i = 1,3
258 1608 : vecsum(i) = vecsum(i) + vec(i)
259 2144 : vecsum2(i) = vecsum2(i) + vec2(i)
260 : END DO
261 : END DO ! end operator loop
262 :
263 : ! Transform from internal to cart. coordinates
264 2808 : starsum=MATMUL(cell%amat,vecsum)
265 2808 : starsum2=MATMUL(cell%amat,vecsum2)
266 992 : DO i = 1,3
267 648 : forc_a21(i) = forc_a21(i) + starsum(i)/sym%invarind(natrun)
268 864 : forc_b4(i) = forc_b4(i) + starsum2(i)/sym%invarind(natrun)
269 : END DO
270 : END DO ! natrun
271 :
272 : ! Add onto existing forces.
273 :
274 : ! NOTE: results%force is real and therefore only the real part of
275 : ! forc_a21 is added. In general, force must be real after the k-star
276 : ! summation. Now, we put the proper operations into real space.
277 : ! Problem: What happens if in real space there is no inversion anymore?
278 : ! But we have inversion in k-space due to time reversal symmetry:
279 : ! E(k)=E(-k)
280 : ! We argue that k-space inversion is automatically taken into account
281 : ! if force = (1/2)(forc_a21+conjg(forc_a21)), because time reversal
282 : ! symmetry means that conjg(PSI) is also a solution of Schrödinger eq.
283 : ! if PSI is one.
284 :
285 512 : DO i = 1, 3
286 384 : results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a21(i) + forc_b4(i))
287 384 : f_a21(i,n) = f_a21(i,n) + forc_a21(i)
288 512 : f_b4(i,n) = f_b4(i,n) + forc_b4(i)
289 : END DO
290 :
291 : END IF ! l_geo(n)
292 : END DO
293 :
294 : ! The result is written in force_a8.
295 :
296 62 : CALL timestop("force_a21")
297 :
298 62 : END SUBROUTINE force_a21
299 : END MODULE m_forcea21
|
