Branch data 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 : : !!***
|
