Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 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_force_sf
8 : USE m_constants
9 : USE m_types
10 : !-----------------------------------------------------------------------------
11 : ! This routine calculates a contribution to the forces stemming
12 : ! from the discontinuity of density and potential at the muffin-tin
13 : ! boundary oint n [ rho V (IS) - rho V (MT) ] dS
14 : ! Klueppelberg May 13
15 : !
16 : ! The surface integral is found in the first two lines of equation (14) in
17 : !
18 : ! Klüppelberg et al., PRB 91 035105 (2015).
19 : !
20 : ! Klueppelberg (force level 3)
21 : !-----------------------------------------------------------------------------
22 :
23 : IMPLICIT NONE
24 :
25 : COMPLEX, PRIVATE, SAVE, ALLOCATABLE :: force_mt(:,:)
26 : COMPLEX, PRIVATE, SAVE, ALLOCATABLE :: force_is(:,:)
27 : LOGICAL, PRIVATE, SAVE :: isdone=.false.,mtdone=.false.
28 :
29 : CONTAINS
30 1 : SUBROUTINE force_sf_is(atoms_in,stars,sym,jsp,cell,qpw,vpw,excpw,vxcpw)
31 : !--------------------------------------------------------------------------
32 : ! This subroutine calculates oint n rho V dS evaluated with quantities
33 : ! from the interstital. The Fourier transform of density and potential is
34 : ! first calculated as an expansion in spherical harmonics. Then the normal
35 : ! vector, which is proportional to Y_1t, connects the l component of rho
36 : ! with the l+-1 components of V. This is done up to a cutoff lmax.
37 : ! It is called in a spin loop at the end of vgen.F90-
38 : !--------------------------------------------------------------------------
39 :
40 : USE m_sphbes
41 : USE m_phasy1
42 : USE m_gaunt
43 :
44 : TYPE(t_sym),INTENT(IN) :: sym
45 : TYPE(t_stars),INTENT(IN) :: stars
46 : TYPE(t_cell),INTENT(IN) :: cell
47 : TYPE(t_atoms),INTENT(IN) :: atoms_in
48 :
49 : INTEGER, INTENT (IN) :: jsp
50 :
51 : COMPLEX, INTENT (IN) :: qpw(:,:) !(stars%ng3,input%jspins)
52 : COMPLEX, INTENT (IN) :: vpw(:,:) !(n3d,jspd)
53 : COMPLEX, INTENT (IN) :: excpw(stars%ng3)
54 : COMPLEX, INTENT (IN) :: vxcpw(:,:) !(stars%ng3,input%jspins)
55 :
56 1 : TYPE(t_atoms) :: atoms
57 :
58 : INTEGER :: n,j,itype,s,l ,lm,t,lp,mp,lmp,jp,natom,m
59 : REAL :: r,r2
60 : COMPLEX :: img,rhoprep,Vprep
61 : LOGICAL :: isthere
62 :
63 1 : INTEGER :: lmaxb(atoms_in%ntype)
64 1 : COMPLEX :: coeff(3,-1:1),qpw2(stars%ng3,size(qpw,2)),qpwcalc(stars%ng3,size(qpw,2))
65 1 : REAL , ALLOCATABLE :: bsl(:,:,:)
66 1 : COMPLEX, ALLOCATABLE :: pylm(:,:,:),rho(:),V(:),pylm2(:,:)
67 :
68 1 : CALL timestart("Force level 3 (IS)")
69 :
70 4 : atoms=atoms_in
71 5 : atoms%lmax = 2*atoms_in%lmax
72 1 : atoms%lmaxd = 2*atoms_in%lmaxd
73 4 : lmaxb = atoms%lmax
74 1 : img = cmplx(0.0,1.0)
75 :
76 1 : CALL init_sf(sym,cell,atoms)
77 1 : isdone = .true.
78 :
79 5 : ALLOCATE ( bsl(0:atoms%lmaxd,stars%ng3,atoms%ntype) )
80 :
81 4 : ALLOCATE ( pylm2((atoms%lmaxd+1)**2,atoms%ntype ))
82 4 : ALLOCATE ( rho((atoms%lmaxd+1)**2),V((atoms%lmaxd+1)**2) )
83 :
84 1 : coeff(:, :) = cmplx(0.0,0.0)
85 1 : coeff(1,-1) = sqrt(tpi_const/3.)
86 1 : coeff(1, 1) = -sqrt(tpi_const/3.)
87 1 : coeff(2,-1) = img*sqrt(tpi_const/3.)
88 1 : coeff(2, 1) = img*sqrt(tpi_const/3.)
89 1 : coeff(3, 0) = sqrt(2.*tpi_const/3.)
90 :
91 3556 : qpwcalc = qpw
92 :
93 4 : DO itype = 1,atoms%ntype
94 3 : CALL sphbes(atoms%lmax(itype),0.0,bsl(:,1,itype))
95 10663 : DO s = 2,stars%ng3
96 : ! Only call sphbes if the length of the star changed.
97 10662 : IF (abs(stars%sk3(s)-stars%sk3(s-1)).gt.1.0e-14) THEN
98 10659 : r = stars%sk3(s)*atoms%rmt(itype)
99 10659 : CALL sphbes(atoms%lmax(itype),r,bsl(:,s,itype))
100 : ELSE
101 0 : bsl(:,s,itype) = bsl(:,s-1,itype)
102 : END IF
103 : END DO ! s
104 : END DO ! itype
105 :
106 :
107 13 : force_is = 0.0
108 :
109 4 : DO itype = 1,atoms%ntype
110 3 : natom = atoms%firstAtom(itype)
111 3 : r2 = atoms%rmt(itype)**2
112 870 : rho = 0.0
113 870 : V = 0.0
114 :
115 10665 : DO s = 1,stars%ng3
116 10662 : CALL phasy1(atoms,stars,sym,cell,s,pylm2(:,:))
117 :
118 181257 : DO l = 0,atoms%lmax(itype)-1
119 170592 : rhoprep = stars%nstr(s) * bsl(l,s,itype) * qpwcalc(s,jsp)
120 170592 : IF (l.eq.0) THEN
121 : Vprep = stars%nstr(s) * bsl(l,s,itype) * &
122 10662 : (vpw(s,jsp)-vxcpw(s,jsp)+excpw(s))
123 : ! Switching between Veff and VCoul + exc
124 : END IF
125 : ! For l = 0 we calculate rho_00 and V_00:
126 2900064 : DO m = -l,l
127 2729472 : lm = l*(l+1) + m + 1
128 2729472 : rho(lm) = rho(lm) + rhoprep * pylm2(lm,itype)!pylm(lm,s,itype)!
129 2729472 : IF (l.gt.0) CYCLE
130 2900064 : V(lm) = V(lm) + Vprep * pylm2(lm,itype)!pylm(lm,s,itype)!
131 : END DO ! m
132 : ! And V_1mp, for l > 0, we calculate rho_lm, V_l+1,mp:
133 : Vprep = stars%nstr(s) * bsl(l+1,s,itype) * &
134 170592 : (vpw(s,jsp)-vxcpw(s,jsp)+excpw(s))
135 3251910 : DO m = -l-1,l+1
136 3070656 : lm = (l+1)*(l+2) + m + 1
137 3241248 : V(lm) = V(lm) + Vprep * pylm2(lm,itype)!pylm(lm,s,itype)!
138 : END DO ! m
139 : END DO ! l
140 : END DO ! s
141 :
142 52 : DO l = 0,atoms%lmax(itype)-1 ! new: altered s and l loop above
143 819 : DO m = -l,l
144 768 : lm = l*(l+1) + m + 1
145 : ! Because rho_lm occurs with V_l-1,mp and V_l+1,mp:
146 2349 : DO lp = abs(l-1),l+1,2
147 6900 : DO t = -1,1
148 4599 : mp = t-m
149 4599 : IF (lp.lt.abs(mp)) CYCLE
150 4329 : lmp = lp*(lp+1) + mp + 1
151 : force_is(:,itype) = force_is(:,itype) + r2 * rho(lm) * &
152 : V(lmp) * conjg(coeff(:,t)) * &
153 18849 : gaunt2(1,l,lp,t,m,mp,atoms%lmax(itype))
154 : END DO ! t
155 : END DO ! lp
156 : END DO ! m
157 : END DO ! l
158 : END DO ! itype
159 :
160 1 : DEALLOCATE ( bsl,rho,V )
161 1 : DEALLOCATE ( pylm2 )
162 :
163 1 : CALL timestop("Force level 3 (IS)")
164 :
165 4 : END SUBROUTINE force_sf_is
166 :
167 1 : SUBROUTINE force_sf_mt(atoms,sphhar,jspin,ispin,fmpi,vr,excr,vxcr,rho,sym,cell )
168 : !--------------------------------------------------------------------------
169 : ! This subroutine calculates the contribution evaluated with quantities
170 : ! from the Muffin Tins.
171 : !
172 : ! n rho V = sum(nu,nup) rho(nu)V(nup)
173 : ! * sum(m,mu,mup) c_1m* Y_1m* c_lnu Y_numu c_lnup Y_nupmup
174 : !
175 : ! It is called in a spin loop at the end of cdngen.F90
176 : !--------------------------------------------------------------------------
177 :
178 : USE m_gaunt
179 : USE m_ylm
180 :
181 : TYPE(t_mpi),INTENT(IN) :: fmpi
182 : TYPE(t_sym),INTENT(IN) :: sym
183 : TYPE(t_cell),INTENT(IN) :: cell
184 : TYPE(t_sphhar),INTENT(IN):: sphhar
185 : TYPE(t_atoms),INTENT(IN) :: atoms
186 :
187 : INTEGER, INTENT (IN) :: jspin
188 : INTEGER, INTENT (IN) :: ispin
189 :
190 : REAL , INTENT (IN) :: vr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype) !
191 : REAL , INTENT (IN) :: rho(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
192 : REAL , INTENT (IN) :: excr(atoms%jmtd,0:sphhar%nlhd,atoms%ntype)
193 : REAL , INTENT (IN) :: vxcr(:,0:,:,:)!(atoms%jmtd,0:sphhar%nlhd,atoms%ntype,input%jspins)
194 :
195 : INTEGER :: natom,itype,nd,lh,l,lhp,lp,mem,m,memp,mp,t,i,lmp
196 : REAL :: pot,den
197 : COMPLEX :: img,factor
198 :
199 : COMPLEX :: coeff(3,-1:1)
200 1 : COMPLEX :: d1((atoms%lmaxd+1)**2,atoms%ntype ),d2((atoms%lmaxd+1)**2,atoms%ntype )
201 :
202 1 : CALL timestart("Force level 3 (MT)")
203 :
204 1 : CALL init_sf(sym,cell,atoms)
205 1 : mtdone = .true.
206 :
207 1 : img = cmplx(0.0,1.0)
208 13 : force_mt = 0.0
209 :
210 1 : coeff(:, :) = cmplx(0.0,0.0)
211 1 : coeff(1,-1) = sqrt(tpi_const/3.)
212 1 : coeff(1, 1) = -sqrt(tpi_const/3.)
213 1 : coeff(2,-1) = img*sqrt(tpi_const/3.)
214 1 : coeff(2, 1) = img*sqrt(tpi_const/3.)
215 1 : coeff(3, 0) = sqrt(2.*tpi_const/3.)
216 :
217 247 : d1 = 0
218 247 : d2 = 0
219 :
220 : ! Calculate forces: For each atom, loop over all lattice harmonics.
221 4 : DO itype = 1,atoms%ntype
222 3 : natom = atoms%firstAtom(itype)
223 3 : nd = sym%ntypsy(natom)
224 :
225 247 : DO lh = 0,sphhar%nlh(nd)
226 243 : l = sphhar%llh(lh,nd)
227 :
228 : ! The l=0 component of the potential array is saved with an additional
229 : ! factor r/sfp. For this calculation, we need the pure potential.
230 243 : pot = vr(atoms%jri(itype),lh,itype)
231 243 : IF (lh.eq.0) THEN
232 3 : pot = pot*sfp_const/atoms%rmt(itype)
233 : END IF
234 :
235 243 : pot = excr(atoms%jri(itype),lh,itype)-vxcr(atoms%jri(itype),lh,itype,ispin)+pot
236 :
237 702 : DO mem = 1,sphhar%nmem(lh,nd)
238 459 : m = sphhar%mlh(mem,lh,nd)
239 459 : lmp = l*(l+1) + m + 1
240 : d1(lmp,itype) = d1(lmp,itype) + sphhar%clnu(mem,lh,nd) * &
241 459 : rho(atoms%jri(itype),lh,itype,ispin)/atoms%rmt(itype)**2
242 702 : d2(lmp,itype) = d2(lmp,itype) + sphhar%clnu(mem,lh,nd) * pot
243 : END DO ! mem
244 :
245 19929 : DO lhp = 0,sphhar%nlh(nd)
246 19683 : lp = sphhar%llh(lhp,nd)
247 19683 : IF (abs(l-lp).ne.1) CYCLE
248 :
249 4848 : den = rho(atoms%jri(itype),lhp,itype,ispin)
250 :
251 14355 : DO mem = 1,sphhar%nmem(lh,nd)
252 9264 : m = sphhar%mlh(mem,lh,nd)
253 46659 : DO memp = 1,sphhar%nmem(lhp,nd)
254 17712 : mp = sphhar%mlh(memp,lhp,nd)
255 17712 : IF (abs(m+mp).gt.1) CYCLE
256 :
257 : ! Due to the normal vector n, the lattice harmonics form a
258 : ! Gaunt coefficient with the Y_1m.
259 : factor = pot * den * sphhar%clnu(mem,lh,nd) * sphhar%clnu(memp,lhp,nd)&
260 4104 : * gaunt1(1,l,lp,m+mp,m,mp,atoms%lmaxd)
261 :
262 25680 : force_mt(:,itype) = force_mt(:,itype) + factor * conjg(coeff(:,m+mp))
263 :
264 : END DO ! memp
265 : END DO ! mem
266 :
267 : END DO ! lhp
268 :
269 : END DO ! lh
270 : END DO ! itype
271 :
272 1 : CALL timestop("Force level 3 (MT)")
273 :
274 1 : END SUBROUTINE force_sf_mt
275 :
276 4 : SUBROUTINE init_sf(sym,cell,atoms)
277 : ! Initialize force arrays if neither force_sf_is nor force_sf_mt were
278 : ! executed yet.
279 : ! Called at the beginning of cdnval.F90 to once fill the wigner array
280 : ! Also called in force_sf_is/mt to guarantee that all arrays are allocated.
281 :
282 : TYPE(t_sym),INTENT(IN) :: sym
283 : TYPE(t_cell),INTENT(IN) :: cell
284 : TYPE(t_atoms),INTENT(IN) :: atoms
285 :
286 4 : IF (isdone.OR.mtdone) RETURN
287 2 : IF (.NOT.ALLOCATED(force_is)) THEN
288 8 : ALLOCATE ( force_is(3,atoms%ntype),force_mt(3,atoms%ntype) )
289 : END IF
290 26 : force_is = 0.0
291 26 : force_mt = 0.0
292 :
293 : END SUBROUTINE init_sf
294 :
295 2 : SUBROUTINE exit_sf(isp,atoms,force)
296 : ! Write out the force contribution from the surface terms and deallocate
297 : ! arrays if all force_sf_is and force_sf_mt were executed.
298 : ! Called at the end of totale.f90.
299 :
300 : INTEGER,INTENT(IN) :: isp
301 : TYPE(t_atoms),INTENT(IN) :: atoms
302 : REAL,INTENT(INOUT) :: force(:,:,:)
303 :
304 : INTEGER :: itype,dir
305 2 : COMPLEX :: force_sf(3,atoms%ntype)
306 :
307 2 : IF (isdone.AND.mtdone) THEN
308 13 : force_sf(:,:) = force_is(:,:) - force_mt(:,:)
309 13 : force(:,:,isp) = force(:,:,isp) + real(force_sf(:,:))
310 1 : WRITE (oUnit,*)
311 4 : DO itype = 1,atoms%ntype
312 3 : WRITE (oUnit,FMT=8010) itype
313 4 : WRITE (oUnit,FMT=8020) (force_sf(dir,itype),dir=1,3)
314 : END DO ! itype
315 1 : isdone = .false.
316 1 : mtdone = .false.
317 1 : DEALLOCATE ( force_is,force_mt )
318 : END IF
319 :
320 : 8010 FORMAT (' FORCES: SURFACE CORRECTION FOR ATOM TYPE',i4)
321 : 8020 FORMAT (' FX_SF=',2f10.6,' FY_SF=',2f10.6,' FZ_SF=',2f10.6)
322 :
323 2 : END SUBROUTINE exit_sf
324 :
325 : END MODULE m_force_sf
|