Branch data Line data Source code
1 : : !!****m* ABINIT/m_pawrad
2 : : !! NAME
3 : : !! m_pawrad
4 : : !!
5 : : !! FUNCTION
6 : : !! Module containing all the functions related to the PAW radial meshes
7 : : !!
8 : : !! COPYRIGHT
9 : : !! Copyright (C) 2013-2026 ABINIT group (MT,FJ,MG)
10 : : !! This file is distributed under the terms of the
11 : : !! GNU General Public License, see ~abinit/COPYING
12 : : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : : !!
14 : : !! NOTES
15 : : !! * Routines tagged with "@type_name" are strongly connected to the definition of the data type.
16 : : !! Strongly connected means that the proper functioning of the implementation relies on the
17 : : !! assumption that the tagged procedure is consistent with the type declaration.
18 : : !! Every time a developer changes the structure "type_name" adding new entries, he/she has to make sure
19 : : !! that all the strongly connected routines are changed accordingly to accommodate the modification of the data type
20 : : !! Typical examples of strongly connected routines are creation, destruction or reset methods.
21 : : !!
22 : : !! * FOR DEVELOPERS: in order to preserve the portability of libPAW library,
23 : : !! please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt
24 : : !!
25 : : !! SOURCE
26 : :
27 : : #include "libpaw.h"
28 : :
29 : : MODULE m_pawrad
30 : :
31 : : USE_DEFS
32 : : USE_MSG_HANDLING
33 : : USE_MPI_WRAPPERS
34 : : USE_MEMORY_PROFILING
35 : :
36 : : use m_paw_numeric, only : paw_derfc
37 : :
38 : : implicit none
39 : :
40 : : private
41 : :
42 : : !public procedures.
43 : : public :: pawrad_init ! Main creation method
44 : : public :: pawrad_free ! Free the allocated memory
45 : : public :: pawrad_print ! Printout of the basic info
46 : : public :: pawrad_isame ! Checks whether two meshes are equivalent or have the same equation.
47 : : public :: pawrad_copy ! Returns a copy of the mesh.
48 : : public :: pawrad_ifromr ! Retrieve the Index FROM a given R value in a radial grid.
49 : : public :: pawrad_deducer0 ! Extrapolate r=0 value of a function from values near r=0.
50 : : public :: pawrad_bcast ! Broadcast pawrad datastructure over a given MPI communicator
51 : : public :: simp_gen ! Performs integral on a given (generalized) grid using Simpson rule.
52 : : public :: nderiv_gen ! Do corrected first (and 2nd) derivation on a given (generalized) grid.
53 : : public :: nderiv_lin ! Do corrected first (and 2nd) derivation on a given linear grid.
54 : : public :: bound_deriv ! Computes derivatives at boundaries of the mesh
55 : : public :: poisson ! Solves Poisson eq. for angularly dependent charge distribution of angular momentum l
56 : : public :: screened_coul_kernel ! Kernel used to compute short-range screened Coulomb integrals
57 : : public :: calc_slatradl ! Calculates the radial part of Slater integrals.
58 : :
59 : : interface pawrad_free
60 : : module procedure pawrad_free_0D
61 : : module procedure pawrad_free_1D
62 : : end interface pawrad_free
63 : :
64 : : ! TODO: Might use bit flags, but all radmesh stuff should be encapsulated here!
65 : : integer,private,parameter :: RMESH_LINEAR = 1
66 : : integer,private,parameter :: RMESH_LOG1 = 2
67 : : integer,private,parameter :: RMESH_LOG2 = 3
68 : : integer,private,parameter :: RMESH_LOG3 = 4
69 : : integer,private,parameter :: RMESH_NL = 5
70 : : !!***
71 : :
72 : : !-------------------------------------------------------------------------
73 : :
74 : : !!****t* m_pawrad/pawrad_type
75 : : !! NAME
76 : : !! pawrad_type
77 : : !!
78 : : !! FUNCTION
79 : : !! For PAW, RADial mesh discretization and related data
80 : : !!
81 : : !! SOURCE
82 : :
83 : : type, public :: pawrad_type
84 : :
85 : : !Integer scalars
86 : :
87 : : integer :: int_meshsz=0
88 : : ! Mesh size used in integrals computation
89 : : ! Integrals will be computed up to r(int_meshsz)
90 : :
91 : : integer :: mesh_size=0
92 : : ! Dimension of radial mesh
93 : :
94 : : integer :: mesh_type=-1
95 : : ! Type of mesh
96 : : ! 1=regular grid: r(i)=(i-1)*AA
97 : : ! 2=logarithmic grid: r(i)=AA*(exp[BB*(i-1)]-1)
98 : : ! 3=logarithmic grid: r(i>1)=AA*exp[BB*(i-1)] and r(1)=0
99 : : ! 4=logarithmic grid: r(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
100 : :
101 : : !Real (real(dp)) scalars
102 : :
103 : : real(dp) :: lstep=zero
104 : : ! Exponential step of the mesh (BB parameter above)
105 : : ! Defined only if mesh type is logarithmic
106 : :
107 : : real(dp) :: rmax=zero
108 : : ! Max. value of r = rad(mesh_size)
109 : :
110 : : real(dp) :: rstep=zero
111 : : ! Radial step of the mesh (AA parameter above)
112 : :
113 : : real(dp) :: stepint=zero
114 : : ! Radial step used to convert any function from the
115 : : ! present grid onto a regular grid in order to
116 : : ! integrate it using trapeze method
117 : :
118 : : !Real (real(dp)) arrays
119 : :
120 : : real(dp), allocatable :: rad(:)
121 : : ! rad(mesh_size)
122 : : ! Coordinates of all the points of the mesh
123 : :
124 : : real(dp), allocatable :: radfact(:)
125 : : ! radfact(mesh_size)
126 : : ! Factor used to compute radial integrals
127 : : ! Before being integrated on the present mesh,
128 : : ! any function is multiplied by this factor
129 : :
130 : : real(dp), allocatable :: simfact(:)
131 : : ! simfact(mesh_size)
132 : : ! Factor used to compute radial integrals by the a Simpson scheme
133 : : ! Integral[f] = Sum_i [simfact(i)*f(i)]
134 : :
135 : : end type pawrad_type
136 : : !!***
137 : :
138 : : CONTAINS
139 : : !===========================================================
140 : : !!***
141 : :
142 : : !!****f* m_pawrad/pawrad_init
143 : : !! NAME
144 : : !! pawrad_init
145 : : !!
146 : : !! FUNCTION
147 : : !! Creation method for radial meshes.
148 : : !! Compute all points (and related weights) of a radial mesh.
149 : : !! Grid can be regular or logarithimc.
150 : : !!
151 : : !! INPUTS
152 : : !! [mesh_size]=Dimension of the radial mesh
153 : : !! If not present, take mesh%mesh_size
154 : : !! [mesh_type]=Type of mesh
155 : : !! If not present, take mesh%mesh_type
156 : : !! [rstep]=Radial step of the mesh (AA parameter above)
157 : : !! If not present, take mesh%rstep
158 : : !! [lstep]=Exponential step of the mesh (BB parameter above)
159 : : !! If not present, take mesh%lstep
160 : : !! Needed only if mesh type is logarithmic.
161 : : !! [r_for_intg]=Mesh size used in integrals computation
162 : : !! If not present, take mesh%r_for_intg
163 : : !! Integrals will be computed up to rr(r_for_intg)
164 : : !! (can be negative for an integration over the whole grid)
165 : : !!
166 : : !! OUTPUT
167 : : !!
168 : : !! SIDE EFFECTS
169 : : !! mesh<pawrad_type>=The object completely initialized (containing radial grid information).
170 : : !! The following quantities are calculated inside the routine:
171 : : !! %stepint = Radial step used to convert any function from the
172 : : !! present grid onto a regular grid in order to integrate it using trapeze method
173 : : !! %rad(mesh_size) = Coordinates of all the points of the mesh.
174 : : !! %radfact(mesh_size) = Factors used to compute radial integrals.
175 : : !! %int_meshsz = Integrals will be computed up to r(int_meshsz)
176 : : !! %simfact(mesh_size) = Factor used to compute radial integrals by the a Simpson scheme
177 : : !! Integral[f] = Sum_i [simfact(i)*f(i)]
178 : : !! %rmax = Max. value of r = rad(mesh_size)
179 : : !!
180 : : !! NOTES
181 : : !! Possible mesh types (mesh%mesh_type)
182 : : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
183 : : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
184 : : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
185 : : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
186 : : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
187 : : !!
188 : : !! SOURCE
189 : :
190 : 5212 : subroutine pawrad_init(mesh,mesh_size,mesh_type,rstep,lstep,r_for_intg)
191 : :
192 : : !Arguments ------------------------------------
193 : : !scalars
194 : : integer,intent(in),optional :: mesh_size,mesh_type
195 : : real(dp),intent(in),optional :: rstep,lstep
196 : : real(dp),intent(in),optional :: r_for_intg
197 : : type(pawrad_type),intent(inout) :: mesh
198 : :
199 : : !Local variables-------------------------------
200 : : !scalars
201 : : integer :: ir,ir_last,isim,mesh_size_,mesh_type_
202 : : real(dp) :: hh,lstep_,rstep_,r_for_intg_
203 : : character(len=500) :: msg
204 : :
205 : : ! *************************************************************************
206 : :
207 : : !@pawrad_type
208 : :
209 : : !Retrieve mesh data
210 [ + + ]: 5212 : mesh_size_ = mesh%mesh_size ; if (present(mesh_size)) mesh_size_ = mesh_size
211 [ + + ]: 5212 : mesh_type_ = mesh%mesh_type ; if (present(mesh_type)) mesh_type_ = mesh_type
212 [ + + ]: 5212 : rstep_ = mesh%rstep ; if (present(rstep)) rstep_ = rstep
213 [ + + ]: 5212 : lstep_ = mesh%lstep ; if (present(lstep)) lstep_ = lstep
214 : :
215 [ + + ]: 5212 : r_for_intg_=-1._dp;if (present(r_for_intg)) r_for_intg_=r_for_intg
216 : :
217 : 5212 : mesh%mesh_size = mesh_size_
218 : 5212 : mesh%mesh_type = mesh_type_
219 : 5212 : mesh%rstep = rstep_
220 : 5212 : mesh%lstep = lstep_
221 [ - + + - : 15636 : LIBPAW_ALLOCATE(mesh%rad ,(mesh%mesh_size))
+ - + - -
+ - + -
+ ]
222 [ - + + - : 10424 : LIBPAW_ALLOCATE(mesh%radfact,(mesh%mesh_size))
- + - + -
+ ]
223 [ - + + - : 10424 : LIBPAW_ALLOCATE(mesh%simfact,(mesh%mesh_size))
- + - + -
+ ]
224 [ + + ]: 7939779 : mesh%simfact=zero
225 : : if (mesh%mesh_type==1) then
226 : 1636 : isim=3
227 : 1636 : mesh%stepint=mesh%rstep
228 : 1636 : mesh%rad(1)=zero;mesh%radfact(1)=one
229 [ + + ]: 3038160 : do ir=2,mesh%mesh_size
230 : 3036524 : mesh%rad(ir) =mesh%rstep*dble(ir-1)
231 : 3038160 : mesh%radfact(ir)=one
232 : : end do
233 : : else if (mesh%mesh_type==2) then
234 : 3427 : isim=3
235 : 3427 : mesh%stepint=mesh%lstep
236 : 3427 : mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep
237 [ + + ]: 4803340 : do ir=2,mesh%mesh_size
238 : 4799913 : mesh%rad(ir) =mesh%rstep*(exp(mesh%lstep*dble(ir-1))-one)
239 : 4803340 : mesh%radfact(ir)=mesh%rad(ir)+mesh%rstep
240 : : end do
241 : : else if (mesh%mesh_type==3) then
242 : 133 : isim=4
243 : 133 : mesh%stepint=mesh%lstep
244 : 133 : mesh%rad(1)=zero;mesh%radfact(1)=zero
245 [ + + ]: 85883 : do ir=2,mesh%mesh_size
246 : 85750 : mesh%rad(ir) =mesh%rstep*exp(mesh%lstep*dble(ir-2))
247 : 85883 : mesh%radfact(ir)=mesh%rad(ir)
248 : : end do
249 : : else if (mesh%mesh_type==4) then
250 : 0 : isim=3
251 : 0 : mesh%lstep=one/dble(mesh%mesh_size)
252 : 0 : mesh%stepint=mesh%lstep
253 : 0 : mesh%rad(1)=zero;mesh%radfact(1)=mesh%rstep
254 [ # # ]: 0 : do ir=2,mesh%mesh_size
255 : 0 : mesh%rad(ir) =-mesh%rstep*log(one-mesh%lstep*dble(ir-1))
256 : 0 : mesh%radfact(ir)=mesh%rstep/(one-mesh%lstep*dble(ir-1))
257 : : end do
258 : : else if (mesh%mesh_type==5) then
259 : 16 : isim=3
260 : 16 : mesh%stepint=mesh%rstep
261 : 16 : mesh%rad(1)=zero;mesh%radfact(1)=1/mesh%lstep
262 [ + + ]: 7184 : do ir=2,mesh%mesh_size
263 : 7168 : mesh%rad(ir) =mesh%rstep*dble(ir-1)/(mesh%lstep-dble(ir-1))
264 : 7184 : mesh%radfact(ir)=(mesh%rstep+mesh%rad(ir))/(mesh%lstep-dble(ir-1))/mesh%rstep
265 : : end do
266 : :
267 : : else ! Other values of mesh_type are not allowed (see psp7in.F90)
268 : 0 : write(msg,'(a,i0)')" Unknown value of mesh_type: ",mesh%mesh_type
269 : 0 : LIBPAW_ERROR(msg)
270 : : end if
271 : :
272 : 5212 : mesh%int_meshsz=mesh%mesh_size
273 [ + + ]: 5212 : if (r_for_intg_>0.d0) then
274 : 2156 : ir=min(pawrad_ifromr(mesh,r_for_intg_),mesh%mesh_size)
275 [ + + ]: 2156 : if (ir<mesh%mesh_size) then
276 [ + + ]: 1826 : if (abs(mesh%rad(ir+1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir+1
277 : : end if
278 [ - + ]: 2156 : if (ir>1) then
279 [ + - ]: 2156 : if (abs(mesh%rad(ir-1)-r_for_intg_)<abs(mesh%rad(ir)-r_for_intg_)) ir=ir-1
280 : : end if
281 : 2156 : mesh%int_meshsz=ir
282 : : end if
283 : :
284 : 5212 : hh=mesh%stepint/3.d0
285 : 5212 : mesh%simfact(mesh%int_meshsz)=hh*mesh%radfact(mesh%int_meshsz)
286 [ + + ]: 10557 : mesh%simfact(1:isim-2)=zero
287 : 5212 : ir_last=1
288 [ - + ]: 5212 : do ir=mesh%int_meshsz,isim,-2
289 : 3529925 : mesh%simfact(ir-1)=4.d0*hh*mesh%radfact(ir-1)
290 : 3529925 : mesh%simfact(ir-2)=2.d0*hh*mesh%radfact(ir-2)
291 [ + + ]: 3529925 : ir_last=ir-2
292 : : end do
293 : 5212 : mesh%simfact(ir_last)=half*mesh%simfact(ir_last)
294 [ + + + + ]: 872652 : if (mesh%int_meshsz<mesh%mesh_size) mesh%simfact(mesh%int_meshsz+1:mesh%mesh_size)=zero
295 : :
296 : 5212 : mesh%rmax=mesh%rad(mesh%mesh_size)
297 : :
298 : 5212 : end subroutine pawrad_init
299 : : !!***
300 : :
301 : : !----------------------------------------------------------------------
302 : :
303 : : !!****f* m_pawrad/pawrad_free_0D
304 : : !! NAME
305 : : !! pawrad_free_0D
306 : : !!
307 : : !! FUNCTION
308 : : !! Frees all memory allocated in the object
309 : : !!
310 : : !! SOURCE
311 : :
312 : 9737 : subroutine pawrad_free_0D(Rmesh)
313 : :
314 : : !Arguments ------------------------------------
315 : : !arrays
316 : : type(pawrad_type),intent(inout) :: Rmesh
317 : :
318 : : !Local variables-------------------------------
319 : :
320 : : ! *************************************************************************
321 : :
322 : : !@Pawrad_type
323 : :
324 [ + + ]: 9737 : if (allocated(Rmesh%rad )) then
325 : 8781 : LIBPAW_DEALLOCATE(Rmesh%rad)
326 : : end if
327 [ + + ]: 9737 : if (allocated(Rmesh%radfact)) then
328 : 8781 : LIBPAW_DEALLOCATE(Rmesh%radfact)
329 : : end if
330 [ + + ]: 9737 : if (allocated(Rmesh%simfact)) then
331 : 8781 : LIBPAW_DEALLOCATE(Rmesh%simfact)
332 : : end if
333 : 9737 : Rmesh%int_meshsz=0
334 : 9737 : Rmesh%mesh_size=0
335 : 9737 : Rmesh%mesh_type=-1
336 : :
337 : 9737 : end subroutine pawrad_free_0D
338 : : !!***
339 : :
340 : : !----------------------------------------------------------------------
341 : :
342 : : !!****f* m_pawrad/pawrad_free_1D
343 : : !! NAME
344 : : !! pawrad_free_1D
345 : : !!
346 : : !! FUNCTION
347 : : !! Destroy all objects in an array of pawrad data structures
348 : : !!
349 : : !! SOURCE
350 : :
351 [ - + ]: 1993 : subroutine pawrad_free_1D(Rmesh)
352 : :
353 : : !Arguments ------------------------------------
354 : : type(pawrad_type),intent(inout) :: Rmesh(:)
355 : :
356 : : !Local variables-------------------------------
357 : : integer :: ii,nn
358 : :
359 : : ! *************************************************************************
360 : :
361 : : !@pawrad_type
362 : :
363 : 1993 : nn=size(Rmesh)
364 [ + + ]: 1993 : if (nn==0) return
365 : :
366 [ + + ]: 2948 : do ii=1,nn
367 : 2948 : call pawrad_free_0D(Rmesh(ii))
368 : : end do
369 : :
370 : : end subroutine pawrad_free_1D
371 : : !!***
372 : :
373 : : !----------------------------------------------------------------------
374 : :
375 : : !!****f* m_pawrad/pawrad_print
376 : : !! NAME
377 : : !! pawrad_print
378 : : !!
379 : : !! FUNCTION
380 : : !! Reports basic info on the object.
381 : : !!
382 : : !! INPUTS
383 : : !! Rmesh<pawrad_type>=Object defining the radial mesh
384 : : !! header=String for the header provided by the user.
385 : : !! [unit]=Unit number for output, defaults to std_out
386 : : !! [prtvol]=Verbosity level, minimal if not specified.
387 : : !! [mode_paral]=Either "COLL" or "PERS". Passed to wrtout. Defaults to "COLL"
388 : : !!
389 : : !! OUTPUT
390 : : !! Only writing.
391 : : !!
392 : : !! SOURCE
393 : :
394 : 0 : subroutine pawrad_print(Rmesh,header,unit,prtvol,mode_paral)
395 : :
396 : : !Arguments ------------------------------------
397 : : integer,intent(in),optional :: prtvol,unit
398 : : character(len=4),intent(in),optional :: mode_paral
399 : : character(len=*),intent(in),optional :: header
400 : : type(pawrad_type),intent(in) :: Rmesh
401 : :
402 : : !Local variables-------------------------------
403 : : !scalars
404 : : integer :: my_unt,my_prtvol
405 : : character(len=4) :: my_mode
406 : : character(len=500) :: msg
407 : :
408 : : ! *************************************************************************
409 : :
410 : : !@pawrad_type
411 : 0 : my_unt =std_out; if (PRESENT(unit )) my_unt =unit
412 [ # # ]: 0 : my_prtvol=0 ; if (PRESENT(prtvol )) my_prtvol=prtvol
413 [ # # ]: 0 : my_mode ='COLL' ; if (PRESENT(mode_paral)) my_mode =mode_paral
414 : :
415 : 0 : msg=ch10//' ==== Info on the Radial Mesh ==== '
416 [ # # # # : 0 : if (PRESENT(header)) msg=ch10//' ==== '//TRIM(ADJUSTL(header))//' ==== '
# # ]
417 : 0 : call wrtout(my_unt,msg,my_mode)
418 : :
419 : 0 : SELECT CASE (Rmesh%mesh_type)
420 : :
421 : : CASE (RMESH_LINEAR)
422 : : write(msg,'(a,i4,a,g12.5)')&
423 : 0 : & ' - Linear mesh: r(i)=step*(i-1), size=',Rmesh%mesh_size,', step=',Rmesh%rstep
424 : :
425 : : CASE (RMESH_LOG1)
426 : : write(msg,'(a,i4,2(a,g12.5))')&
427 : 0 : & ' - Logarithimc mesh: r(i)=AA*[exp(BB*(i-1))-1], size=',Rmesh%mesh_size,', AA=',Rmesh%rstep,' BB=',Rmesh%lstep
428 : :
429 : : CASE (RMESH_LOG2)
430 : : write(msg,'(a,i4,2(a,g12.5))')&
431 : 0 : & ' - Logarithimc mesh: r(i)=AA*exp(BB*(i-2)), size=',Rmesh%mesh_size,', AA=',Rmesh%rstep,' BB=',Rmesh%lstep
432 : :
433 : : CASE (RMESH_LOG3)
434 : : write(msg,'(a,i1,a,i4,a,g12.5)')&
435 : 0 : & ' - Logarithimc mesh: r(i)=-AA*ln(1-(i-1)/n), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep
436 : :
437 : : CASE (RMESH_NL)
438 : : write(msg,'(a,i1,a,i4,a,g12.5)')&
439 : 0 : & ' - Non-linear mesh: r(i)=-AA*i/(n-i), n=size=',Rmesh%mesh_size,', AA=',Rmesh%rstep
440 : :
441 : : CASE DEFAULT
442 : 0 : msg = ' Unknown mesh type! Action : check your pseudopotential or input file.'
443 [ # # # # : 0 : LIBPAW_ERROR(msg)
# # ]
444 : : END SELECT
445 : :
446 : 0 : call wrtout(my_unt,msg,my_mode)
447 : :
448 [ # # ]: 0 : if (my_prtvol>1) then
449 : 0 : write(msg,'(a,i4)')' Mesh size for integrals = ',Rmesh%int_meshsz
450 : 0 : call wrtout(my_unt,msg,my_mode)
451 : 0 : write(msg,'(a,g12.5)')' rmax=rad(mesh_size) = ',Rmesh%rmax
452 : 0 : call wrtout(my_unt,msg,my_mode)
453 : 0 : write(msg,'(a,g12.5)')' Value of stepint = ',Rmesh%stepint
454 : 0 : call wrtout(my_unt,msg,my_mode)
455 : : end if
456 : :
457 [ # # ]: 0 : end subroutine pawrad_print
458 : : !!***
459 : :
460 : : !----------------------------------------------------------------------
461 : :
462 : : !!****f* m_pawrad/pawrad_isame
463 : : !! NAME
464 : : !! pawrad_isame
465 : : !!
466 : : !! FUNCTION
467 : : !! Check two radial meshes, returns a logical flag defining
468 : : !! whether the meshes have the same equation and an integer
469 : : !! flag
470 : : !!
471 : : !! INPUTS
472 : : !! Rmesh1,Rmesh2<pawrad_type>=The two radial meshes.
473 : : !!
474 : : !! OUTPUT
475 : : !! hasameq=.true. if the two meshes are defined by the same equation.
476 : : !! whichdenser=
477 : : !! * 0 if meshes are not compatible
478 : : !! * 1 if Rmesh1 is denser than Rmesh2
479 : : !! * 2 if Rmesh2 is denser than Rmesh1
480 : : !!
481 : : !! SOURCE
482 : :
483 : 2 : subroutine pawrad_isame(Rmesh1,Rmesh2,hasameq,whichdenser)
484 : :
485 : : !Arguments ------------------------------------
486 : : integer,intent(out) :: whichdenser
487 : : logical,intent(out) :: hasameq
488 : : type(pawrad_type),intent(in) :: Rmesh1
489 : : type(pawrad_type),intent(in) :: Rmesh2
490 : :
491 : : !Local variables-------------------------------
492 : : character(len=50) :: msg
493 : :
494 : : ! *************************************************************************
495 : :
496 : : !@pawrad_type
497 : :
498 : 2 : whichdenser =0 ; hasameq=.FALSE.
499 : :
500 [ - + ]: 2 : if (Rmesh1%mesh_type /= Rmesh2%mesh_type) RETURN
501 : :
502 : 0 : SELECT CASE (Rmesh1%mesh_type)
503 : :
504 : : CASE (RMESH_LINEAR) !check for linear meshes
505 : 0 : hasameq = (Rmesh1%rstep == Rmesh2%rstep)
506 : :
507 : : CASE (RMESH_LOG1,& !check for logarithmic meshes
508 : : & RMESH_LOG2,&
509 : : & RMESH_LOG3)
510 : :
511 : : hasameq = ( Rmesh1%rstep == Rmesh2%rstep &
512 [ - + - - ]: 4 : & .and.Rmesh1%lstep == Rmesh2%lstep )
513 : :
514 : : CASE (RMESH_NL) !check for linear meshes
515 : 0 : hasameq = (Rmesh1%rstep == Rmesh2%rstep)
516 : :
517 : : CASE DEFAULT
518 : 0 : msg='Unknown mesh type'
519 [ - + - - ]: 2 : LIBPAW_ERROR(msg)
520 : :
521 : : END SELECT
522 : :
523 : : ! === If meshes have same equation, check whether they are equal ===
524 : : ! * Note that also int_meshsz must be equal
525 [ - + ]: 2 : if (hasameq) then
526 : 0 : whichdenser= 1
527 [ # # ]: 0 : if (Rmesh2%mesh_size > Rmesh1%mesh_size ) whichdenser = 2
528 : : !if (Rmesh1%int_meshsz == Rmesh2%int_meshsz) whichdenser = 2
529 : : end if
530 : :
531 : : end subroutine pawrad_isame
532 : : !!***
533 : :
534 : : !----------------------------------------------------------------------
535 : :
536 : : !!****f* m_pawrad/pawrad_copy
537 : : !! NAME
538 : : !! pawrad_copy
539 : : !!
540 : : !! FUNCTION
541 : : !! Copy one radial mesh (in a generalized format) to another
542 : : !!
543 : : !! INPUTS
544 : : !! mesh1 <type(pawrad_type)>=data containing radial grid information of input mesh
545 : : !!
546 : : !! OUTPUT
547 : : !! mesh2 <type(pawrad_type)>=data containing radial grid information of output mesh
548 : : !!
549 : : !! NOTES
550 : : !! Possible mesh types (mesh%mesh_type)
551 : : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
552 : : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
553 : : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
554 : : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
555 : : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
556 : : !!
557 : : !! SOURCE
558 : :
559 [ - + - + : 3382 : subroutine pawrad_copy(mesh1,mesh2)
- + ]
560 : :
561 : : !Arguments ------------------------------------
562 : : !scalars
563 : : type(pawrad_type),intent(in) :: mesh1
564 : : type(pawrad_type),intent(out) :: mesh2
565 : :
566 : : !Local variables-------------------------------
567 : : !scalars
568 : : integer :: ir
569 : :
570 : : ! *************************************************************************
571 : :
572 : 3382 : mesh2%mesh_type =mesh1%mesh_type
573 : 3382 : mesh2%mesh_size =mesh1%mesh_size
574 : 3382 : mesh2%int_meshsz=mesh1%int_meshsz
575 : 3382 : mesh2%lstep =mesh1%lstep
576 : 3382 : mesh2%rstep =mesh1%rstep
577 : 3382 : mesh2%stepint =mesh1%stepint
578 : 3382 : mesh2%rmax =mesh1%rmax
579 : :
580 [ - + + - : 10146 : LIBPAW_ALLOCATE(mesh2%rad,(mesh1%mesh_size))
+ - + - -
+ - + ]
581 [ - + + - : 6764 : LIBPAW_ALLOCATE(mesh2%radfact,(mesh1%mesh_size))
- + - + ]
582 [ - + + - : 6764 : LIBPAW_ALLOCATE(mesh2%simfact,(mesh1%mesh_size))
- + - + ]
583 [ + + ]: 3805553 : do ir=1,mesh1%mesh_size
584 : 3802171 : mesh2%rad(ir) =mesh1%rad(ir)
585 : 3802171 : mesh2%radfact(ir)=mesh1%radfact(ir)
586 : 3805553 : mesh2%simfact(ir)=mesh1%simfact(ir)
587 : : end do
588 : :
589 : 3382 : end subroutine pawrad_copy
590 : : !!***
591 : :
592 : : !----------------------------------------------------------------------
593 : :
594 : : !!****f* m_pawrad/pawrad_deducer0
595 : : !! NAME
596 : : !! pawrad_deducer0
597 : : !!
598 : : !! FUNCTION
599 : : !! Extrapolate r=0 value of a function from values near r=0
600 : : !! using a 3 points formula
601 : : !!
602 : : !! INPUTS
603 : : !! funcsz=size of array func
604 : : !! radmesh <type(pawrad_type)>=data containing radial grid information
605 : : !!
606 : : !! SIDE EFFECTS
607 : : !! func(funcsz)=array containing values of function to extrapolate
608 : : !!
609 : : !! SOURCE
610 : :
611 : 11394323 : subroutine pawrad_deducer0(func,funcsz,radmesh)
612 : :
613 : : !Arguments ------------------------------------
614 : : !scalars
615 : : integer,intent(in) :: funcsz
616 : : type(pawrad_type),intent(in) :: radmesh
617 : : !arrays
618 : : real(dp),intent(inout) :: func(funcsz)
619 : :
620 : : ! *************************************************************************
621 : :
622 [ + + + + ]: 11394323 : if (radmesh%mesh_type==1.or.radmesh%mesh_type==2.or.radmesh%mesh_type==4.or.radmesh%mesh_type==5) then
623 : 11284952 : func(1)=func(4)+3*(func(2)-func(3))
624 [ + - ]: 109371 : else if (radmesh%mesh_type==3) then
625 : 109371 : func(1)=func(4)+exp(two*radmesh%lstep)/(exp(radmesh%lstep)-one)*(func(2)-func(3))
626 : : end if
627 : :
628 : 11394323 : end subroutine pawrad_deducer0
629 : : !!***
630 : :
631 : : !----------------------------------------------------------------------
632 : :
633 : : !!****f* m_pawrad/pawrad_bcast
634 : : !! NAME
635 : : !! pawrad_bcast
636 : : !!
637 : : !! FUNCTION
638 : : !! Communicate pawrad data over a given MPI communicator
639 : : !!
640 : : !! INPUTS
641 : : !! comm_mpi= communicator used to broadcast data
642 : : !!
643 : : !! SIDE EFFECTS
644 : : !! pawrad=<type pawrad_type>= a radial mesh datastructure for PAW
645 : : !!
646 : : !! SOURCE
647 : :
648 : 284 : subroutine pawrad_bcast(pawrad,comm_mpi)
649 : :
650 : : !Arguments ------------------------------------
651 : : !scalars
652 : : integer,intent(in) :: comm_mpi
653 : : type(pawrad_type),intent(inout) :: pawrad
654 : :
655 : : !Local variables-------------------------------
656 : : !scalars
657 : : integer :: ierr,indx,me,nn,isz1
658 : : integer :: if_rad,if_radfact,if_simfact !flags used to communicate
659 : : character(len=500) :: msg
660 : : !arrays
661 : 284 : integer,allocatable :: list_int(:)
662 : 284 : real(dp),allocatable :: list_dpr(:)
663 : :
664 : : !*************************************************************************
665 : :
666 : 284 : me=xmpi_comm_rank(comm_mpi)
667 : :
668 : : !Initializations
669 : 284 : if_rad=0; if_radfact=0; if_simfact=0
670 : :
671 : : !calculate the size of the reals
672 [ + + ]: 284 : if(me==0) then
673 [ + - ]: 72 : if (allocated(pawrad%rad)) then
674 : 72 : if_rad=1 !communicate rad
675 : 72 : isz1=size(pawrad%rad)
676 [ - + ]: 72 : if(isz1/=pawrad%mesh_size) then
677 : 0 : msg='rad: sz1 /= pawrad%mesh_size (1)'
678 : 0 : LIBPAW_BUG(msg)
679 : : end if
680 : : end if
681 [ + - ]: 72 : if (allocated(pawrad%radfact)) then
682 : 72 : if_radfact=1 !communicate radfact
683 : 72 : isz1=size(pawrad%radfact)
684 [ - + ]: 72 : if(isz1/=pawrad%mesh_size) then
685 : 0 : msg='radfact: sz1 /= pawrad%mesh_size (2)'
686 : 0 : LIBPAW_BUG(msg)
687 : : end if
688 : : end if
689 [ + - ]: 72 : if (allocated(pawrad%simfact)) then
690 : 72 : if_simfact=1 !communicate simfact
691 : 72 : isz1=size(pawrad%simfact)
692 [ - + ]: 72 : if(isz1/=pawrad%mesh_size) then
693 : 0 : msg='simfact: sz1 /= pawrad%mesh_size (3)'
694 : 0 : LIBPAW_BUG(msg)
695 : : end if
696 : : end if
697 : : end if
698 : :
699 : : !Brodcast the integers
700 [ - + ]: 284 : LIBPAW_ALLOCATE(list_int,(6))
701 [ + + ]: 284 : if(me==0) then
702 : 72 : list_int(1)=pawrad%int_meshsz
703 : 72 : list_int(2)=pawrad%mesh_size
704 : 72 : list_int(3)=pawrad%mesh_type
705 : 72 : list_int(4)=if_rad
706 : 72 : list_int(5)=if_radfact
707 : 72 : list_int(6)=if_simfact
708 : : end if
709 : 284 : call xmpi_bcast(list_int,0,comm_mpi,ierr)
710 [ + + ]: 284 : if(me/=0) then
711 : 212 : pawrad%int_meshsz =list_int(1)
712 : 212 : pawrad%mesh_size =list_int(2)
713 : 212 : pawrad%mesh_type =list_int(3)
714 : 212 : if_rad=list_int(4)
715 : 212 : if_radfact=list_int(5)
716 : 212 : if_simfact=list_int(6)
717 : : end if
718 [ - + ]: 284 : LIBPAW_DEALLOCATE(list_int)
719 : :
720 : : !Broadcast the reals
721 : 284 : nn=4+pawrad%mesh_size*(if_rad+if_radfact+if_simfact)
722 [ - + + - : 852 : LIBPAW_ALLOCATE(list_dpr,(nn))
+ - + - -
+ - + ]
723 [ + + ]: 284 : if(me==0) then
724 : 72 : list_dpr(1)=pawrad%lstep
725 : 72 : list_dpr(2)=pawrad%rmax
726 : 72 : list_dpr(3)=pawrad%rstep
727 : 72 : list_dpr(4)=pawrad%stepint
728 : 72 : indx=5
729 [ + - ]: 72 : if (if_rad==1) then
730 : 72 : isz1=pawrad%mesh_size
731 [ + + ]: 83774 : list_dpr(indx:indx+isz1-1)=pawrad%rad(1:isz1)
732 : 72 : indx=indx+isz1
733 : : end if
734 [ + - ]: 72 : if (if_radfact==1) then
735 : 72 : isz1=pawrad%mesh_size
736 [ + + ]: 83774 : list_dpr(indx:indx+isz1-1)=pawrad%radfact(1:isz1)
737 : : indx=indx+isz1
738 : : end if
739 [ + - ]: 72 : if (if_simfact==1) then
740 : 72 : isz1=pawrad%mesh_size
741 [ + + ]: 83774 : list_dpr(indx:indx+isz1-1)=pawrad%simfact(1:isz1)
742 : 284 : indx=indx+isz1
743 : : end if
744 : : end if
745 : 284 : call xmpi_bcast(list_dpr,0,comm_mpi,ierr)
746 [ + + ]: 284 : if(me/=0) then
747 : 212 : pawrad%lstep=list_dpr(1)
748 : 212 : pawrad%rmax=list_dpr(2)
749 : 212 : pawrad%rstep=list_dpr(3)
750 : 212 : pawrad%stepint=list_dpr(4)
751 : 212 : indx=5
752 : : ! Deallocate all arrays:
753 [ + + ]: 212 : if (allocated(pawrad%rad)) then
754 : 24 : LIBPAW_DEALLOCATE(pawrad%rad)
755 : : end if
756 [ + + ]: 212 : if (allocated(pawrad%radfact)) then
757 : 24 : LIBPAW_DEALLOCATE(pawrad%radfact)
758 : : end if
759 [ + + ]: 212 : if (allocated(pawrad%simfact)) then
760 : 24 : LIBPAW_DEALLOCATE(pawrad%simfact)
761 : : end if
762 : : ! Communicate if flag is set to 1:
763 [ + - ]: 212 : if(if_rad==1) then
764 : 212 : isz1=pawrad%mesh_size
765 [ - + + - : 636 : LIBPAW_ALLOCATE(pawrad%rad,(isz1))
+ - + - -
+ - + -
+ ]
766 [ + + ]: 255450 : pawrad%rad(1:isz1)=list_dpr(indx:indx+isz1-1)
767 : 212 : indx=indx+isz1
768 : : end if
769 [ + - ]: 212 : if(if_radfact==1) then
770 : 212 : isz1=pawrad%mesh_size
771 [ - + + - : 636 : LIBPAW_ALLOCATE(pawrad%radfact,(isz1))
+ - + - -
+ - + -
+ ]
772 [ + + ]: 255450 : pawrad%radfact(1:isz1)=list_dpr(indx:indx+isz1-1)
773 : : indx=indx+isz1
774 : : end if
775 [ + - ]: 212 : if(if_simfact==1) then
776 : 212 : isz1=pawrad%mesh_size
777 [ - + + - : 636 : LIBPAW_ALLOCATE(pawrad%simfact,(isz1))
+ - + - -
+ - + -
+ ]
778 [ + + ]: 255450 : pawrad%simfact(1:isz1)=list_dpr(indx:indx+isz1-1)
779 : 284 : indx=indx+isz1
780 : : end if
781 : : end if
782 [ - + ]: 284 : LIBPAW_DEALLOCATE(list_dpr)
783 : :
784 : 284 : end subroutine pawrad_bcast
785 : : !!***
786 : :
787 : : !----------------------------------------------------------------------
788 : :
789 : : !!****f* m_pawrad/simp_gen
790 : : !! NAME
791 : : !! simp_gen
792 : : !!
793 : : !! FUNCTION
794 : : !! Do integral on a given (generalized) grid using Simpson rule
795 : : !!
796 : : !! INPUTS
797 : : !! radmesh <type(pawrad_type)>=data containing radial grid information
798 : : !! func(:)=integrand values
799 : : !! r_for_intg=upper boundary for (future) integration over the radial grid
800 : : !! (can be negative for an integration over the whole grid)
801 : : !!
802 : : !! OUTPUT
803 : : !! intg=resulting integral by Simpson rule
804 : : !!
805 : : !! NOTES
806 : : !! Possible mesh types (radmesh%mesh_type)
807 : : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
808 : : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
809 : : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
810 : : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
811 : : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
812 : : !!
813 : : !! SOURCE
814 : :
815 [ - + ]: 89436765 : subroutine simp_gen(intg,func,radmesh,r_for_intg)
816 : :
817 : : #if defined HAVE_AVX_SAFE_MODE
818 : : !DEC$ NOOPTIMIZE
819 : : #endif
820 : :
821 : : !Arguments ------------------------------------
822 : : !scalars
823 : : real(dp),intent(out) :: intg
824 : : real(dp),intent(in),optional :: r_for_intg
825 : : type(pawrad_type),intent(in) :: radmesh
826 : : !arrays
827 : : real(dp),intent(in) :: func(:)
828 : :
829 : : !Local variables-------------------------------
830 : : !scalars
831 : : integer :: ii,int_meshsz,ir,ir_last,isim,nn
832 : : real(dp) :: hh,resid,simp
833 : 89436765 : real(dp),allocatable :: simfact(:)
834 : : character(len=500) :: msg
835 : :
836 : : ! *************************************************************************
837 : :
838 [ + + ]: 89436765 : if (present(r_for_intg)) then
839 [ + - ]: 3333 : if (r_for_intg>0.d0) then
840 : 3333 : ir=min(pawrad_ifromr(radmesh,r_for_intg),radmesh%mesh_size)
841 [ + + ]: 3333 : if (ir<radmesh%mesh_size) then
842 [ + + ]: 2708 : if (abs(radmesh%rad(ir+1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir+1
843 : : end if
844 [ - + ]: 3333 : if (ir>1) then
845 [ + - ]: 3333 : if (abs(radmesh%rad(ir-1)-r_for_intg)<abs(radmesh%rad(ir)-r_for_intg)) ir=ir-1
846 : : end if
847 : 3333 : int_meshsz=ir
848 : : else
849 : 0 : int_meshsz=radmesh%mesh_size
850 : : end if
851 [ + - - + ]: 3333 : if (int_meshsz>radmesh%mesh_size.or.int_meshsz>size(func)) then
852 : 0 : write(msg,'(3(a,i4))')"int_meshsz= ",int_meshsz," > mesh_size=",radmesh%mesh_size,&
853 : 0 : & ", size(func)=",size(func)
854 : 0 : LIBPAW_BUG(msg)
855 : : end if
856 [ + + ]: 3333 : isim=3; if (radmesh%mesh_type==3)isim=4
857 [ - + + - : 9999 : LIBPAW_ALLOCATE(simfact,(radmesh%mesh_size))
+ - + - -
+ - + ]
858 : 3333 : hh=radmesh%stepint/3.d0
859 : 3333 : simfact(int_meshsz)=hh*radmesh%radfact(int_meshsz)
860 [ + + ]: 6859 : simfact(1:isim-2)=zero
861 : 3333 : ir_last=1
862 [ + - ]: 3333 : do ir=int_meshsz,isim,-2
863 : 1576195 : simfact(ir-1)=4.d0*hh*radmesh%radfact(ir-1)
864 : 1576195 : simfact(ir-2)=2.d0*hh*radmesh%radfact(ir-2)
865 [ + + ]: 1576195 : ir_last=ir-2
866 : : end do
867 : 3333 : simfact(ir_last)=half*simfact(ir_last)
868 [ + + + + ]: 1036370 : if (int_meshsz<radmesh%mesh_size) simfact(int_meshsz+1:radmesh%mesh_size)=zero
869 : :
870 : 3333 : nn=int_meshsz
871 : 3333 : simp=zero
872 [ + + ]: 3160958 : do ii=1,nn
873 : 3160958 : simp=simp+func(ii)*simfact(ii)
874 : : end do
875 : 3333 : LIBPAW_DEALLOCATE(simfact)
876 : :
877 : : else
878 [ - + ]: 89433432 : if (radmesh%int_meshsz>size(func)) then
879 : 0 : write(msg,'(2(a,i4))')"int_meshsz= ",radmesh%int_meshsz," > size(func)=",size(func)
880 : 0 : LIBPAW_BUG(msg)
881 : : end if
882 : 89433432 : nn=radmesh%int_meshsz
883 : 89433432 : simp=zero
884 [ + + ]: >11033*10^7 : do ii=1,nn
885 : >11033*10^7 : simp=simp+func(ii)*radmesh%simfact(ii)
886 : : end do
887 : : end if
888 : :
889 : 89436765 : resid=zero
890 [ + + ]: 89436765 : if (radmesh%mesh_type==3) then
891 : 1302684 : resid=half*(func(2)+func(1))*(radmesh%rad(2)-radmesh%rad(1))
892 [ + + ]: 1302684 : if (mod(nn,2)==1) resid=resid+radmesh%stepint/3.d0*(1.25d0*func(2)*radmesh%radfact(2) &
893 : 99204 : & +2.d0*func(3)*radmesh%radfact(3)-0.25d0*func(4)*radmesh%radfact(4))
894 [ + + ]: 88134081 : else if (mod(nn,2)==0) then
895 : : resid=radmesh%stepint/3.d0*(1.25d0*func(1)*radmesh%radfact(1)+2.d0*func(2)*radmesh%radfact(2) &
896 : 51926117 : & -0.25d0*func(3)*radmesh%radfact(3))
897 : : end if
898 : :
899 : 89436765 : intg=simp+resid
900 : :
901 : 89436765 : end subroutine simp_gen
902 : : !!***
903 : :
904 : : !----------------------------------------------------------------------
905 : :
906 : : !!****f* m_pawrad/nderiv_gen
907 : : !! NAME
908 : : !! nderiv_gen
909 : : !!
910 : : !! FUNCTION
911 : : !! Do corrected first (and -if requested- second) derivation on a given (generalized) grid.
912 : : !! This routine interfaces nderiv_lin (derivation on a linear grid).
913 : : !!
914 : : !! INPUTS
915 : : !! func(:)=input function
916 : : !! radmesh <type(pawrad_type)>=data containing radial grid information
917 : : !!
918 : : !! OUTPUT
919 : : !! der(:)= 1st derivative of input function
920 : : !! [der2(:)]= -- optional -- 2nd derivative of input function
921 : : !!
922 : : !! NOTES
923 : : !! Possible mesh types (radmesh%mesh_type)
924 : : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
925 : : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
926 : : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
927 : : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
928 : : !!
929 : : !! SOURCE
930 : :
931 [ - + - + : 7652406 : subroutine nderiv_gen(der,func,radmesh,der2)
- + - - -
- ]
932 : :
933 : : !Arguments ------------------------------------
934 : : !scalars
935 : : type(pawrad_type),intent(in) :: radmesh
936 : : !arrays
937 : : real(dp),intent(in) :: func(:)
938 : : real(dp),intent(out) :: der(:)
939 : : real(dp),optional,intent(out) :: der2(:)
940 : :
941 : : !Local variables-------------------------------
942 : : !scalars
943 : : integer :: msz
944 : : logical :: compute_2der
945 : : character(len=500) :: msg
946 : :
947 : : ! *************************************************************************
948 : :
949 : 7652406 : msz=size(func)
950 [ + - - + ]: 7652406 : if (size(der)/=msz.or.msz>radmesh%mesh_size) then
951 : 0 : msg='wrong sizes for in/out arrays!'
952 : 0 : LIBPAW_BUG(msg)
953 : : end if
954 : :
955 : 7652406 : compute_2der=(present(der2))
956 : :
957 [ + + + - : 7652406 : if (radmesh%mesh_type==1) then
- - ]
958 : :
959 [ - + - - : 138 : call nderiv_lin(radmesh%rstep,func,der,msz,1)
- - - + -
- - + - +
- - ]
960 [ - + ]: 138 : if (compute_2der) then
961 [ # # # # : 0 : call nderiv_lin(radmesh%rstep,func,der2,msz,2)
# # # # #
# # # # #
# # ]
962 : : end if
963 : :
964 : : else if (radmesh%mesh_type==2) then
965 : :
966 [ - + - - : 7652132 : call nderiv_lin(radmesh%lstep,func,der,msz,1)
- - - + -
- - + - +
- - ]
967 [ + + ]: 9632392065 : der(1:msz)=der(1:msz)/radmesh%radfact(1:msz)
968 [ - + ]: 7652132 : if (compute_2der)then
969 [ # # # # : 0 : call nderiv_lin(radmesh%lstep,func,der2,msz,2)
# # # # #
# # # # #
# # ]
970 [ # # ]: 0 : der2(1:msz)=(der2(1:msz)/radmesh%radfact(1:msz)-der(1:msz))/radmesh%radfact(1:msz)
971 : : end if
972 : :
973 : : else if (radmesh%mesh_type==3) then
974 : :
975 [ - + - - : 136 : call nderiv_lin(radmesh%lstep,func(2:msz),der(2:msz),msz-1,1)
- - - + -
- - + - +
- - ]
976 [ + + ]: 47600 : der(2:msz)=der(2:msz)/radmesh%radfact(2:msz)
977 [ - + - - : 136 : call pawrad_deducer0(der,msz,radmesh)
- - - + -
- ]
978 [ - + ]: 136 : if (compute_2der)then
979 [ # # # # : 0 : call nderiv_lin(radmesh%lstep,func(2:msz),der2(2:msz),msz-1,2)
# # # # #
# # # # #
# # ]
980 [ # # ]: 0 : der2(2:msz)=(der2(2:msz)/radmesh%radfact(2:msz)-der(2:msz))/radmesh%radfact(2:msz)
981 [ # # # # : 0 : call pawrad_deducer0(der2,msz,radmesh)
# # # # #
# ]
982 : : end if
983 : :
984 : : else if (radmesh%mesh_type==4) then
985 : :
986 [ # # # # : 0 : call nderiv_lin(radmesh%lstep,func,der,msz,1)
# # # # #
# # # # #
# # ]
987 [ # # ]: 0 : der(1:msz)=der(1:msz)/radmesh%radfact(1:msz)
988 [ # # ]: 0 : if (compute_2der)then
989 [ # # # # : 0 : call nderiv_lin(radmesh%lstep,func,der2,msz,2)
# # # # #
# # # # #
# # ]
990 [ # # ]: 0 : der2(1:msz)=der2(1:msz)/radmesh%radfact(1:msz)**2-der(1:msz)/radmesh%rstep
991 : : end if
992 : :
993 : : else if (radmesh%mesh_type==5) then
994 : :
995 [ # # # # : 0 : call nderiv_lin(one,func,der,msz,1)
# # # # #
# # # # #
# # ]
996 [ # # ]: 0 : der(1:msz)=der(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep)
997 [ # # ]: 0 : if (compute_2der)then
998 [ # # # # : 0 : call nderiv_lin(one,func,der2,msz,2)
# # # # #
# # # # #
# # ]
999 : : der2(1:msz)=der2(1:msz)/(radmesh%radfact(1:msz)*radmesh%rstep)**2-two*der(1:msz)/&
1000 [ # # ]: 0 : & (radmesh%rstep+radmesh%rad(1:msz))
1001 : : end if
1002 : : end if
1003 : :
1004 : 7652406 : end subroutine nderiv_gen
1005 : : !!***
1006 : :
1007 : : !----------------------------------------------------------------------
1008 : :
1009 : : !!****f* m_pawrad/nderiv_lin
1010 : : !! NAME
1011 : : !! nderiv_lin
1012 : : !!
1013 : : !! FUNCTION
1014 : : !! Do corrected first (and -if requested- second) derivation on a given LINEAR grid.
1015 : : !!
1016 : : !! INPUTS
1017 : : !! hh= radial step
1018 : : !! ndim= radial mesh size
1019 : : !! yy(ndim)= input function
1020 : : !! norder= order of derivation (1 or 2)
1021 : : !!
1022 : : !! OUTPUT
1023 : : !! zz(ndim)= first or second derivative of y
1024 : : !!
1025 : : !! SOURCE
1026 : :
1027 : 7652406 : subroutine nderiv_lin(hh,yy,zz,ndim,norder)
1028 : :
1029 : : !Arguments ---------------------------------------------
1030 : : !scalars
1031 : : integer,intent(in) :: ndim,norder
1032 : : real(dp),intent(in) :: hh
1033 : : !arrays
1034 : : real(dp),intent(in) :: yy(ndim)
1035 : : real(dp),intent(out) :: zz(ndim)
1036 : :
1037 : : !Local variables ---------------------------------------
1038 : : !scalars
1039 : : integer :: ier,ii
1040 : : real(dp) :: aa,bb,cc,h1,y1
1041 : :
1042 : : ! *************************************************************************
1043 : :
1044 : : !Initialization (common to 1st and 2nd derivative)
1045 : 7652406 : h1=one/(12.d0*hh)
1046 : 7652406 : y1=yy(ndim-4)
1047 : :
1048 : : !FIRST DERIVATIVE
1049 : : !================
1050 [ + - ]: 7652406 : if (norder==1) then
1051 : :
1052 : : ! Prepare differentiation loop
1053 : 7652406 : bb=h1*(-25.d0*yy(1)+48.d0*yy(2)-36.d0*yy(3)+16.d0*yy(4)-3.d0*yy(5))
1054 : 7652406 : cc=h1*(-3.d0*yy(1)-10.d0*yy(2)+18.d0*yy(3)-6.d0*yy(4)+yy(5))
1055 : : ! Start differentiation loop
1056 [ + + ]: 9601898407 : do ii=5,ndim
1057 : 9594246001 : aa=bb;bb=cc
1058 : 9594246001 : cc=h1*(yy(ii-4)-yy(ii)+8.d0*(yy(ii-1)-yy(ii-3)))
1059 : 9601898407 : zz(ii-4)=aa
1060 : : end do
1061 : : ! Normal exit
1062 : 7652406 : ier=0
1063 : 7652406 : aa=h1*(-y1+6.d0*yy(ndim-3)-18.d0*yy(ndim-2)+10.d0*yy(ndim-1)+3.d0*yy(ndim))
1064 : 7652406 : zz(ndim)=h1*(3.d0*y1-16.d0*yy(ndim-3)+36.d0*yy(ndim-2) -48.d0*yy(ndim-1)+25.d0*yy(ndim))
1065 : 7652406 : zz(ndim-1)=aa
1066 : 7652406 : zz(ndim-2)=cc
1067 : 7652406 : zz(ndim-3)=bb
1068 : :
1069 : : ! SECOND DERIVATIVE
1070 : : ! =================
1071 : : else
1072 : 0 : h1=h1/hh
1073 : : ! Prepare differentiation loop
1074 : 0 : bb=h1*(35.d0*yy(1)-104.d0*yy(2)+114.d0*yy(3)-56.d0*yy(4)+11.d0*yy(5))
1075 : 0 : cc=h1*(11.d0*yy(1)-20.d0*yy(2)+6.d0*yy(3)+4.d0*yy(4)-yy(5))
1076 : : ! Start differentiation loop
1077 [ # # ]: 0 : do ii=5,ndim
1078 : 0 : aa=bb;bb=cc
1079 : 0 : cc=h1*(-yy(ii-4)-yy(ii)+16.d0*(yy(ii-1)+yy(ii-3))-30.d0*yy(ii-2))
1080 : 0 : zz(ii-4)=aa
1081 : : end do
1082 : : ! Normal exit
1083 : 0 : ier=0
1084 : 0 : aa=h1*(-y1+4.d0*yy(ndim-3)+6.d0*yy(ndim-2)-20.d0*yy(ndim-1)+11.d0*yy(ndim))
1085 : 0 : zz(ndim)=h1*(11.d0*y1-56.d0*yy(ndim-3)+114.d0*yy(ndim-2) -104.d0*yy(ndim-1)+35.d0*yy(ndim))
1086 : 0 : zz(ndim-1)=aa
1087 : 0 : zz(ndim-2)=cc
1088 : 0 : zz(ndim-3)=bb
1089 : :
1090 : : end if !norder
1091 : :
1092 : 7652406 : end subroutine nderiv_lin
1093 : : !!***
1094 : :
1095 : : !----------------------------------------------------------------------
1096 : :
1097 : : !!****f* m_pawrad/bound_deriv
1098 : : !! NAME
1099 : : !! bound_deriv
1100 : : !!
1101 : : !! FUNCTION
1102 : : !! Computes derivatives of a function a boundaries of interval (first and last derivative)
1103 : : !!
1104 : : !! INPUTS
1105 : : !! func(n)= array containing function
1106 : : !! mesh <type(pawrad_type)>= radial mesh and related data
1107 : : !! nn= size of intervall
1108 : : !!
1109 : : !! OUTPUT
1110 : : !! yp1,ypn= derivatives of func at r(1) and r(n)
1111 : : !!
1112 : : !! SOURCE
1113 : :
1114 : 104 : subroutine bound_deriv(func,mesh,nn,yp1,ypn)
1115 : :
1116 : : !Arguments----------------------
1117 : : integer, intent(in) :: nn
1118 : : real(dp), intent(in) :: func(nn)
1119 : : real(dp), intent(out) :: yp1,ypn
1120 : : type(pawrad_type),intent(in) :: mesh
1121 : :
1122 : : !*************************************************************************
1123 : :
1124 [ + - ]: 104 : if (mesh%radfact(1)>zero) then
1125 : : yp1=1._dp/12._dp/mesh%stepint/mesh%radfact(1) &
1126 : 104 : & *(-25._dp*func(1)+48._dp*func(2)-36._dp*func(3)+16._dp*func(4)-3._dp*func(5))
1127 : : else
1128 : 0 : yp1=(func(2)-func(1))/(mesh%rad(2)-mesh%rad(1))
1129 : : end if
1130 : : ypn=1._dp/12._dp/mesh%stepint &
1131 : : & *( 3._dp*func(nn-4)-16._dp*func(nn-3)+36._dp*func(nn-2)-48._dp*func(nn-1) &
1132 : 104 : & +25._dp*func(nn))/mesh%radfact(nn)
1133 : :
1134 : 104 : end subroutine bound_deriv
1135 : :
1136 : : !!***
1137 : :
1138 : : !----------------------------------------------------------------------
1139 : :
1140 : : !!****f* m_pawrad/poisson
1141 : : !! NAME
1142 : : !! poisson
1143 : : !!
1144 : : !! FUNCTION
1145 : : !! Solve poisson equation for angularly dependent charge
1146 : : !! distribution of angular momentum l
1147 : : !! Densities and potentials are given on a (generalized) radial grid
1148 : : !!
1149 : : !! INPUTS
1150 : : !! den(:)= electron density * (4*pi*r**2) appropriate for l
1151 : : !! ll= l quantum number
1152 : : !! radmesh <type(pawrad_type)>=data containing radial grid information
1153 : : !! [screened_sr_separation]= --optional-- separation for screened short-range kernel
1154 : : !! no screening by default
1155 : : !!
1156 : : !! OUTPUT
1157 : : !! [qq]= --optional-- lth moment of the charge; not compatible with screened Coulomb interaction
1158 : : !! rv(:)= electrostatic potential * r in (Hartree*Bohr) units
1159 : : !! where v(r)=\frac{1}{2l+1}(\frac{int[(r''^(l+2))g(r'')dr'']} {r^(l+1)}
1160 : : !! +(r^l) int[r''^(1-l)g(r'')dr''])
1161 : : !!
1162 : : !! SOURCE
1163 : :
1164 [ - + - + ]: 298693 : subroutine poisson(den,ll,radmesh,rv,screened_sr_separation,qq)
1165 : :
1166 : : !Arguments ---------------------------------------------
1167 : : !scalars
1168 : : integer,intent(in) :: ll
1169 : : real(dp),intent(out),optional :: qq
1170 : : real(dp),intent(in),optional :: screened_sr_separation
1171 : : type(pawrad_type),intent(in) :: radmesh
1172 : : !arrays
1173 : : real(dp),intent(in) :: den(:)
1174 : : real(dp),intent(out) :: rv(:)
1175 : :
1176 : : !Local variables ---------------------------------------
1177 : : !scalars
1178 : : integer :: ir,jr,mesh_size,mm,nn
1179 : : logical :: use_numerov,use_screened
1180 : : real(dp) :: angm,hh,intg,qq_,ri,rj,rr,sr_omega
1181 : : !arrays
1182 : : real(dp) :: ee(4)
1183 : 298693 : real(dp),allocatable :: aa(:),bb(:),cc(:),dd(:),radl(:),radl1(:)
1184 : :
1185 : : ! ***************************************************************************
1186 : :
1187 : 298693 : mesh_size=size(den)
1188 [ + - - + ]: 298693 : if (size(rv)/=mesh_size.or.mesh_size>radmesh%mesh_size) then
1189 : 0 : LIBPAW_BUG('wrong sizes!')
1190 : : end if
1191 : :
1192 : 298693 : use_numerov=(radmesh%mesh_type==1)
1193 : :
1194 : 298693 : use_screened=.false.;sr_omega=zero
1195 [ + + ]: 298693 : if (present(screened_sr_separation)) then
1196 : 358 : sr_omega=screened_sr_separation
1197 : 358 : use_screened=(sr_omega>tol8)
1198 : : ! Numerov method not coded for screened Coulomb interaction
1199 [ + - ]: 358 : if (use_numerov.and.use_screened) use_numerov=.false.
1200 : : end if
1201 : :
1202 : : !==============================================
1203 : : !UNIFORM GRID - NUMEROV ALGORITHM
1204 : : !==============================================
1205 [ + + ]: 298693 : if (use_numerov) then
1206 : 7449 : hh=radmesh%rstep
1207 : 7449 : nn=radmesh%int_meshsz-1
1208 [ - + + - : 22347 : LIBPAW_ALLOCATE(aa,(nn))
+ - + - -
+ - + ]
1209 [ - + + - : 14898 : LIBPAW_ALLOCATE(bb,(nn))
- + - + ]
1210 [ + - - + : 29796 : LIBPAW_ALLOCATE(cc,(nn+1))
+ - + - +
- - + -
+ ]
1211 [ + + ]: 5207879 : do ir=1,nn
1212 : 5200430 : aa(ir)=two*hh*den(ir+1)/(ir)
1213 : 5207879 : bb(ir)=den(ir+1)*((ir*hh)**ll)
1214 : : end do
1215 : 7449 : cc(1)=zero
1216 [ + + ]: 5207879 : cc(2:nn+1)=bb(1:nn)
1217 : 7449 : call simp_gen(qq_,cc,radmesh)
1218 : 7449 : qq_=qq_/dble(2*ll+1)
1219 : 7449 : rv(1)=aa(1)+0.1_dp*aa(2)
1220 [ + + ]: 5192981 : do ir=2,nn-1
1221 : 5192981 : rv(ir)=aa(ir)+0.1_dp*(aa(ir+1)+aa(ir-1))
1222 : : end do
1223 : 7449 : rv(nn)=aa(nn)+0.1_dp*aa(nn-1)
1224 : 7449 : angm=dble(ll*(ll+1))
1225 : 7449 : rr=(nn+1)*hh
1226 : 7449 : rv(nn)=rv(nn)+(2.4_dp-0.2_dp*angm/((nn+1)**2))*qq_/(rr**ll)
1227 [ + + ]: 5207879 : do ir=1,nn
1228 : 5200430 : aa(ir)=angm/(ir*ir)
1229 : 5207879 : bb(ir)=2.4_dp+aa(ir)
1230 : : end do
1231 [ + + ]: 5200430 : do ir=1,nn-1
1232 : 5200430 : cc(ir)=-1.2_dp+0.1_dp*aa(ir+1)
1233 : : end do
1234 [ + + ]: 5200430 : do ir=nn,2,-1
1235 : 5200430 : aa(ir)=-1.2_dp+0.1_dp*aa(ir-1)
1236 : : end do
1237 [ - + ]: 7449 : if (nn.eq.1) then
1238 : 0 : rv(2)=rv(1)/bb(1)
1239 : 0 : rv(1)=zero
1240 : : else
1241 [ + + ]: 5200430 : do ir=2,nn
1242 : 5200430 : bb(ir)=bb(ir)-aa(ir)*cc(ir-1)/bb(ir-1)
1243 : : end do
1244 : 7449 : rv(1)=rv(1)/bb(1)
1245 [ + + ]: 5200430 : do ir=2,nn
1246 : 5200430 : rv(ir)=(rv(ir)-aa(ir)*rv(ir-1))/bb(ir)
1247 : : end do
1248 [ + + ]: 5200430 : do ir=nn-1,1,-1
1249 : 5200430 : rv(ir)=rv(ir)-cc(ir)*rv(ir+1)/bb(ir)
1250 : : end do
1251 [ + + ]: 5207879 : do ir=nn+1,2,-1
1252 : 5207879 : rv(ir)=rv(ir-1)
1253 : : end do
1254 : 7449 : rv(1)=zero
1255 [ + + + + ]: 44544 : if (nn+1<mesh_size) rv(nn+2:mesh_size)=zero
1256 : : end if
1257 [ + + ]: 5252423 : rv(:)=half*rv(:)
1258 [ - + ]: 7449 : if (present(qq)) qq=qq_
1259 : 7449 : LIBPAW_DEALLOCATE(aa)
1260 [ - + ]: 7449 : LIBPAW_DEALLOCATE(bb)
1261 : 14898 : LIBPAW_DEALLOCATE(cc)
1262 : :
1263 : : ! ==============================================
1264 : : ! ANY OTHER GRID - SIMPSON ALGORITHM
1265 : : ! ==============================================
1266 : : else
1267 : 291244 : nn=mesh_size;hh=third*radmesh%stepint
1268 [ + + + + ]: 1487921 : do while (abs(den(nn))<tol16.and.nn>radmesh%int_meshsz)
1269 : 1201010 : nn=nn-1
1270 : : end do
1271 [ + + ]: 291244 : mm=nn;if (radmesh%mesh_type==3) mm=mm-1
1272 [ - + + - : 873732 : LIBPAW_ALLOCATE(aa,(nn))
+ - + - -
+ - + ]
1273 [ - + + - : 582488 : LIBPAW_ALLOCATE(bb,(nn))
- + - + ]
1274 [ - + + - : 582488 : LIBPAW_ALLOCATE(cc,(nn))
- + - + ]
1275 [ - + + - : 582488 : LIBPAW_ALLOCATE(dd,(nn))
- + - + ]
1276 : :
1277 [ + + ]: 291244 : if (.not.use_screened) then
1278 : :
1279 : : ! Standard Coulomb integral
1280 [ - + + - : 581772 : LIBPAW_ALLOCATE(radl,(nn))
- + - + ]
1281 [ - + + - : 581772 : LIBPAW_ALLOCATE(radl1,(nn))
- + - + ]
1282 [ + + ]: 281532282 : do jr=nn,2,-1
1283 : 281241396 : ir=nn-jr+1
1284 : 281241396 : radl(jr) =radmesh%rad(jr)**ll
1285 : 281241396 : radl1(jr)=radmesh%rad(jr)*radl(jr)
1286 : 281241396 : aa(ir)=den(jr)*radmesh%radfact(jr)*radl(jr)
1287 : 281532282 : bb(ir)=den(jr)*radmesh%radfact(jr)/radl1(jr)
1288 : : end do
1289 : 290886 : radl(1)=zero;radl1(1)=zero
1290 : 290886 : ee(2)=aa(nn-1);ee(3)=aa(nn-2);ee(4)=aa(nn-3)
1291 : 290886 : call pawrad_deducer0(ee,4,radmesh)
1292 : 290886 : aa(nn)=ee(1)
1293 : 290886 : ee(2)=bb(nn-1);ee(3)=bb(nn-2);ee(4)=bb(nn-3)
1294 : 290886 : call pawrad_deducer0(ee,4,radmesh)
1295 : 290886 : bb(nn)=ee(1)
1296 : 290886 : cc(1)=zero;dd(1)=zero
1297 [ - + ]: 290886 : do ir=3,mm,2
1298 : 140531037 : cc(ir) =cc(ir-2)+hh*(aa(ir-2)+four*aa(ir-1)+aa(ir))
1299 : 140531037 : cc(ir-1)=cc(ir-2)+hh*(1.25_dp*aa(ir-2)+two*aa(ir-1)-quarter*aa(ir))
1300 : 140531037 : dd(ir) =dd(ir-2)+hh*(bb(ir-2)+four*bb(ir-1)+bb(ir))
1301 [ + + ]: 140531037 : dd(ir-1)=dd(ir-2)+hh*(1.25_dp*bb(ir-2)+two*bb(ir-1)-quarter*bb(ir))
1302 : : end do
1303 [ + + ]: 290886 : if (mod(mm,2)==0) then
1304 : : ! cc(mm)=cc(mm-2)+hh*(aa(mm-2)+four*aa(mm-1)+aa(mm))
1305 : : ! dd(mm)=dd(mm-2)+hh*(bb(mm-2)+four*bb(mm-1)+bb(mm))
1306 : 168293 : cc(mm)=cc(mm-1)+hh*(1.25_dp*aa(mm-2)+two*aa(mm-1)-quarter*aa(mm))
1307 : 168293 : dd(mm)=dd(mm-1)+hh*(1.25_dp*bb(mm-2)+two*bb(mm-1)-quarter*bb(mm))
1308 : : end if
1309 [ + + ]: 290886 : if (mm<nn) then
1310 : 11029 : cc(nn)=cc(mm)+half*(aa(mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1))
1311 : 11029 : dd(nn)=dd(mm)+half*(bb(mm)+bb(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1))
1312 : : end if
1313 : 290886 : rv(1)=zero
1314 [ + + ]: 281532282 : do ir=2,nn
1315 : 281241396 : jr=nn-ir+1
1316 : 281532282 : rv(ir)=(dd(jr)*radl1(ir)+(cc(nn)-cc(jr))/radl(ir))/(two*ll+one)
1317 : : end do
1318 [ + + + + ]: 1485773 : if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn)
1319 [ - + ]: 290886 : if (present(qq)) qq=cc(nn)/(two*ll+one)
1320 : 290886 : LIBPAW_DEALLOCATE(radl)
1321 : 290886 : LIBPAW_DEALLOCATE(radl1)
1322 : : else
1323 : :
1324 : : ! Short-range screened Coulomb integral
1325 : 358 : rv(1)=zero
1326 [ + + ]: 278133 : do ir=2,nn
1327 : 277775 : ri=sr_omega*radmesh%rad(ir)
1328 [ + + ]: 264402876 : do jr=2,nn
1329 : 264125101 : rj=sr_omega*radmesh%rad(jr)
1330 : 264402876 : aa(jr)=den(jr)*radmesh%radfact(jr)*sr_omega*screened_coul_kernel(ll,ri,rj)
1331 : : end do
1332 : 277775 : call pawrad_deducer0(aa(1:4),4,radmesh)
1333 : 277775 : intg=zero
1334 : : ! Compute the integral in 2 parts because the function is not differentiable at r(ir)
1335 : : ! Integral from zero to r(ir)
1336 [ + + ]: 277775 : if (ir>2) then
1337 [ - + ]: 277417 : do jr=ir-2,1,-2
1338 [ + + ]: 66031246 : intg=intg+hh*(aa(jr)+four*aa(jr+1)+aa(jr+2))
1339 : : end do
1340 [ + + ]: 277417 : if (mod(ir,2)==0) intg=intg+hh*(-quarter*aa(1)+two*aa(2)+1.25_dp*aa(3))
1341 : : else if (ir==2) then
1342 [ + - ]: 358 : if (nn >2) intg=intg+hh*(+1.25_dp*aa(1)+two*aa(2)-quarter*aa(3))
1343 [ - + ]: 358 : if (nn==2) intg=intg+(aa(1)+aa(2))*hh*half
1344 : : end if
1345 [ - + ]: 277775 : if (mm<nn) intg=intg+half*(aa(1+nn-mm)+aa(nn))*(radmesh%rad(1+nn-mm)-radmesh%rad(1))
1346 : : ! Integral from r(ir) to rmax
1347 [ + + ]: 277775 : if (ir<nn-2) then
1348 [ - + ]: 276701 : do jr=ir+2,nn,2
1349 [ + + ]: 65892059 : intg=intg+hh*(aa(jr-2)+four*aa(jr-1)+aa(jr))
1350 : : end do
1351 [ + + ]: 276701 : if (mod(nn-ir,2)==1) then
1352 [ + + ]: 138471 : if (ir<=nn-4) then
1353 : : !Use a high-order formula because points are spaced
1354 : : intg=intg+hh*(251._dp*aa(nn)+646._dp*aa(nn-1)-264._dp*aa(nn-2) &
1355 : 138113 : & +106._dp*aa(nn-3)-19._dp*aa(nn-4))/240._dp
1356 : : else
1357 : 358 : intg=intg+hh*(1.25_dp*aa(nn)+two*aa(nn-1)-quarter*aa(nn-2))
1358 : : endif
1359 : : end if
1360 [ + + ]: 1074 : else if (ir==nn-1) then
1361 : 358 : intg=intg+(aa(nn-1)+aa(nn))*hh*half
1362 : : end if
1363 : 278133 : rv(ir)=intg/(two*ll+one)*radmesh%rad(ir)
1364 : : end do
1365 [ + - + + ]: 2148 : if (nn<mesh_size) rv(nn+1:mesh_size)=rv(nn)
1366 [ - + ]: 358 : if (present(qq)) qq=zero ! Not relevant
1367 : :
1368 : : end if
1369 [ - + ]: 291244 : LIBPAW_DEALLOCATE(aa)
1370 [ - + ]: 291244 : LIBPAW_DEALLOCATE(bb)
1371 [ - + ]: 291244 : LIBPAW_DEALLOCATE(cc)
1372 [ - + ]: 291244 : LIBPAW_DEALLOCATE(dd)
1373 : :
1374 : : end if
1375 : :
1376 [ - + ]: 298693 : end subroutine poisson
1377 : : !!***
1378 : :
1379 : : !----------------------------------------------------------------------
1380 : :
1381 : : !!****f* m_pawrad/pawrad_ifromr
1382 : : !! NAME
1383 : : !! pawrad_ifromr
1384 : : !!
1385 : : !! FUNCTION
1386 : : !! Retreive Index FROM a given R value in a radial grid
1387 : : !! Grid can be regular or logarithimc
1388 : : !!
1389 : : !! INPUTS
1390 : : !! rr=given input r value
1391 : : !! radmesh <type(pawrad_type)>=data containing radial grid information
1392 : : !!
1393 : : !! OUTPUT
1394 : : !! pawrad_ifromr=index of rr in radial grid
1395 : : !!
1396 : : !! NOTES
1397 : : !! Possible mesh types (radmesh%mesh_type)
1398 : : !! mesh_type=1 (regular grid): rad(i)=(i-1)*AA
1399 : : !! mesh_type=2 (logari. grid): rad(i)=AA*(exp[BB*(i-1)]-1)
1400 : : !! mesh_type=3 (logari. grid): rad(i>1)=AA*exp[BB*(i-2)] and rad(1)=0
1401 : : !! mesh_type=4 (logari. grid): rad(i)=-AA*ln[1-BB*(i-1)] with BB=1/n
1402 : : !! mesh_type=5 ( grid): rad(i)=AA*i/(n-i)
1403 : : !!
1404 : : !! SOURCE
1405 : :
1406 : 11341 : function pawrad_ifromr(radmesh,rr)
1407 : :
1408 : : !Arguments ------------------------------------
1409 : : !scalars
1410 : : integer :: pawrad_ifromr
1411 : : real(dp),intent(in) :: rr
1412 : : type(pawrad_type),intent(in) :: radmesh
1413 : :
1414 : : !Local variables-------------------------------
1415 : : character(len=500) :: msg
1416 : :
1417 : : ! *************************************************************************
1418 : :
1419 [ + + + - : 11341 : if (radmesh%mesh_type==1) then
+ - ]
1420 : 158 : pawrad_ifromr=int(tol8+rr/radmesh%rstep)+1
1421 : : else if (radmesh%mesh_type==2) then
1422 : 10659 : pawrad_ifromr=int(tol8+log(1.d0+rr/radmesh%rstep)/radmesh%lstep)+1
1423 : : else if (radmesh%mesh_type==3) then
1424 [ + - ]: 484 : if (rr<radmesh%rstep) then
1425 : : pawrad_ifromr=1
1426 : : else
1427 : 484 : pawrad_ifromr=int(tol8+log(rr/radmesh%rstep)/radmesh%lstep)+2
1428 : : end if
1429 : : else if (radmesh%mesh_type==4) then
1430 : 0 : pawrad_ifromr=int(tol8+(1.d0-exp(-rr/radmesh%rstep))/radmesh%lstep)+1
1431 : : else if (radmesh%mesh_type==5) then
1432 : 40 : pawrad_ifromr=int(tol8+(radmesh%lstep*rr)/(radmesh%rstep+rr))+1
1433 : : else
1434 : : ! Other values of mesh_type are not allowed (see psp7in.F90)
1435 : 0 : write(msg,'(a,i0)')" Unknown value of %mesh_type ",radmesh%mesh_type
1436 : 0 : LIBPAW_ERROR(msg)
1437 : : end if
1438 : :
1439 : 11341 : end function pawrad_ifromr
1440 : : !!***
1441 : :
1442 : : !----------------------------------------------------------------------
1443 : :
1444 : : !!****f* m_pawrad/screened_coul_kernel
1445 : : !! NAME
1446 : : !! screened_coul_kernel
1447 : : !!
1448 : : !! FUNCTION
1449 : : !! Evaluates the kernel function used to compute the integral
1450 : : !! of a screened Coulomb potential (with erfc function)
1451 : : !! See Angyan, Gerber, Marsman, J. Phys. A: Math. Gen. 39, 8613 (2006) [[cite:Angyan2006]]
1452 : : !!
1453 : : !! INPUTS
1454 : : !! [formula]=optional; used to force formula (1=full; 2=dev. near r2=0; 3=dev near r1*r2=0)
1455 : : !! order= order of the function (typically l quantum number)
1456 : : !! r1,r2=input arguments
1457 : : !!
1458 : : !! OUTPUT
1459 : : !! screened_coul_kernel=output radial function
1460 : : !!
1461 : : !! SOURCE
1462 : :
1463 : 264125101 : function screened_coul_kernel(order,r1,r2,formula)
1464 : :
1465 : : !Arguments ------------------------------------
1466 : : !scalars
1467 : : integer,intent(in) :: order
1468 : : integer,optional :: formula
1469 : : real(dp),intent(in) :: r1,r2
1470 : : real(dp) :: screened_coul_kernel
1471 : : !Local variables-------------------------------
1472 : : integer :: formula_
1473 : : real(dp) :: difexp,erfcp,erfcm,f0,f1,f2,f3,f4,f5,f6,hh,sqrtpi,sumexp,xx,xy,yy
1474 : : character(len=100) :: msg
1475 : :
1476 : : !******************************************************************
1477 : :
1478 [ - + ]: 264125101 : if (order>6) then
1479 : 0 : msg='PAW screened exchange not coded for l>2!'
1480 : 0 : LIBPAW_ERROR(msg)
1481 : : end if
1482 : :
1483 : : !Use max and min of arguments
1484 : 264125101 : xx=max(r1,r2) ; yy=min(r1,r2)
1485 : :
1486 : : !Unscreened Coulomb interaction for very small (or negative) arguments
1487 [ - + ]: 264125101 : if (xx<tol8) then
1488 : 0 : screened_coul_kernel=yy**order/xx**(order+1) ; return
1489 : : end if
1490 : :
1491 : : !Choice of formula
1492 : : !Empirical criterion, inspired by J. Phys. A 39, pp8624 [[cite:Angyan2006]] and adjusted
1493 [ - + - + ]: 264125101 : formula_=1;if (xx<0.25_dp.and.yy<0.25_dp) formula_=2
1494 [ - + ]: 264125101 : if (present(formula)) formula_=formula
1495 : 0 : select case (formula_)
1496 : :
1497 : : !Full formula
1498 : : !J. Phys. A 39, 8613 (2006) - Eq. (21), (22), (24) [[cite:Angyan2006]]
1499 : : ! Note: typo in the paper: in eq (22), the sum begins at p=0
1500 : : !------------------------------------------------------------------
1501 : : case(1)
1502 : 0 : xy=xx*yy ; sqrtpi=sqrt(pi)
1503 : 0 : sumexp=exp(-(xx+yy)**2)+exp(-(xx-yy)**2)
1504 : 0 : difexp=exp(-(xx+yy)**2)-exp(-(xx-yy)**2)
1505 [ # # # # : 0 : erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy)
# # # # ]
1506 : 264125101 : if (order==0) then
1507 : 0 : f0=-difexp/(two*sqrtpi*xy)
1508 : 0 : hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy)
1509 : 0 : screened_coul_kernel = hh + f0
1510 : : else if(order==1) then
1511 : 0 : f0=-difexp/(two*sqrtpi*xy)
1512 : 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1513 : 0 : hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2)
1514 : 0 : screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy
1515 : : else if (order==2) then
1516 : 0 : f0=-difexp/(two*sqrtpi*xy)
1517 : 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1518 : 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1519 : 0 : hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3)
1520 : 0 : screened_coul_kernel = hh + f2 + f1*(xx**2+yy**2)/xy + f0*(xx**4+yy**4)/xy**2
1521 : : else if (order==3) then
1522 : 0 : f0=-difexp/(two*sqrtpi*xy)
1523 : 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1524 : 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1525 : 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1526 : 0 : hh=((xx**7+yy**7)*erfcp-(xx**7-yy**7)*erfcm)/(two*xy**4)
1527 : : screened_coul_kernel = hh + f3 + f2*(xx**2+yy**2)/xy + f1*(xx**4+yy**4)/xy**2 &
1528 : 0 : & + f0*(xx**6+yy**6)/xy**3
1529 : : else if (order==4) then
1530 : 0 : f0=-difexp/(two*sqrtpi*xy)
1531 : 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1532 : 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1533 : 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1534 : : f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp &
1535 : 0 : & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5)
1536 : 0 : hh=((xx**9+yy**9)*erfcp-(xx**9-yy**9)*erfcm)/(two*xy**5)
1537 : : screened_coul_kernel = hh + f4 + f3*(xx**2+yy**2)/xy + f2*(xx**4+yy**4)/xy**2 &
1538 : 0 : & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4
1539 : : else if (order==5) then
1540 : 0 : f0=-difexp/(two*sqrtpi*xy)
1541 : 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1542 : 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1543 : 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1544 : : f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp &
1545 : 0 : & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5)
1546 : : f5=((240._dp*xy**4+1680._dp*xy**2+945._dp)*difexp &
1547 : 0 : & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6)
1548 : 0 : hh=((xx**11+yy**11)*erfcp-(xx**11-yy**11)*erfcm)/(two*xy**6)
1549 : : screened_coul_kernel = hh + f5 + f4*(xx**2+yy**2) /xy + f3*(xx**4+yy**4)/xy**2 &
1550 : : & + f2*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8)/xy**4 &
1551 : 0 : & + f0*(xx**10+yy**10)/xy**5
1552 : : else if (order==6) then
1553 : 0 : f0=-difexp/(two*sqrtpi*xy)
1554 : 0 : f1=(difexp+two*xy*sumexp)/(sqrtpi*(two*xy)**2)
1555 : 0 : f2= -((four*xy**2+three)*difexp+6._dp*xy*sumexp)/(sqrtpi*(two*xy)**3)
1556 : 0 : f3=((24._dp*xy**2+15._dp)*difexp+(8._dp*xy**3+30._dp*xy)*sumexp)/(sqrtpi*(two*xy)**4)
1557 : : f4=-((16._dp*xy**4+180._dp*xy**2+105._dp)*difexp &
1558 : 0 : & +(80._dp*xy**3+210._dp*xy)*sumexp)/(sqrtpi*(two*xy)**5)
1559 : : f5=((240._dp*xy**4+1680._dp*xy**2+945._dp)*difexp &
1560 : 0 : & +(32._dp*xy**5+840._dp*xy**3+1890._dp*xy)*sumexp)/(sqrtpi*(two*xy)**6)
1561 : : f6=-((64._dp*xy**6+3360._dp*xy**4+18900._dp*xy**2+10395._dp)*difexp+ &
1562 : 0 : & (672._dp*xy**5+10080._dp*xy**3+20790._dp*xy)*sumexp)/(sqrtpi*(two*xy)**7)
1563 : 0 : hh=((xx**13+yy**13)*erfcp-(xx**13-yy**13)*erfcm)/(two*xy**7)
1564 : : screened_coul_kernel = hh + f6 + f5*(xx**2+yy**2) /xy + f4*(xx**4+yy**4) /xy**2 &
1565 : : & + f3*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8) /xy**4 &
1566 : 0 : & + f1*(xx**10+yy**10)/xy**5 + f0*(xx**12+yy**12)/xy**6
1567 : : end if
1568 : :
1569 : : !Development for yy->0
1570 : : !J. Phys. A 39, 8613 (2006) - Eq. (28), (29), (30) [[cite:Angyan2006]]
1571 : : !------------------------------------------------------------------
1572 : : case(2)
1573 [ - + + + : 264125101 : sqrtpi=sqrt(pi)
- - - - ]
1574 : 0 : if (order==0) then
1575 : : screened_coul_kernel = paw_derfc(xx)/xx + exp(-xx**2)/sqrtpi * &
1576 : : & (two/three *yy**2 &
1577 : 113142729 : & +(two*xx**2-three)/15._dp*yy**4)
1578 : : else if(order==1) then
1579 : : screened_coul_kernel = paw_derfc(xx)*yy/xx**2 + exp(-xx**2)/sqrtpi * &
1580 : : & (two/xx *yy &
1581 : : & +four*xx/5._dp *yy**3 &
1582 : 55508225 : & +two*xx*(two*xx**2-5._dp)/35._dp*yy**5)
1583 : : else if (order==2) then
1584 : : screened_coul_kernel = paw_derfc(xx)*yy**2/xx**3 + exp(-xx**2)/sqrtpi * &
1585 : : & (two*(two*xx**2+three)/(three*xx**2) *yy**2 &
1586 : : & +8._dp*xx**2/21._dp *yy**4 &
1587 : 95474147 : & +four*xx**2*(two*xx**2-7._dp)/189._dp*yy**6)
1588 : : else if (order==3) then
1589 : : screened_coul_kernel = paw_derfc(xx)*yy**3/xx**4 + exp(-xx**2)/sqrtpi * &
1590 : : & (two*(four*xx**4+10._dp*xx**2+15._dp)/(15._dp*xx**3)*yy**3 &
1591 : : & +16._dp*xx**3/135._dp *yy**5 &
1592 : 0 : & +8._dp*xx**3*(two*xx**2-9._dp)/1485._dp *yy**7)
1593 : : else if (order==4) then
1594 : : screened_coul_kernel = paw_derfc(xx)*yy**4/xx**5 + exp(-xx**2)/sqrtpi * &
1595 : : & (two*(8._dp*xx**6+28._dp*xx**4+70._dp*xx**2+105._dp)/(105._dp*xx**4)*yy**4 &
1596 : : & +32._dp*xx**4/1155._dp *yy**6 &
1597 : 0 : & +16._dp*xx**4*(two*xx**2-11._dp)/15015._dp *yy**8)
1598 : : else if (order==5) then
1599 : : screened_coul_kernel = paw_derfc(xx)*yy**5/xx**6 + exp(-xx**2)/sqrtpi * &
1600 : : & (two*(16._dp*xx**8+72._dp*xx**6+252._dp*xx**4+630._dp*xx**2+945._dp)/(945._dp*xx**5)*yy**5 &
1601 : : & +64._dp*xx**5/12285._dp *yy**7 &
1602 : 0 : & +32._dp*xx**5*(two*xx**2-13._dp)/184275._dp *yy**9)
1603 : : else if (order==6) then
1604 : : screened_coul_kernel = paw_derfc(xx)*yy**6/xx**7 + exp(-xx**2)/sqrtpi * &
1605 : : & (two*(32._dp*xx**10+176._dp*xx**8+792._dp*xx**6+2772._dp*xx**4+6930._dp*xx**2+10395._dp)/(10395._dp*xx**6)*yy**6 &
1606 : : & +128._dp*xx**6/155925._dp *yy**8 &
1607 : 0 : & +64._dp*xx**6*(two*xx**2-15._dp)/2650725._dp *yy**10)
1608 : : end if
1609 : :
1610 : : !Development for xx*yy->0
1611 : : !J. Phys. A 39, 8613 (2006) - Eq. (24), (26) [[cite:Angyan2006]]
1612 : : ! Note: typo in the paper: (2n+3)! should be (2n+3)!!
1613 : : !------------------------------------------------------------------
1614 : : case(3)
1615 : 0 : xy=xx*yy ; sqrtpi=sqrt(pi)
1616 : 0 : sumexp=exp(-(xx**2+yy**2))
1617 [ # # # # : 0 : erfcp=paw_derfc(xx+yy) ; erfcm=paw_derfc(xx-yy)
# # # # ]
1618 [ - - + - ]: 264125101 : if (order==0) then
1619 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1620 : 0 : hh=((xx+yy)*erfcp-(xx-yy)*erfcm)/(two*xy)
1621 : 0 : screened_coul_kernel = hh + f0
1622 : : else if(order==1) then
1623 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1624 : 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1625 : 0 : hh=((xx**3+yy**3)*erfcp-(xx**3-yy**3)*erfcm)/(two*xy**2)
1626 : 0 : screened_coul_kernel = hh + f1 + f0*(xx**2+yy**2)/xy
1627 : : else if (order==2) then
1628 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1629 : 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1630 : 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1631 : 0 : hh=((xx**5+yy**5)*erfcp-(xx**5-yy**5)*erfcm)/(two*xy**3)
1632 : 0 : screened_coul_kernel = hh + f2 + f1*(xx**2+yy**2)/xy + f0*(xx**4+yy**4)/xy**2
1633 : : else if (order==3) then
1634 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1635 : 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1636 : 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1637 : 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1638 : 0 : hh=((xx**7+yy**7)*erfcp-(xx**7-yy**7)*erfcm)/(two*xy**4)
1639 : : screened_coul_kernel = hh + f3 + f2*(xx**2+yy**2)/xy + f1*(xx**4+yy**4)/xy**2 &
1640 : 0 : & + f0*(xx**6+yy**6)/xy**3
1641 : : else if (order==4) then
1642 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1643 : 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1644 : 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1645 : 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1646 : 0 : f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp
1647 : 0 : hh=((xx**9+yy**9)*erfcp-(xx**9-yy**9)*erfcm)/(two*xy**5)
1648 : : screened_coul_kernel = hh + f4 + f3*(xx**2+yy**2)/xy + f2*(xx**4+yy**4)/xy**2 &
1649 : 0 : & + f1*(xx**6+yy**6)/xy**3 + f0*(xx**8+yy**8)/xy**4
1650 : : else if (order==5) then
1651 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1652 : 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1653 : 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1654 : 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1655 : 0 : f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp
1656 : 0 : f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp
1657 : 0 : hh=((xx**11+yy**11)*erfcp-(xx**11-yy**11)*erfcm)/(two*xy**6)
1658 : : screened_coul_kernel = hh + f5 + f4*(xx**2+yy**2) /xy + f3*(xx**4+yy**4)/xy**2 &
1659 : : & + f2*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8)/xy**4 &
1660 : 0 : & + f0*(xx**10+yy**10)/xy**5
1661 : : else if (order==6) then
1662 : 0 : f0 = two *(three +two*xy**2)/(three *sqrtpi)*sumexp
1663 : 0 : f1 = four *xy *(5._dp +two*xy**2)/(15._dp *sqrtpi)*sumexp
1664 : 0 : f2 = 8._dp *xy**2*(7._dp +two*xy**2)/(105._dp *sqrtpi)*sumexp
1665 : 0 : f3 = 16._dp *xy**2*(9._dp +two*xy**2)/(945._dp *sqrtpi)*sumexp
1666 : 0 : f4 = 32._dp *xy**2*(11._dp+two*xy**2)/(10395._dp *sqrtpi)*sumexp
1667 : 0 : f5 = 64._dp *xy**2*(13._dp+two*xy**2)/(135135._dp *sqrtpi)*sumexp
1668 : 0 : f6 = 128._dp*xy**2*(15._dp+two*xy**2)/(1027025._dp*sqrtpi)*sumexp
1669 : 0 : hh=((xx**13+yy**13)*erfcp-(xx**13-yy**13)*erfcm)/(two*xy**7)
1670 : : screened_coul_kernel = hh + f6 + f5*(xx**2+yy**2) /xy + f4*(xx**4+yy**4) /xy**2 &
1671 : : & + f3*(xx**6+yy**6) /xy**3 + f1*(xx**8+yy**8) /xy**4 &
1672 : 0 : & + f1*(xx**10+yy**10)/xy**5 + f0*(xx**12+yy**12)/xy**6
1673 : : end if
1674 : :
1675 : : end select
1676 : :
1677 : : end function screened_coul_kernel
1678 : : !!***
1679 : :
1680 : : !----------------------------------------------------------------------
1681 : :
1682 : : !!****f* m_pawrad/calc_slatradl
1683 : : !! NAME
1684 : : !! calc_slatradl
1685 : : !!
1686 : : !! FUNCTION
1687 : : !! Calculate the radial part of Slater integrals. See below.
1688 : : !!
1689 : : !! INPUTS
1690 : : !! ll= l quantum number in the expansion of the Coulomb term.
1691 : : !! mesh_size=Number of points on the radial mesh.
1692 : : !! ff1(radmesh), ff2(radmesh)= The two functions to be integrated.
1693 : : !! Pawrad <type(pawrad_type)>=Structure containing radial grid information.
1694 : : !!
1695 : : !! OUTPUT
1696 : : !! integral =
1697 : : !! $ \dfrac{4\pi}{2L+1} \int ff1(r1) \dfrac{r_<^L}{r_>^{L+1}} ff2(r2) dr1 dr2 $
1698 : : !! where $r_< = min(r1,r2)$ and $r_> = Max(r1,r2)$.
1699 : : !!
1700 : : !! SOURCE
1701 : :
1702 : 162 : subroutine calc_slatradl(ll,mesh_size,ff1,ff2,Pawrad,integral)
1703 : :
1704 : : !scalars
1705 : : integer,intent(in) :: mesh_size,ll
1706 : : real(dp),intent(out) :: integral
1707 : : !arrays
1708 : : real(dp),intent(in) :: ff1(mesh_size),ff2(mesh_size)
1709 : : type(pawrad_type),intent(in) :: Pawrad
1710 : :
1711 : : !Local variables ---------------------------------------
1712 : : !scalars
1713 : : integer :: int_meshsz
1714 : : character(len=100) :: msg
1715 : : !arrays
1716 : 162 : real(dp),allocatable :: hh(:),gg(:)
1717 : :
1718 : : ! *************************************************************************
1719 : :
1720 [ - + ]: 162 : if (mesh_size > Pawrad%mesh_size) then
1721 : 0 : msg='mesh_size > pawrad%mesh_size!'
1722 : 0 : LIBPAW_BUG(msg)
1723 : : end if
1724 : :
1725 [ - + + - : 486 : LIBPAW_ALLOCATE(hh,(mesh_size))
+ - + - -
+ - + ]
1726 [ - + + - : 324 : LIBPAW_ALLOCATE(gg,(mesh_size))
- + - + ]
1727 : : !hh = zero
1728 : : !gg = zero
1729 : :
1730 : 162 : int_meshsz=Pawrad%int_meshsz
1731 : : !$int_meshsz=Pawrad%mesh_size
1732 : :
1733 : : ! the line below requires hh as work array.
1734 [ - + + + ]: 224208 : hh = ff2
1735 : :
1736 : : ! TODO find where int_meshsz is calculated and if it can affects the results.
1737 [ - + - - ]: 162 : if (int_meshsz<mesh_size) hh(int_meshsz+1:mesh_size)=zero
1738 : :
1739 : 162 : call poisson(hh,ll,Pawrad,gg)
1740 : :
1741 [ + + ]: 223884 : gg(2:mesh_size) = gg(2:mesh_size)/Pawrad%rad(2:mesh_size)
1742 : :
1743 : 162 : hh(1) = zero
1744 [ + + ]: 223884 : hh(2:mesh_size)= ff1(2:mesh_size) * gg(2:mesh_size)
1745 [ - + ]: 162 : LIBPAW_DEALLOCATE(gg)
1746 : :
1747 : 162 : call simp_gen(integral,hh,Pawrad)
1748 : 162 : integral = four_pi * integral
1749 : :
1750 : 162 : LIBPAW_DEALLOCATE(hh)
1751 : :
1752 : 162 : end subroutine calc_slatradl
1753 : : !!***
1754 : :
1755 : : !----------------------------------------------------------------------
1756 : :
1757 [ # # # # : 0 : END MODULE m_pawrad
# # # # #
# # # # #
# # # # ]
1758 : : !!***
|
