Line data Source code
1 : !!****m* m_paw_yukawa/m_paw_yukawa
2 : !! NAME
3 : !! m_paw_yukawa
4 : !!
5 : !! FUNCTION
6 : !! This module contains several routines related to the Yukawa parametrization
7 : !! of Coulomb interactions in the PAW approach.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2025-2026 ABINIT group
11 : !! These routines are inspired by K. Haule routines in embedded DMFT.
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public License, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !!
16 : !! SOURCE
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : MODULE m_paw_yukawa
25 :
26 : use defs_basis
27 : use m_abicore
28 : use m_errors
29 : use m_pawrad, only : pawrad_type
30 :
31 : implicit none
32 :
33 : private
34 :
35 : public :: compute_slater
36 : public :: get_lambda
37 :
38 : CONTAINS !========================================================================================
39 : !!***
40 :
41 : !!****f* m_paw_yukawa/compute_slater
42 : !! NAME
43 : !! compute_slater
44 : !!
45 : !! FUNCTION
46 : !!
47 : !! Compute Slater integrals for a screened Yukawa potential v(r,r') = exp(-lambda*(r-r'))/(epsilon*(r-r'))
48 : !! This is eq. 36 in the supplementary of Physical review letters, Haule, K. (2015), 115(19), 196403
49 : !!
50 : !! INPUTS
51 : !! lpawu = angular momentum
52 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data:
53 : !! %mesh_size=Dimension of radial mesh
54 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
55 : !! proj2 = u(r)**2 where u(r) is the atomic orbital multiplied by r
56 : !! meshsz = size of the radial mesh
57 : !! lambda, eps = parameters of the Yukawa potential
58 : !!
59 : !! OUTPUT
60 : !! fk(lpawu+1)= Slater integrals
61 : !!
62 : !! SOURCE
63 :
64 0 : subroutine compute_slater(lpawu,pawrad,proj2,meshsz,lambda,eps,fk)
65 :
66 : use m_pawrad, only : pawrad_type,simp_gen
67 : use m_bessel2, only : bessel_iv,bessel_kv
68 :
69 : !Arguments ------------------------------------
70 : integer, intent(in) :: lpawu,meshsz
71 : real(dp), intent(in) :: lambda,eps
72 : real(dp), intent(in) :: proj2(meshsz)
73 : real(dp), intent(inout) :: fk(lpawu+1)
74 : type(pawrad_type), intent(in) :: pawrad
75 : !Local variables ------------------------------
76 : integer :: ir,k,mesh_type
77 : real(dp) :: dum,r_for_intg,y0
78 0 : real(dp), allocatable :: r_k(:),u_inside(:),u_outside(:),y1(:),y2(:)
79 : !************************************************************************
80 :
81 0 : mesh_type = pawrad%mesh_type
82 0 : r_for_intg = pawrad%rad(meshsz)
83 0 : ABI_MALLOC(r_k,(meshsz))
84 0 : ABI_MALLOC(u_inside,(meshsz))
85 0 : ABI_MALLOC(u_outside,(meshsz))
86 :
87 0 : if (lambda == zero) then
88 :
89 0 : do k=0,2*lpawu+1,2
90 0 : u_inside(1) = zero
91 0 : r_k(:) = pawrad%rad(1:meshsz)**k
92 0 : do ir=2,meshsz
93 0 : if (ir == 2) then
94 : ! Use a trapezoidal rule
95 : u_inside(ir) = half * (proj2(1)*r_k(1)+proj2(2)*r_k(2)) * &
96 0 : & (pawrad%rad(2)-pawrad%rad(1))
97 0 : else if (ir == 3 .and. mesh_type == 3) then
98 : ! simp_gen doesn't handle this case, so we use a trapezoidal rule instead
99 : u_inside(ir) = u_inside(2) + half*(proj2(2)*r_k(2)+proj2(3)*r_k(3))* &
100 0 : & (pawrad%rad(3)-pawrad%rad(2))
101 : else
102 : ! Use Simpson rule when enough points are available
103 0 : call simp_gen(u_inside(ir),proj2(1:ir)*r_k(1:ir),pawrad,r_for_intg=pawrad%rad(ir))
104 : end if ! ir=2
105 : end do ! ir
106 0 : u_outside(1) = zero
107 0 : u_outside(2:meshsz) = two * u_inside(2:meshsz) * proj2(2:meshsz) / (pawrad%rad(2:meshsz)*r_k(2:meshsz))
108 0 : call simp_gen(fk(k/2+1),u_outside(:),pawrad,r_for_intg=r_for_intg)
109 : end do ! k
110 :
111 : else
112 :
113 0 : ABI_MALLOC(y1,(meshsz))
114 0 : ABI_MALLOC(y2,(meshsz))
115 :
116 0 : r_k(:) = sqrt(pawrad%rad(1:meshsz))
117 :
118 0 : do k=0,2*lpawu+1,2
119 :
120 0 : y0 = zero
121 0 : if (k == 0) y0 = lambda * sqrt(two/pi)
122 0 : y1(1) = y0
123 0 : u_inside(1) = zero
124 :
125 0 : do ir=2,meshsz
126 :
127 0 : call bessel_iv(half+dble(k),lambda*pawrad%rad(ir),zero,y1(ir),dum)
128 0 : y1(ir) = y1(ir) / r_k(ir)
129 :
130 0 : if (ir == 2) then
131 : ! Use a trapezoidal rule
132 : u_inside(ir) = half * (proj2(1)*y1(1)+proj2(2)*y1(2)) * &
133 0 : & (pawrad%rad(2)-pawrad%rad(1))
134 0 : else if (ir == 3 .and. mesh_type == 3) then
135 : ! simp_gen doesn't handle this case, so we use a trapezoidal rule instead
136 : u_inside(ir) = u_inside(2) + half*(proj2(2)*y1(2)+proj2(3)*y1(3)) * &
137 0 : & (pawrad%rad(3)-pawrad%rad(2))
138 : else
139 : ! Use Simpson rule when enough points are available
140 0 : call simp_gen(u_inside(ir),proj2(1:ir)*y1(1:ir),pawrad,r_for_intg=pawrad%rad(ir))
141 : end if ! ir
142 :
143 0 : call bessel_kv(half+dble(k),lambda*pawrad%rad(ir),zero,y2(ir),dum)
144 0 : y2(ir) = y2(ir) / r_k(ir)
145 : end do ! ir
146 :
147 0 : u_outside(1) = zero
148 0 : u_outside(2:meshsz) = two * (two*dble(k)+one)*u_inside(2:meshsz)*proj2(2:meshsz)*y2(2:meshsz)
149 0 : call simp_gen(fk(k/2+1),u_outside(:),pawrad,r_for_intg=r_for_intg)
150 :
151 : end do ! k
152 :
153 0 : ABI_FREE(y1)
154 0 : ABI_FREE(y2)
155 :
156 : end if ! lambda
157 :
158 0 : fk(1:lpawu+1) = fk(1:lpawu+1) / eps
159 :
160 0 : ABI_FREE(r_k)
161 0 : ABI_FREE(u_inside)
162 0 : ABI_FREE(u_outside)
163 :
164 0 : end subroutine compute_slater
165 : !!***
166 :
167 : !----------------------------------------------------------------------
168 :
169 : !!****f* m_paw_yukawa/get_lambda
170 : !! NAME
171 : !! get_lambda
172 : !!
173 : !! FUNCTION
174 : !!
175 : !! Conversion from U,J,f4/f2,f6/f2 parametrization of Slater integrals
176 : !! to lambda,epsilon parametrization, with lambda and epsilon the parameters
177 : !! of the screened Yukawa potential v(r,r') = exp(-lambda*(r-r'))/(epsilon*(r-r')).
178 : !!
179 : !! CAREFUL: this routine does not handle custom f4/f2 and f6/f2 values, and set them
180 : !! to their default values f4/f2=0.625 for l=2 and f4/f2=0.6681, f6/f2=0.4943 for l=3.
181 : !!
182 : !! CAREFUL: For l>=2 we lose information since we have more input parameters than output
183 : !! parameters. In this case, a compromise has to be made, and you will no longer
184 : !! have the exact same Slater integrals as before.
185 : !!
186 : !! INPUTS
187 : !! lpawu = angular momentum
188 : !! pawrad <type(pawrad_type)>=paw radial mesh and related data:
189 : !! %mesh_size=Dimension of radial mesh
190 : !! %rad(mesh_size)=The coordinates of all the points of the radial mesh
191 : !! proj2 = u(r)**2 where u(r) is the atomic orbital multiplied by r
192 : !! meshsz = size of the radial mesh
193 : !! upawu,jpawu = parameters for Slater integrals
194 : !! yukawa_param = if set to 1, search for lambda and epsilon yielding the values closest to u and j
195 : !! if set to 2, search for lambda yielding u, and set epsilon to 1
196 : !!
197 : !! OUTPUT
198 : !! lambda,epsilon = parameters of the corresponding Yukawa potential
199 : !!
200 : !! SOURCE
201 :
202 0 : subroutine get_lambda(lpawu,pawrad,proj2,meshsz,upawu,jpawu,lambda,eps,yukawa_param)
203 :
204 : use m_brentq, only : brentq
205 : use m_hybrd, only : hybrd
206 :
207 : !Arguments ------------------------------------
208 : integer, intent(in) :: lpawu,meshsz,yukawa_param
209 : real(dp), intent(in) :: upawu,jpawu
210 : real(dp), intent(out) :: lambda,eps
211 : real(dp), intent(in) :: proj2(meshsz)
212 : type(pawrad_type), intent(in) :: pawrad
213 : !Local variables ------------------------------
214 : integer :: i,ierr,info,ldfjac,lr,maxfev,ml,mode,mu,n,nfev,nprint
215 : real(dp) :: epsfcn,fac,lmb_temp,upbound,xtol
216 0 : real(dp) :: diag(2),fjac(2,2),fkk(lpawu+1),fvec(2),lmb_eps(2)
217 : real(dp) :: r(3),qtf(2),wa1(2),wa2(2),wa3(2),wa4(2)
218 : character(len=500) :: message
219 : !************************************************************************
220 :
221 : ! Find suitable upper bound for brentq routine
222 0 : upbound = five
223 0 : do i=1,10
224 :
225 0 : call compute_slater(lpawu,pawrad,proj2(:),meshsz,upbound,one,fkk(:))
226 0 : if (fkk(1) < upawu) exit
227 0 : upbound = two * upbound
228 :
229 : end do ! i
230 :
231 0 : write(message,'(4a)') "An error occurred when trying to find a suitable lambda and ", &
232 0 : & "epsilon for your input values of upawu and jpawu.", ch10, &
233 0 : & "Either try different values or use dmft_yukawa_lambda and dmft_yukawa_epsilon."
234 :
235 0 : if (fkk(1) > upawu) ABI_ERROR(message)
236 :
237 : ! First, set epsilon to 1, and find lambda which yields the correct F0=upawu, to have a good starting point
238 0 : call brentq(get_coulomb_u,zero,upbound,two*tol12,four*epsilon(one),100,lmb_temp,ierr)
239 :
240 0 : if (ierr == 0) ABI_ERROR(message)
241 :
242 : ! Initial values for lambda and epsilon
243 0 : lmb_eps(1) = lmb_temp
244 0 : lmb_eps(2) = one
245 :
246 0 : lambda = lmb_temp
247 0 : eps = one
248 :
249 0 : if (yukawa_param == 2) return
250 :
251 0 : if (lpawu > 0) then
252 :
253 : ! Default values from scipy
254 0 : epsfcn = epsilon(one) ; fac = dble(100.) ; n = 2
255 0 : ldfjac = n ; lr = n * (n+1) / 2
256 0 : maxfev = 200 * (n+1) ; ml = n - 1 ; mode = 1
257 0 : mu = n - 1 ; nprint = 0 ; xtol = dble(1.49012e-8)
258 :
259 : ! Now find lambda and epsilon
260 : call hybrd(get_coulomb_uj,2,lmb_eps(:),fvec(:),xtol,maxfev,ml,mu,epsfcn,diag(:),mode, &
261 0 : & fac,nprint,info,nfev,fjac(:,:),ldfjac,r(:),lr,qtf(:),wa1(:),wa2(:),wa3(:),wa4(:))
262 :
263 0 : if (info /= 1) ABI_ERROR(message)
264 :
265 : end if ! lpawu > 0
266 :
267 0 : lambda = lmb_eps(1)
268 0 : eps = lmb_eps(2)
269 :
270 : contains
271 :
272 0 : subroutine get_coulomb_u(lmb,uu)
273 :
274 : !Arguments ------------------------------------
275 : real(dp), intent(in) :: lmb
276 : real(dp), intent(out) :: uu
277 : !Local variables ------------------------------
278 0 : real(dp) :: fk(lpawu+1)
279 : !************************************************************************
280 :
281 0 : call compute_slater(lpawu,pawrad,proj2(:),meshsz,lmb,one,fk(:))
282 0 : uu = fk(1) - upawu
283 :
284 0 : end subroutine get_coulomb_u
285 :
286 0 : subroutine get_coulomb_uj(n,lmb_eps,uj,iflag)
287 :
288 : !Arguments ------------------------------------
289 : integer, intent(in) :: iflag,n
290 : real(dp), intent(in) :: lmb_eps(n)
291 : real(dp), intent(inout) :: uj(n)
292 : !Local variables ------------------------------
293 : real(dp) :: eps,f4of2,f6of2,factor,j2,j4,j6,jh,lmb
294 0 : real(dp) :: fk(lpawu+1)
295 : character(len=500) :: message
296 : !************************************************************************
297 :
298 : ABI_UNUSED(iflag)
299 :
300 0 : lmb = lmb_eps(1)
301 0 : eps = lmb_eps(2)
302 0 : call compute_slater(lpawu,pawrad,proj2(:),meshsz,lmb,eps,fk(:))
303 0 : uj(1) = fk(1) - upawu
304 :
305 0 : if (lpawu == 1) then
306 0 : j2 = fk(2) * fifth
307 0 : jh = j2
308 0 : else if (lpawu == 2) then
309 0 : f4of2 = dble(0.625)
310 0 : factor = (one+f4of2) / dble(14)
311 0 : j2 = fk(2)
312 0 : j4 = fk(3) / f4of2
313 0 : jh = (j2+j4) * factor * half
314 0 : else if (lpawu == 3) then
315 0 : f4of2 = dble(0.6681)
316 0 : f6of2 = dble(0.4943)
317 0 : factor = (dble(286.)+dble(195.)*f4of2+dble(250.)*f6of2) / dble(6435.)
318 0 : j2 = fk(2)
319 0 : j4 = fk(3) / f4of2
320 0 : j6 = fk(4) / f6of2
321 0 : jh = (j2+j4+j6) * factor * third
322 : else
323 0 : write(message,'(a,i0,2a)') ' lpawu=',lpawu,ch10,' lpawu not equal to 0, 1, 2 or 3 is not allowed'
324 0 : ABI_ERROR(message)
325 : end if ! lpawu
326 :
327 0 : uj(2) = jh - jpawu
328 :
329 0 : end subroutine get_coulomb_uj
330 :
331 : end subroutine get_lambda
332 : !!***
333 :
334 : !----------------------------------------------------------------------
335 :
336 : END MODULE m_paw_yukawa
337 : !!***
|
