Line data Source code
1 : MODULE m_forcea8
2 : CONTAINS
3 29 : SUBROUTINE force_a8(input,atoms,sym,sphhar,jsp,vr,rho,force,fmpi,results)
4 : !--------------------------------------------------------------------------
5 : ! Pulay 1st term force contribution à la Rici et al.
6 : !
7 : ! Equation A8, Phys. Rev. B 43, 6411
8 : !--------------------------------------------------------------------------
9 : USE m_intgr, ONLY : intgr3
10 : USE m_constants
11 : USE m_gaunt, ONLY :gaunt1
12 : USE m_differentiate,ONLY: difcub
13 : USE m_types
14 : USE m_juDFT
15 :
16 : IMPLICIT NONE
17 :
18 : TYPE(t_input), INTENT(IN) :: input
19 : TYPE(t_atoms), INTENT(IN) :: atoms
20 : TYPE(t_sym), INTENT(IN) :: sym
21 : TYPE(t_sphhar), INTENT(IN) :: sphhar
22 : TYPE(t_force), INTENT(IN) :: force
23 : TYPE(t_mpi), INTENT(IN) :: fmpi
24 : TYPE(t_results), INTENT(INOUT) :: results
25 :
26 : INTEGER, INTENT(IN) :: jsp
27 :
28 : REAL, INTENT(IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
29 : REAL, INTENT(IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
30 :
31 : ! Local scalars
32 : COMPLEX, PARAMETER :: czero=CMPLX(0.0,0.0)
33 : COMPLEX aaa, bbb, ccc, ddd, eee, fff
34 : REAL a8_1, a8_2, qval, rr, truerho, truev, xi
35 : INTEGER i, ir, j, l, l1, l2, lh1, lh2, m1, m2, mem1, mem2, n, nd, na, m
36 :
37 : ! Local arrays
38 : COMPLEX forc_a8(3), gv(3), f_sum(3)
39 29 : REAL rhoaux(atoms%jmtd), rhodif(atoms%jmtd)
40 :
41 : ! Statement functions
42 : REAL alpha, beta, delta, epslon, gamma, phi
43 : INTEGER krondel
44 :
45 : ! Kronecker delta for arguments >=0 AND <0
46 : krondel(i,j) = MIN(ABS(i)+1,ABS(j)+1)/MAX(ABS(i)+1,ABS(j)+1)* (1+SIGN(1,i)*SIGN(1,j))/2
47 : alpha(l,m) = (l+1)*0.5e0*SQRT(REAL((l-m)* (l-m-1))/ REAL((2*l-1)* (2*l+1)))
48 : beta(l,m) = l*0.5e0*SQRT(REAL((l+m+2)* (l+m+1))/ REAL((2*l+1)* (2*l+3)))
49 : GAMMA(l,m) = (l+1)*0.5e0*SQRT(REAL((l+m)* (l+m-1))/ REAL((2*l-1)* (2*l+1)))
50 : delta(l,m) = l*0.5e0*SQRT(REAL((l-m+2)* (l-m+1))/ REAL((2*l+1)* (2*l+3)))
51 : epslon(l,m) = (l+1)*SQRT(REAL((l-m)* (l+m))/ REAL((2*l-1)* (2*l+1)))
52 : phi(l,m) = l*SQRT(REAL((l-m+1)* (l+m+1))/REAL((2*l+1)* (2*l+3)))
53 :
54 29 : CALL timestart("force_a8")
55 :
56 29 : WRITE (oUnit,*)
57 :
58 89 : DO n = 1, atoms%ntype
59 60 : na = atoms%firstAtom(n)
60 89 : IF (atoms%l_geo(n)) THEN
61 60 : nd = sym%ntypsy(na)
62 :
63 240 : DO i = 1,3
64 240 : forc_a8(i) = czero
65 : END DO
66 :
67 : ! TODO: There is no output for this. Do we want some?
68 :
69 60 : CALL intgr3(rho(:,0,n,jsp),atoms%rmsh(:,n),atoms%dx(n),atoms%jri(n),qval)
70 :
71 : ! Check if the l=0 density is correct.
72 : ! Note that in general also all l>0 components of the density have
73 : ! been multiplied by r**2 (see for example checkdop which constructs
74 : ! the true MT-density at rmt).
75 : ! The factor sqrt(4*pi) comes from Y_00*\int d\Omega = 1/sqrt(4*pi)*4*pi
76 : 8000 FORMAT (' FORCE_A8: valence charge=',1p,e16.8)
77 :
78 : ! PART I of FORCE_A8
79 1726 : DO lh1 = 0, sphhar%nlh(nd)
80 1666 : l1 = sphhar%llh(lh1,nd)
81 68780 : DO lh2 = 0, sphhar%nlh(nd)
82 67054 : l2 = sphhar%llh(lh2,nd)
83 :
84 268216 : DO i = 1, 3
85 268216 : gv(i) = czero
86 : END DO
87 :
88 : ! Sum over all m for a particular lattice harmonic.
89 194876 : DO mem1 = 1, sphhar%nmem(lh1,nd)
90 127822 : m1 = sphhar%mlh(mem1,lh1,nd)
91 450218 : DO mem2 = 1, sphhar%nmem(lh2,nd)
92 255342 : m2 = sphhar%mlh(mem2,lh2,nd)
93 : gv(1) = gv(1) + SQRT(2.e0*pi_const/3.e0) * &
94 : sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
95 : (gaunt1(1,l1,l2,-1,m1,m2,atoms%lmaxd) - &
96 255342 : gaunt1(1,l1,l2,1,m1,m2,atoms%lmaxd))
97 : gv(2) = gv(2) - ImagUnit*SQRT(2.e0*pi_const/3.e0) * &
98 : sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
99 : (gaunt1(1,l1,l2,-1,m1,m2,atoms%lmaxd) + &
100 255342 : gaunt1(1,l1,l2,1,m1,m2,atoms%lmaxd))
101 : gv(3) = gv(3) + SQRT(4.e0*pi_const/3.e0) * &
102 : sphhar%clnu(mem1,lh1,nd)*sphhar%clnu(mem2,lh2,nd)*&
103 383164 : gaunt1(1,l1,l2,0,m1,m2,atoms%lmaxd)
104 : END DO
105 : END DO
106 :
107 : ! Note that in general also all l>0 components of the density
108 : ! have been multiplied by r**2.
109 : ! Here we need the true radial denisity for performing the
110 : ! derivative. Rherefore we divide by r**2.
111 29326232 : DO ir = 1, atoms%jri(n)
112 29326232 : rhoaux(ir) = rho(ir,lh2,n,jsp)/ (atoms%rmsh(ir,n)**2)
113 : END DO
114 :
115 : ! NOTE: Here we should have: vr = vtrue
116 : ! difcub performs the analytic derivative of Lagrangian of 3rd order
117 67054 : xi = atoms%rmsh(1,n)
118 67054 : rr = xi*xi
119 67054 : rhodif(1) = difcub(atoms%rmsh(1,n),rhoaux(1),xi)*vr(1,lh1,n)*rr
120 29125070 : DO ir = 2, atoms%jri(n) - 2
121 29058016 : xi = atoms%rmsh(ir,n)
122 29058016 : rr = xi*xi
123 29125070 : rhodif(ir) = difcub(atoms%rmsh(ir-1,n),rhoaux(ir-1),xi)*vr(ir,lh1,n)*rr
124 : END DO
125 :
126 67054 : xi = atoms%rmsh(atoms%jri(n)-1,n)
127 67054 : rr = xi*xi
128 : rhodif(atoms%jri(n)-1) = difcub(atoms%rmsh(atoms%jri(n)-3,n), &
129 : rhoaux(atoms%jri(n)-3),xi) * &
130 67054 : vr(atoms%jri(n)-1,lh1,n)*rr
131 :
132 67054 : xi = atoms%rmsh(atoms%jri(n),n)
133 67054 : rr = xi*xi
134 : rhodif(atoms%jri(n)) = difcub(atoms%rmsh(atoms%jri(n)-3,n), &
135 : rhoaux(atoms%jri(n)-3),xi) * &
136 67054 : vr(atoms%jri(n),lh1,n)*rr
137 :
138 : ! NOTE: vr(l=0) is EXPLICITELY multiplied by r/sqrt(4pi) to be
139 : ! the TRUE r*V which is needed for the radial Schrödinger eq.
140 : ! Here, we need the l=0 component, v_{l=0}(r) which will be
141 : ! multiplied by Y_00 in the lm expansion; therefore we MUST
142 : ! recorrect vr(l=0) by the inverse factor sqrt(4pi)/r. We do
143 : ! the correction for the product array rhodif.
144 67054 : IF (lh1.EQ.0) THEN
145 683288 : DO ir = 1,atoms%jri(n)
146 683288 : rhodif(ir) = rhodif(ir)/atoms%rmsh(ir,n)*sfp_const
147 : END DO
148 : END IF
149 :
150 67054 : CALL intgr3(rhodif,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),a8_1)
151 :
152 269882 : DO i = 1,3
153 268216 : forc_a8(i) = forc_a8(i) + a8_1*gv(i)
154 : END DO
155 :
156 : END DO !lh2 (0:sphhar%nlh(nd))
157 : END DO ! lh1 (0:sphhar%nlh(nd))
158 :
159 : ! PART II of FORCE_A8
160 1726 : DO lh1 = 0, sphhar%nlh(nd)
161 1666 : l1 = sphhar%llh(lh1,nd)
162 68780 : DO lh2 = 0, sphhar%nlh(nd)
163 67054 : l2 = sphhar%llh(lh2,nd)
164 29326232 : DO ir = 1, atoms%jri(n)
165 29259178 : truev = vr(ir,lh1,n)
166 29259178 : truerho = rho(ir,lh2,n,jsp)/ (atoms%rmsh(ir,n)**2)
167 29326232 : rhoaux(ir) = truev*truerho*atoms%rmsh(ir,n)
168 : END DO
169 :
170 67054 : IF (lh1.EQ.0) THEN
171 683288 : DO ir = 1,atoms%jri(n)
172 683288 : rhoaux(ir) = rhoaux(ir)/atoms%rmsh(ir,n)*sfp_const
173 : END DO
174 : END IF
175 :
176 67054 : CALL intgr3(rhoaux,atoms%rmsh(1,n),atoms%dx(n),atoms%jri(n),a8_2)
177 :
178 268216 : DO i = 1, 3
179 268216 : gv(i) = (0.0,0.0)
180 : END DO
181 :
182 : ! Sum over all m for a particular lattice harmonic.
183 194876 : DO mem1 = 1, sphhar%nmem(lh1,nd)
184 127822 : m1 = sphhar%mlh(mem1,lh1,nd)
185 450218 : DO mem2 = 1, sphhar%nmem(lh2,nd)
186 255342 : m2 = sphhar%mlh(mem2,lh2,nd)
187 :
188 : ! NOTE: delta(-m,m')(-1)^m was applied because we have the integrand
189 : ! Y nabla Y instead of Y* nabla Y.
190 : ! We also know that Y*(lm) = (-1)^m Y(l,-m)
191 : ! and Y(lm) = (-1)^m Y*(l,-m).
192 : ! Therefore (-1)^m delta(-m,m') appears.
193 :
194 : aaa = alpha(l2,m2)*sphhar%clnu(mem1,lh1,nd) * &
195 : sphhar%clnu(mem2,lh2,nd) * &
196 : (-1)**m1*krondel(l1,l2-1)* &
197 255342 : krondel(-m1,m2+1)
198 : bbb = beta(l2,m2)*sphhar%clnu(mem1,lh1,nd) * &
199 : sphhar%clnu(mem2,lh2,nd) * &
200 : (-1)**m1*krondel(l1,l2+1)* &
201 255342 : krondel(-m1,m2+1)
202 : ccc = GAMMA(l2,m2)*sphhar%clnu(mem1,lh1,nd) * &
203 : sphhar%clnu(mem2,lh2,nd) * &
204 : (-1)**m1*krondel(l1,l2-1)* &
205 255342 : krondel(-m1,m2-1)
206 : ddd = delta(l2,m2)*sphhar%clnu(mem1,lh1,nd) * &
207 : sphhar%clnu(mem2,lh2,nd) * &
208 : (-1)**m1*krondel(l1,l2+1)* &
209 255342 : krondel(-m1,m2-1)
210 : eee = epslon(l2,m2)*sphhar%clnu(mem1,lh1,nd)* &
211 : sphhar%clnu(mem2,lh2,nd)* &
212 : (-1)**m1*krondel(l1,l2-1)*&
213 255342 : krondel(-m1,m2)
214 : fff = phi(l2,m2)*sphhar%clnu(mem1,lh1,nd) * &
215 : sphhar%clnu(mem2,lh2,nd) * &
216 : (-1)**m1*krondel(l1,l2+1)* &
217 255342 : krondel(-m1,m2)
218 :
219 255342 : gv(1) = gv(1) + aaa + bbb - ccc - ddd
220 255342 : gv(2) = gv(2) - ImagUnit* (aaa+bbb+ccc+ddd)
221 383164 : gv(3) = gv(3) + eee - fff
222 : END DO
223 : END DO
224 :
225 336936 : DO i = 1, 3
226 268216 : forc_a8(i) = forc_a8(i) + a8_2*gv(i)
227 : END DO
228 :
229 : END DO ! lh2 (0:sphhar%nlh(nd))
230 : END DO ! lh1 (0:sphhar%nlh(nd))
231 :
232 : ! Add onto existing forces.
233 :
234 240 : DO i = 1, 3
235 240 : results%force(i,n,jsp) = results%force(i,n,jsp) + REAL(forc_a8(i))
236 : END DO
237 :
238 : ! Write out result.
239 60 : WRITE (oUnit,FMT=8010) n
240 60 : WRITE (oUnit,FMT=8020) (forc_a8(i),i=1,3)
241 : 8010 FORMAT (' FORCES: EQUATION A8 FOR ATOM TYPE',i4)
242 : 8020 FORMAT (' FX_A8=',2f10.6,' FY_A8=',2f10.6,' FZ_A8=',2f10.6)
243 :
244 : END IF
245 : END DO
246 :
247 : ! Write out the result of a12, a21, b4 and b8
248 : ! here as well.
249 :
250 29 : IF (.NOT.input%l_useapw) THEN
251 :
252 29 : WRITE (oUnit,*)
253 :
254 29 : IF (fmpi%isize.EQ.1) THEN
255 0 : DO n=1, atoms%ntype
256 0 : IF (atoms%l_geo(n)) THEN
257 0 : WRITE (oUnit,FMT=8030) n
258 0 : WRITE (oUnit,FMT=8040) (force%f_a12(i,n),i=1,3)
259 : END IF
260 : 8030 FORMAT (' FORCES: EQUATION A12 FOR ATOM TYPE',i4)
261 : 8040 FORMAT (' FX_A12=',2f10.6,' FY_A12=',2f10.6,' FZ_A12=',2f10.6)
262 : END DO
263 : ELSE
264 29 : WRITE (oUnit,*) "If this was a serial calculation, the A12 force component would be written out here. In parallel it holds no meaning."
265 : END IF
266 : ELSE
267 :
268 0 : WRITE (oUnit,*)
269 :
270 0 : DO n=1, atoms%ntype
271 0 : IF (atoms%l_geo(n)) THEN
272 0 : WRITE (oUnit,FMT=8070) n
273 0 : WRITE (oUnit,FMT=8080) (force%f_b4(i,n),i=1,3)
274 : END IF
275 : 8070 FORMAT (' FORCES: EQUATION B4 FOR ATOM TYPE',i4)
276 : 8080 FORMAT (' FX_B4=',2f10.6,' FY_B4=',2f10.6,' FZ_B4=',2f10.6)
277 : END DO
278 :
279 0 : WRITE (oUnit,*)
280 :
281 0 : DO n=1,atoms%ntype
282 0 : IF (atoms%l_geo(n)) THEN
283 0 : WRITE (oUnit,FMT=8090) n
284 0 : WRITE (oUnit,FMT=8100) (force%f_b8(i,n),i=1,3)
285 : END IF
286 : 8090 FORMAT (' FORCES: EQUATION B8 FOR ATOM TYPE',i4)
287 : 8100 FORMAT (' FX_B8=',2f10.6,' FY_B8=',2f10.6,' FZ_B8=',2f10.6)
288 : END DO
289 : END IF
290 :
291 29 : WRITE (oUnit,*)
292 :
293 29 : IF (fmpi%isize.EQ.1) THEN
294 0 : DO n=1,atoms%ntype
295 0 : IF (atoms%l_geo(n)) THEN
296 0 : WRITE (oUnit,FMT=8050) n
297 0 : WRITE (oUnit,FMT=8060) (force%f_a21(i,n),i=1,3)
298 : END IF
299 : 8050 FORMAT (' FORCES: EQUATION A21 FOR ATOM TYPE',i4)
300 : 8060 FORMAT (' FX_A21=',2f10.6,' FY_A21=',2f10.6,' FZ_A21=',2f10.6)
301 : END DO
302 : ELSE
303 29 : WRITE (oUnit,*) "If this was a serial calculation, the A21 force component would be written out here. In parallel it holds no meaning."
304 : END IF
305 :
306 :
307 29 : CALL timestop("force_a8")
308 :
309 29 : END SUBROUTINE force_a8
310 : END MODULE m_forcea8
|
