Octopus
atomic.F90
Go to the documentation of this file.
1!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
21module atomic_oct_m
22 use accel_oct_m
23 use debug_oct_m
24 use global_oct_m
25 use, intrinsic :: iso_fortran_env
26 use logrid_oct_m
27 use messages_oct_m
28 use namespace_oct_m
29 use profiling_oct_m
30 use sort_oct_m
31 use types_oct_m
32 use xc_f03_lib_m
33#ifdef HAVE_LIBXC_FUNCS
34 use xc_f03_funcs_m
35#endif
36 use xc_oct_m, only: xc_is_using_device
37 use xc_functional_oct_m, only: xc_functional_init_func
38 use iso_c_binding
39 implicit none
40
41 private
42 public :: &
43 atomxc, &
44 atomhxc, &
45 vhrtre, egofv
46
48 public :: &
49 valconf_t, &
50 valconf_copy, &
51 valconf_sort_by_l, &
52 valconf_unpolarized_to_polarized, &
53 valconf_string_length
54
55 integer, parameter :: VALCONF_STRING_LENGTH = 80
56
57 type valconf_t
58 ! Components are public by default
59 integer :: z = 0
60 character(len=3) :: symbol = ""
61 integer :: type = 0
62 integer :: p = 0
63 integer :: n(15) = 0
64 integer :: l(15) = 0
65 real(real64) :: occ(15,2) = m_zero
66 real(real64) :: j(15) = m_zero
67 end type valconf_t
68
69
70contains
71
72 ! ---------------------------------------------------------
73 subroutine valconf_copy(cout, cin)
74 type(valconf_t), intent(out) :: cout
75 type(valconf_t), intent(in) :: cin
76
77 push_sub(valconf_copy)
78
79 cout%z = cin%z
80 cout%symbol = cin%symbol
81 cout%type = cin%type
82 cout%p = cin%p
83 cout%n = cin%n
84 cout%l = cin%l
85 cout%occ = cin%occ
86 cout%j = cin%j
87
88 pop_sub(valconf_copy)
89 end subroutine valconf_copy
90
91 subroutine valconf_sort_by_l(c, lmax)
92 type(valconf_t), intent(inout) :: c
93 integer, intent(in) :: lmax
94
95 integer :: ll
96 integer, allocatable :: order(:)
97 type(valconf_t) :: tmp_c
98
99 push_sub(valconf_sort_by_l)
100
101 call valconf_copy(tmp_c, c)
102
103 ! In order to work in the following, we need to sort the occupations by angular momentum
104 safe_allocate(order(1:lmax+1))
105 call sort(c%l(1:lmax+1), order)
106 do ll = 0, lmax
107 c%occ(ll+1, 1:2) = tmp_c%occ(order(ll+1), 1:2)
108 c%n(ll+1) = tmp_c%n(order(ll+1))
109 c%j(ll+1) = tmp_c%j(order(ll+1))
110 end do
111 safe_deallocate_a(order)
112
113 pop_sub(valconf_sort_by_l)
114 end subroutine valconf_sort_by_l
115
116 !In case of spin-polarized calculations, we properly distribute the electrons
117 subroutine valconf_unpolarized_to_polarized(conf)
118 type(valconf_t), intent(inout) :: conf
119
120 integer :: ii, ll
121 real(real64) :: xx
122
123 push_sub(valconf_unpolarized_to_polarized)
124
125 do ii = 1, conf%p
126 ll = conf%l(ii)
127 xx = conf%occ(ii, 1)
128 conf%occ(ii, 1) = min(xx, real(2*ll+1, real64) )
129 conf%occ(ii, 2) = xx - conf%occ(ii,1)
130 end do
131
132 pop_sub(valconf_unpolarized_to_polarized)
133 end subroutine valconf_unpolarized_to_polarized
134
135
136 ! ---------------------------------------------------------
137 subroutine atomhxc(namespace, functl, g, nspin, dens, v, extra)
138 type(namespace_t), intent(in) :: namespace
139 character(len=*), intent(in) :: functl
140 type(logrid_t), intent(in) :: g
141 integer, intent(in) :: nspin
142 real(real64), intent(in) :: dens(g%nrval, nspin)
143 real(real64), intent(out) :: v(g%nrval, nspin)
144 real(real64), intent(in), optional :: extra(g%nrval)
145
146 character(len=5) :: xcfunc, xcauth
147 integer :: is, ir
148 real(real64), allocatable :: xc(:, :), ve(:, :), rho(:, :)
149 real(real64) :: r2, ex, ec, dx, dc
150
151 push_sub(atomhxc)
152
153 safe_allocate( ve(1:g%nrval, 1:nspin))
154 safe_allocate( xc(1:g%nrval, 1:nspin))
155 safe_allocate(rho(1:g%nrval, 1:nspin))
156 ve = m_zero
157 xc = m_zero
158 rho = m_zero
159
160 ! To calculate the Hartree term, we put all the density in one variable.
161 do is = 1, nspin
162 rho(:, 1) = rho(:, 1) + dens(:, is)
163 end do
164 ve = m_zero
165 do is = 1, nspin
166 call vhrtre(rho(:, 1), ve(:, is), g%rofi, g%drdi, g%s, g%nrval, g%a)
167 end do
168 ve(1, 1:nspin) = ve(2, 1:nspin)
169
170 select case (functl)
171 case ('LDA')
172 xcfunc = 'LDA'
173 xcauth = 'PZ'
174 case ('GGA')
175 xcfunc = 'GGA'
176 xcauth = 'PBE'
177 case default
178 message(1) = 'Internal Error in atomhxc: unknown functl'
179 call messages_fatal(1, namespace=namespace)
180 end select
181
182 rho(:, :) = dens(:, :)
183 do is = 1, nspin
184 if (present(extra)) rho(:, is) = rho(:, is) + extra(:)/nspin
185 rho(2:g%nrval, is) = rho(2:g%nrval, is)/(m_four*m_pi*g%r2ofi(2:g%nrval))
186 end do
187 r2 = g%rofi(2)/(g%rofi(3)-g%rofi(2))
188 rho(1, 1:nspin) = rho(2, 1:nspin) - (rho(3, 1:nspin)-rho(2, 1:nspin))*r2
189
190 do is = 1, nspin
191 do ir = 1, g%nrval
192 if (rho(ir, is) < m_epsilon) rho(ir, is) = m_zero
193 end do
194 end do
195
196 call atomxc(namespace, xcfunc, xcauth, g%nrval, g%nrval, g%rofi, nspin, rho, ex, ec, dx, dc, xc)
197 v = ve + xc
198
199 safe_deallocate_a(ve)
200 safe_deallocate_a(xc)
201 safe_deallocate_a(rho)
202
203 pop_sub(atomhxc)
204 end subroutine atomhxc
205
206
207!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
246 subroutine atomxc(namespace, FUNCTL, AUTHOR, NR, MAXR, RMESH, NSPIN, DENS, EX, EC, DX, DC, V_XC)
247 type(namespace_t), intent(in) :: namespace
248 character(len=*), intent(in) :: functl, author
249 integer, intent(in) :: nr, maxr
250 real(real64), intent(in) :: rmesh(maxr)
251 integer, intent(in) :: nspin
252 real(real64), intent(in) :: dens(maxr,nspin)
253 real(real64), intent(out) :: dc, dx, ec, ex, v_xc(maxr,nspin)
254
255!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
256! Internal parameters !
257! NN : order of the numerical derivatives: the number of radial !
258! points used is 2*NN+1 !
259!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
260
261 integer, parameter :: nn = 5
262
263!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
264! Fix energy unit: EUNIT=1.0 => Hartrees, !
265! EUNIT=0.5 => Rydbergs, !
266! EUNIT=0.03674903 => eV !
267!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
268
269 real(real64), parameter :: eunit = m_half
270
271!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
272! DVMIN is added to differential of volume to avoid division by zero !
273!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
274
275 real(real64), parameter :: dvmin = 1.0e-12_real64
276
277!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
278! Local variables and arrays !
279!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
280
281 logical :: gga
282 integer :: in, in1, in2, ir, is, jn
283 real(real64) :: aux(maxr), dgdm(-nn:nn), dgidfj(-nn:nn), drdm, dvol, f1, f2, gd(3,nspin), &
284 dexdgd(3,nspin), decdgd(3,nspin)
285 real(real64), target :: d(nspin), sigma(3), vxsigma(3), vcsigma(3), &
286 epsx(1), epsc(1), dexdd(nspin), decdd(nspin)
287 type(xc_f03_func_t) :: x_conf, c_conf
288
289 logical :: xc_on_host
290 type(accel_mem_t) :: d_buff, decdd_buff, dexdd_buff, sigma_buff, vxsigma_buff, vcsigma_buff, epsx_buff, epsc_buff
291 real(c_double), pointer :: p_d(:), p_sigma(:), p_epsx(:), p_epsc(:), p_dexdd(:), p_decdd(:), p_vxsigma(:), p_vcsigma(:)
292
293 call profiling_in("ATOMXC")
294 push_sub(atomxc)
295
296 xc_on_host = .not. xc_is_using_device(namespace)
297
298 ! sanity check
299 assert(nspin == 1 .or. nspin == 2)
300
301 ! Set GGA switch
302 if (functl .eq. 'LDA' .or. functl .eq. 'lda' .or. &
303 functl .eq. 'LSD' .or. functl .eq. 'lsd') then
304 gga = .false.
305 else if (functl .eq. 'GGA' .or. functl .eq. 'gga') then
306 gga = .true.
307 else
308 gga = .false.
309 write(message(1),'(a,a)') 'atomxc: Unknown functional ', functl
310 call messages_fatal(1, namespace=namespace)
311 end if
312
313 ! initialize xc functional
314 if (gga) then
315 call xc_functional_init_func(x_conf, xc_gga_x_pbe, nspin, xc_on_host)
316 call xc_functional_init_func(c_conf, xc_gga_c_pbe, nspin, xc_on_host)
317 else
318 call xc_functional_init_func(x_conf, xc_lda_x, nspin, xc_on_host)
319 if (author .eq. 'CA' .or. author .eq. 'ca' .or. &
320 author .eq. 'PZ' .or. author .eq. 'pz') then
321 call xc_functional_init_func(c_conf, xc_lda_c_pz, nspin, xc_on_host)
322 else if (author .eq. 'PW92' .or. author .eq. 'pw92') then
323 call xc_functional_init_func(c_conf, xc_lda_c_pw, nspin, xc_on_host)
324 else
325 write(message(1),'(a,a)') 'LDAXC: Unknown author ', author
326 call messages_fatal(1, namespace=namespace)
327 end if
328 end if
329
330 if (.not. xc_on_host) then
331 call accel_create_buffer(d_buff, accel_mem_read_write, type_float, nspin)
332 call accel_create_buffer(dexdd_buff, accel_mem_read_write, type_float, nspin)
333 call accel_create_buffer(decdd_buff, accel_mem_read_write, type_float, nspin)
334 call accel_create_buffer(epsx_buff, accel_mem_read_write, type_float, 1)
335 call accel_create_buffer(epsc_buff, accel_mem_read_write, type_float, 1)
336
337 if (gga) then
338 call accel_create_buffer(sigma_buff, accel_mem_read_write, type_float, 3)
339 call accel_create_buffer(vxsigma_buff, accel_mem_read_write, type_float, 3)
340 call accel_create_buffer(vcsigma_buff, accel_mem_read_write, type_float, 3)
341 end if
342 end if
343
344!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
345! Initialize output !
346!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
347
348 ex = m_zero
349 ec = m_zero
350 dx = m_zero
351 dc = m_zero
352 do is = 1,nspin
353 do ir = 1,nr
354 v_xc(ir,is) = m_zero
355 end do
356 end do
357
358!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
359! Loop on mesh points !
360!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
361
362 do ir = 1,nr
363
364 ! Find interval of neighbour points to calculate derivatives
365 in1 = max( 1, ir-nn ) - ir
366 in2 = min( nr, ir+nn ) - ir
367
368 ! Find weights of numerical derivative from Lagrange
369 ! interpolation formula
370 do in = in1,in2
371 if (in == 0) then
372 dgdm(in) = 0
373 do jn = in1,in2
374 if (jn /= 0) dgdm(in) = dgdm(in) + m_one / (0 - jn)
375 end do
376 else
377 f1 = 1
378 f2 = 1
379 do jn = in1,in2
380 if (jn /= in .and. jn /= 0) f1 = f1 * (0 - jn)
381 if (jn /= in) f2 = f2 * (in - jn)
382 end do
383 dgdm(in) = f1 / f2
384 end if
385 end do
386
387 ! Find dr/dmesh
388 drdm = 0
389 do in = in1,in2
390 drdm = drdm + rmesh(ir+in) * dgdm(in)
391 end do
392
393 ! Find differential of volume. Use trapezoidal integration rule
394 dvol = 4 * m_pi * rmesh(ir)**2 * drdm
395
396 ! DVMIN is a small number added to avoid a division by zero
397 dvol = safe_tol(dvol, dvmin)
398 if (ir == 1 .or. ir == nr) dvol = dvol / 2
399 if (gga) aux(ir) = dvol
400
401 ! Find the weights for the derivative d(gradF(i))/d(F(j)), of
402 ! the gradient at point i with respect to the value at point j
403 if (gga) then
404 do in = in1,in2
405 dgidfj(in) = dgdm(in) / drdm
406 end do
407 end if
408
409 ! Find density and gradient of density at this point
410 do is = 1,nspin
411 d(is) = dens(ir,is)
412 end do
413 if (gga) then
414 do is = 1,nspin
415 gd(:,is) = m_zero
416 do in = in1,in2
417 gd(3,is) = gd(3,is) + dgidfj(in) * dens(ir+in,is)
418 end do
419 end do
420 else
421 gd = m_zero
422 end if
423
424 ! The xc_f03_gga routines need as input these combinations of
425 ! the gradient of the densities.
426 sigma(1) = gd(3, 1)*gd(3, 1)
427 if (nspin > 1) then
428 sigma(2) = gd(3, 1) * gd(3, 2)
429 sigma(3) = gd(3, 2) * gd(3, 2)
430 end if
431
432 if (.not. xc_on_host) then
433 call accel_write_buffer(d_buff, nspin, d)
434 call accel_get_device_pointer(p_d, d_buff, [1])
435 call accel_get_device_pointer(p_epsx, epsx_buff, [1])
436 call accel_get_device_pointer(p_epsc, epsc_buff, [1])
437 call accel_get_device_pointer(p_dexdd, dexdd_buff, [1])
438 call accel_get_device_pointer(p_decdd, decdd_buff, [1])
439 if (gga) then
440 call accel_write_buffer(sigma_buff, 3, sigma)
441 call accel_get_device_pointer(p_sigma, sigma_buff, [1])
442 call accel_get_device_pointer(p_vxsigma, vxsigma_buff, [1])
443 call accel_get_device_pointer(p_vcsigma, vcsigma_buff, [1])
444 end if
445 else
446 p_d => d
447 p_epsx => epsx
448 p_epsc => epsc
449 p_dexdd => dexdd
450 p_decdd => decdd
451 if (gga) then
452 p_sigma => sigma
453 p_vxsigma => vxsigma
454 p_vcsigma => vcsigma
455 end if
456 end if
457
458 ! Find exchange and correlation energy densities and their
459 ! derivatives with respect to density and density gradient
460 if (gga) then
461 call xc_f03_gga_exc_vxc(x_conf, int(1, c_size_t), p_d, p_sigma, p_epsx, p_dexdd, p_vxsigma)
462 call xc_f03_gga_exc_vxc(c_conf, int(1, c_size_t), p_d, p_sigma, p_epsc, p_decdd, p_vcsigma)
463 else
464 call xc_f03_lda_exc_vxc(x_conf, int(1, c_size_t), p_d, p_epsx, p_dexdd)
465 call xc_f03_lda_exc_vxc(c_conf, int(1, c_size_t), p_d, p_epsc, p_decdd)
466 end if
467
468 if (.not. xc_on_host) then
469 call accel_finish()
470 call accel_read_buffer(dexdd_buff, nspin, dexdd)
471 call accel_read_buffer(decdd_buff, nspin, decdd)
472 call accel_read_buffer(epsx_buff, 1, epsx)
473 call accel_read_buffer(epsc_buff, 1, epsc)
474 if (gga) then
475 call accel_read_buffer(vxsigma_buff, 3, vxsigma)
476 call accel_read_buffer(vcsigma_buff, 3, vcsigma)
477 end if
478 end if
479
480 ! The derivatives of the exchange and correlation energies per particle with
481 ! respect to the density gradient have to be recovered from the derivatives
482 ! with respect to the sigma values defined above.
483 if (gga) then
484 dexdgd(:, 1) = m_two * vxsigma(1) * gd(:, 1)
485 decdgd(:, 1) = m_two * vcsigma(1) * gd(:, 1)
486 if (nspin == 2) then
487 dexdgd(:, 1) = dexdgd(:, 1) + vxsigma(2)*gd(:, 2)
488 dexdgd(:, 2) = m_two * vxsigma(3) * gd(:, 2) + vxsigma(2) * gd(:, 1)
489 decdgd(:, 1) = decdgd(:, 1) + vcsigma(2)*gd(:, 2)
490 decdgd(:, 2) = m_two * vcsigma(3) * gd(:, 2) + vcsigma(2) * gd(:, 1)
491 end if
492 end if
493
494 ! Add contributions to exchange-correlation energy and its
495 ! derivatives with respect to density at all points
496 do is = 1, nspin
497 ex = ex + dvol * d(is) * epsx(1)
498 ec = ec + dvol * d(is) * epsc(1)
499 dx = dx + dvol * d(is) * (epsx(1) - dexdd(is))
500 dc = dc + dvol * d(is) * (epsc(1) - decdd(is))
501 if (gga) then
502 v_xc(ir,is) = v_xc(ir,is) + dvol * (dexdd(is) + decdd(is))
503 do in = in1,in2
504 dx = dx - dvol * dens(ir+in,is) * dexdgd(3,is) * dgidfj(in)
505 dc = dc - dvol * dens(ir+in,is) * decdgd(3,is) * dgidfj(in)
506 v_xc(ir+in,is) = v_xc(ir+in,is) + dvol * &
507 (dexdgd(3,is) + decdgd(3,is)) * dgidfj(in)
508 end do
509 else
510 v_xc(ir,is) = dexdd(is) + decdd(is)
511 end if
512 end do
513
514 end do
515
516 if(.not. xc_on_host) then
517 call accel_free_buffer(d_buff)
518 call accel_free_buffer(dexdd_buff)
519 call accel_free_buffer(decdd_buff)
520 call accel_free_buffer(epsx_buff)
521 call accel_free_buffer(epsc_buff)
522
523 if(gga) then
524 call accel_free_buffer(sigma_buff)
525 call accel_free_buffer(vxsigma_buff)
526 call accel_free_buffer(vcsigma_buff)
527 end if
528 end if
529
530!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
531! Divide by volume element to obtain the potential (per electron) !
532!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
533
534 if (gga) then
535 do is = 1,nspin
536 do ir = 1,nr
537 dvol = aux(ir)
538 v_xc(ir,is) = v_xc(ir,is) / dvol
539 end do
540 end do
541 end if
542
543!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
544! Divide by energy unit !
545!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
546
547 ex = ex / eunit
548 ec = ec / eunit
549 dx = dx / eunit
550 dc = dc / eunit
551 do is = 1,nspin
552 do ir = 1,nr
553 v_xc(ir,is) = v_xc(ir,is) / eunit
554 end do
555 end do
556
557 call xc_f03_func_end(x_conf)
558 call xc_f03_func_end(c_conf)
559
560 pop_sub(atomxc)
561 call profiling_out("ATOMXC")
562
563 end subroutine atomxc
564
565
566!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
585 subroutine vhrtre(rho, v, r, drdi, srdrdi, nr, a)
586 real(real64), intent(in) :: rho(:)
587 real(real64), intent(out) :: v(:)
588 real(real64), intent(in) :: r(:), drdi(:), srdrdi(:)
589 integer, intent(in) :: nr
590 real(real64), intent(in) :: a
591
592 real(real64) :: x,y,q,a2by4,ybyq,qbyy,qpartc,v0,qt,dz,t,beta,dv
593 integer :: nrm1,nrm2,ir
594
595 push_sub(vhrtre)
596
597 nrm1 = nr - 1
598 nrm2 = nr - 2
599 a2by4 = a*a/m_four
600 ybyq = m_one - a*a/48.0_real64
601 qbyy = m_one/ybyq
602
603!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
604! SIMPSON`S RULE IS USED TO PERFORM TWO INTEGRALS OVER THE ELECTRON !
605! DENSITY. THE TOTAL CHARGE QT IS USED TO FIX THE POTENTIAL AT R=R(NR) !
606! AND V0 (THE INTEGRAL OF THE ELECTRON DENSITY DIVIDED BY R) FIXES !
607! THE ELECTROSTATIC POTENTIAL AT THE ORIGIN !
608! Modified to be consistent with pseudopotential generation (use the !
609! trapezoidal rule for integration). A. Rubio. Jan. 2000 !
610!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
611
612 v0 = m_zero
613 qt = m_zero
614 do ir = 2, nrm1, 2
615 dz = drdi(ir)*rho(ir)
616 qt = qt + dz
617 v0 = v0 + dz/r(ir)
618 end do
619
620 do ir=3, nrm2, 2
621 dz = drdi(ir)*rho(ir)
622 qt = qt + dz
623 v0 = v0 + dz/r(ir)
624 end do
625 dz = drdi(nr)*rho(nr)
626
627 qt = (qt + qt + dz)*m_half
628 v0 = (v0 + v0 + dz/r(nr))
629 v(1) = m_two*v0
630
631!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
632! THE ELECTROSTATIC POTENTIAL AT R=0 IS SET EQUAL TO !
633! THE AVERAGE VALUE OF RHO(R)/R !
634! BEGIN CONSTRUCTION OF THE POTENTIAL AT FINITE !
635!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
636
637 ir = 2
638 t = srdrdi(ir)/r(ir)
639 beta = drdi(ir)*t*rho(ir)
640
641!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
642! THE NEXT 4 STATEMENTS INDICATE THAT WE FIRST FIND THE PARTICULAR !
643! SOLUTION TO THE INHOMOGENEOUS EQN. FOR WHICH Q(2)=0, WE THEN !
644! ADD TO THIS PARTICULAR SOLUTION A SOLUTION OF THE HOMOGENEOUS EQN. !
645! (A CONSTANT IN V OR A Q PROPORTIONAL TO R) !
646! WHICH WHEN DIVIDED BY R IN GOING FROM Q TO V GIVES !
647! THE POTENTIAL THE DESIRED COULOMB TAIL OUTSIDE THE ELECTRON DENSITY. !
648! THE SIGNIFICANCE OF THE SOLUTION VANISHING AT THE SECOND RADIAL !
649! MESH POINT IS THAT, SINCE ALL REGULAR SOLUTIONS OF THE EQUATION !
650! FOR Q=R*V VANISH AT THE ORIGIN, THE KNOWLEDGE OF THE SOLUTION !
651! VALUE AT THE SECOND MESH POINT PROVIDES THE TWO SOLUTION VALUES !
652! REQUIRED TO START THE NUMEROV PROCEDURE. !
653!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
654
655 x = m_zero
656 y = m_zero
657 q = (y - beta / 12.0_real64) * qbyy
658 v(ir) = m_two * t * q
659
660!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
661! BEGINNING OF THE NUMEROV ALGORITHM !
662!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
663
664 do ir = 2, nr
665 x = x + a2by4 * q - beta
666 y = y + x
667 t = srdrdi(ir) / r(ir)
668 beta = t * drdi(ir) * rho(ir)
669 q = (y - beta / 12.0_real64) * qbyy
670 v(ir) = m_two * t * q
671 end do
672
673!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
674! END OF THE NUMEROV ALGORITHM !
675! !
676! WE HAVE NOW FOUND A PARTICULAR SOLUTION TO THE INHOMOGENEOUS EQN. !
677! FOR WHICH Q(R) AT THE SECOND RADIAL MESH POINT EQUALS ZERO. !
678! NOTE THAT ALL REGULAR SOLUTIONS TO THE EQUATION FOR Q=R*V !
679! VANISH AT THE ORIGIN. !
680!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
681
682 qpartc = r(nr)*v(nr)/m_two
683 dz = qt - qpartc
684 dv = m_two*dz/r(nr)
685
686!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
687! THE LOOP FOLLOWING ADDS THE CONSTANT SOLUTION OF THE HOMOGENEOUS !
688! EQN TO THE PARTICULAR SOLUTION OF THE INHOMOGENEOUS EQN. !
689! NOTE THAT V(1) IS CONSTRUCTED INDEPENDENTLY !
690!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
691
692 do ir = 2, nr
693 v(ir) = v(ir) + dv
694 end do
695
696 pop_sub(vhrtre)
697 end subroutine vhrtre
698
699!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
700! FROM NOW, ON, ALL SUBROUTINES ARE IN DOUBLE PRECISION, NO MATTER IF THE CODE
701! IS COMPILED WITH SINGLE PRECISION OR NOT. APPARENTLY DOUBLE PRECISION IS
702! NEEDED FOR THESE PROCEDURES.
703!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
704
705!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
723 subroutine egofv(namespace, h, s, n, e, g, y, l, z, a, b, rmax, nprin, nnode, dr, ierr)
724 type(namespace_t), intent(in) :: namespace
725 real(real64), intent(in) :: h(:), s(:)
726 integer, intent(in) :: n
727 real(real64), intent(inout) :: e
728 real(real64), intent(out) :: g(:), y(:)
729 integer, intent(in) :: l
730 real(real64), intent(in) :: z, a, b, rmax
731 integer, intent(in) :: nprin, nnode
732 real(real64), intent(in) :: dr
733 integer, intent(out) :: ierr
734
735 integer :: i,ncor,n1,n2,niter,nt
736 real(real64) :: e1, e2, del, de, et, t
737 real(real64), parameter :: tol = 1.0e-5_real64
738
739 push_sub(egofv)
740
741 ncor = nprin-l-1
742 n1 = nnode
743 n2 = nnode-1
744 e1 = e
745 e2 = e
746 ierr = 1
747
748!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
749! !
750! the labels 1 and 2 refer to the bisection process, defining the !
751! range in which the desired solution is located. the initial !
752! settings of n1, n2, e1 and e2 are not consistent with the bisection !
753! algorithm; they are set to consistent values when the desired !
754! energy interval has been located. !
755! !
756!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
757
758 nt = 0
759 del = m_half
760 de = m_zero
761
762 do niter = 0, 40
763
764 et = e + de
765 ! the following line is the fundamental "bisection"
766 e = m_half*(e1+e2)
767
768!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
769! !
770! the following concatenation of logical ors ensures that node !
771! counting is used unless the previous integration of the radial !
772! eq produced both the correct number of nodes and a sensible !
773! prediction for the energy. !
774! !
775! sensible means that et must be greater than e1 and less than e2 !
776! correct number of nodes means that nt = nnode or nnode-1. !
777! !
778! leaving e set to its bisection value, and transfering to !
779! the call to yofe means that we are performing bisection, !
780! whereas setting e to et is use of the variational estimate. !
781! !
782!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
783
784
785 if (.not. (et <= e1 .or. et >= e2 .or. nt < nnode-1 .or. nt > nnode)) then
786 e = et
787 if (abs(de) < tol) then
788 ierr = 0
789 exit
790 end if
791 end if
792 call yofe(e,de,dr,rmax,h,s,y,n,l,ncor,nt,z,a,b)
793
794
795!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
796! !
797! yofe integrates the schro eq.; now the bisection logic !
798! !
799!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
800
801 if (nt < nnode) then
802 ! too few nodes; set e1 and n1
803 e1 = e
804 n1 = nt
805
806!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
807! !
808! at this point, we have just set the bottom of the bisection range; !
809! if the top is also set, we proceed. if the top of the range has not !
810! been set, it means that we have yet to find an e greater than the !
811! desired energy. the upper end of the range is extended. !
812! !
813!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
814
815 if (n2 >= nnode) cycle
816 del=del*m_two
817 e2=e1+del
818 cycle
819 end if
820
821 ! too many nodes; set e2 and n2
822 e2=e
823 n2=nt
824
825
826!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
827! !
828! at this point, we have just set the top of the bisection range; !
829! if the top is also set, we proceed. if the top of the range has !
830! not been set, it means that we have yet to find an e less than the !
831! desired energy. the lower end of the range is extended. !
832! !
833!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
834
835 if (n1 < nnode) cycle
836 del=del*m_two
837 e1=e2-del
838
839 end do
840
841!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
842! !
843! the Numerov method uses a transformation of the radial wavefunction. !
844! that we call "y". having located the eigenenergy, we transform !
845! y to "g", from which the density is easily constructed. !
846! finally, the call to "nrmlzg" normalizes g to one electron. !
847! !
848!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
849
850 if (ierr /= 1) then
851 g(1) = m_zero
852 do i = 2, n
853 t = h(i) - e * s(i)
854 g(i) = y(i) / (m_one - t / 12._real64)
855 end do
856 call nrmlzg(namespace, g, s, n)
857 end if
858
859 pop_sub(egofv)
860 end subroutine egofv
861
862
863 subroutine yofe(e,de,dr,rmax,h,s,y,nmax,l,ncor,nnode,z,a,b)
864
865
866!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
883
884 real(real64), intent(inout) :: e
885 real(real64), intent(out) :: de
886 real(real64), intent(in) :: dr, rmax
887 real(real64), intent(in) :: h(:), s(:)
888 real(real64), intent(out) :: y(:)
889 integer, intent(in) :: nmax, l
890 integer, intent(in) :: ncor
891 integer, intent(out) :: nnode
892 real(real64), intent(in) :: z, a, b
893
894 real(real64) :: y2, g, gsg, x, gin, gsgin, xin
895 integer :: n,knk,nndin,i
896 real(real64) :: zdr, yn, ratio, t
897
898 ! No PUSH SUB, called too often.
899
900 zdr = z*a*b
901
902 do n = nmax, 1, -1
903 if (h(n)-e*s(n) < m_one) then
904 knk = n
905 exit
906 end if
907 y(n) = m_zero
908 end do
909
910 call bcorgn(e,h,s,l,zdr,y2)
911
912
913!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
914! !
915! bcorgn computes y2, which embodies the boundary condition !
916! satisfied by the radial wavefunction at the origin !
917!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
918
919 call numout(e,h,s,y,ncor,knk,nnode,y2,g,gsg,x)
920
921!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
922! !
923! the outward integration is now complete !
924! !
925! we first decide if the kinetic energy is sufficiently non !
926! negative to permit use of the numerov eq at rmax. if !
927! it is not, then zero-value boundary condition is used !
928! !
929!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
930
931
932 yn = m_zero
933 if (.not. (n < nmax .or. abs(dr) > 1.e3_real64)) then
934 call bcrmax(e,dr,rmax,h,s,n,yn,a,b)
935 end if
936
937
938!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
939! !
940! bcrmax computes yn, which embodies the boundary condition !
941! satisfied by the radial wavefunction at rmax !
942! !
943!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
944
945 ! n should be larger than 1
946 n = max(n, 2)
947 knk = max(knk, 1)
948 call numin(e,h,s,y,n,nndin,yn,gin,gsgin,xin,knk)
949
950
951!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
952! !
953! numin performs the inward integration by the Numerov method !
954! !
955! the energy increment is now evaluated from the kink in psi !
956! !
957!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
958
959
960 ratio = g / gin
961 xin = xin * ratio
962 gsg = gsg + gsgin * ratio * ratio
963 t = h(knk) - e * s(knk)
964 de = g * (x + xin + t * g) / gsg
965 nnode = nnode + nndin
966 if (de < m_zero) nnode = nnode + 1
967
968 do i = knk, n
969 y(i) = y(i) * ratio
970 end do
971
972 end subroutine yofe
973
974
975 subroutine nrmlzg(namespace, g, s, n)
976
977
978!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
997 type(namespace_t), intent(in) :: namespace
998 real(real64), intent(inout) :: g(:)
999 real(real64), intent(in) :: s(:)
1000 integer, intent(in) :: n
1001
1002 integer :: nm1,nm2,i
1003 real(real64) :: norm, srnrm
1004
1005 push_sub(nrmlzg)
1006
1007!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1008! determine the norm of g(i) using Simpson`s rule !
1009!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1010
1011
1012 if (mod(n,2) /= 1) then
1013 write(message(1),'(a,i6)') ' nrmlzg: n should be odd. n =', n
1014 call messages_warning(1, namespace=namespace)
1015 end if
1016
1017 norm = m_zero
1018 nm1 = n - 1
1019 do i = 2,nm1,2
1020 norm=norm + g(i)*s(i)*g(i)
1021 end do
1022 norm = norm * m_two
1023 nm2 = n - 2
1024 do i = 3,nm2,2
1025 norm=norm + g(i)*s(i)*g(i)
1026 end do
1027 norm = norm * m_two
1028 nm2 = n - 2
1029 do i = 1, n, nm1
1030 norm = norm + g(i) * s(i) * g(i)
1031 end do
1032 norm = norm / m_three
1033 srnrm = sqrt(norm)
1034 do i = 1, n
1035 g(i) = g(i) / srnrm
1036 end do
1037
1038 pop_sub(nrmlzg)
1039 end subroutine nrmlzg
1040
1041
1042 subroutine bcorgn(e,h,s,l,zdr,y2)
1043 real(real64), intent(in) :: e
1044 real(real64), intent(in) :: h(:), s(:)
1045 integer, intent(in) :: l
1046 real(real64), intent(in) :: zdr
1047 real(real64), intent(out) :: y2
1048
1049 real(real64) :: t2, t3, d2, c0, c1, c2
1050
1051 ! no PUSH SUB, called too often
1052
1053!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1054! the quantity called d(i) in the program is actually the inverse !
1055! of the diagonal of the tri-diagonal Numerov matrix !
1056!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1057
1058 t2 = h(2) - e * s(2)
1059 d2 = -((24._real64 + 10.0_real64 * t2) / (12._real64 - t2))
1060
1061!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1062! The following section deals with the fact that the independent !
1063! variable "y" in the Numerov equation is not zero at the origin !
1064! for l less than 2. !
1065! The l=0 solution g vanishes, but the first and second !
1066! derivatives are finite, making the Numerov variable y finite. !
1067! The l=1 solution g vanishes, and gprime also vanishes, but !
1068! The second derivative gpp is finite making y finite. For l > 1, !
1069! g and its first two derivatives vanish, making y zero. !
1070!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1071
1072 if (l < 2) then
1073 if (l <= 0) then
1074 c0=zdr/6.0_real64
1075 c0=c0/(m_one-0.75_real64*zdr)
1076 else
1077 c0=m_one/12._real64
1078 c0=(-c0)*8.0_real64/m_three
1079 end if
1080
1081 c1=c0*(12._real64+13._real64*t2)/(12._real64-t2)
1082 t3=h(3)-e*s(3)
1083 c2=(-m_half)*c0*(24._real64-t3)/(12._real64-t3)
1084 d2=(d2-c1)/(m_one-c2)
1085 end if
1086 y2=(-m_one)/d2
1087
1088
1089 if (l < 2) then
1090 if (l <= 0) then
1091 c0=zdr/6.0_real64
1092 c0=c0/(m_one-0.75_real64*zdr)
1093 else
1094 c0=m_one/12._real64
1095 c0=(-c0)*8.0_real64/m_three
1096 end if
1097 c1=c0*(12._real64+13._real64*t2)/(12._real64-t2)
1098 t3=h(3)-e*s(3)
1099 c2=(-m_half)*c0*(24._real64-t3)/(12._real64-t3)
1100 d2=(d2-c1)/(m_one-c2)
1101 end if
1102 y2=(-m_one)/d2
1103
1104 end subroutine bcorgn
1105
1106
1107!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1108! 22.7.85 !
1109!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1110 subroutine bcrmax(e, dr, rmax, h, s, n, yn, a, b)
1111 real(real64), intent(in) :: e, dr, rmax
1112 real(real64), intent(in) :: h(:), s(:)
1113 integer, intent(in) :: n
1114 real(real64), intent(out) :: yn
1115 real(real64), intent(in) :: a, b
1116
1117 real(real64) :: tnm1, tn, tnp1, beta, dg, c1, c2, c3, dn
1118
1119 push_sub(bcrmax)
1120
1121 tnm1=h(n-1)-e*s(n-1)
1122 tn =h(n )-e*s(n )
1123 tnp1=h(n+1)-e*s(n+1)
1124 beta=m_one+b/rmax
1125 dg=a*beta*(dr+m_one-m_half/beta)
1126
1127 c2=24._real64*dg/(12._real64-tn)
1128 dn=-((24._real64+10._real64*tn)/(12._real64-tn))
1129
1130 c1= (m_one-tnm1/6.0_real64)/(m_one-tnm1/12._real64)
1131 c3=-((m_one-tnp1/6.0_real64)/(m_one-tnp1/12._real64))
1132 yn=-((m_one-c1/c3)/(dn-c2/c3))
1133
1134 pop_sub(bcrmax)
1135 end subroutine bcrmax
1136
1137
1138 subroutine numin(e,h,s,y,n,nnode,yn,g,gsg,x,knk)
1139 real(real64), intent(in) :: e, h(:), s(:)
1140 real(real64), intent(inout) :: y(:)
1141 integer, intent(in) :: n
1142 integer, intent(out) :: nnode
1143 real(real64), intent(in) :: yn
1144 real(real64), intent(out) :: g, gsg, x
1145 integer, intent(in) :: knk
1146
1147 integer :: i
1148 real(real64) :: t
1149
1150 ! no PUSH SUB, called too often
1151
1152 y(n)=yn
1153 t=h(n)-e*s(n)
1154 g=y(n)/(m_one-t/12._real64)
1155 gsg=g*s(n)*g
1156 i=n-1
1157 y(i)=m_one
1158 t=h(i)-e*s(i)
1159 g=y(i)/(m_one-t/12._real64)
1160 gsg=gsg+g*s(i)*g
1161 x=y(i)-y(n)
1162 nnode=0
1163
1164
1165!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1166! !
1167! begin the inward integration by the Numerov method !
1168! !
1169!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1170
1171 assert(knk > 0)
1172
1173 do i = n - 2, knk, -1
1174 x = x + t * g
1175 y(i) = y(i+1) + x
1176 ! y(i) * y(i+1) can lead to overflows. This is now prevented
1177 ! by comparing only the signs
1178 if (sign(m_one,y(i)) * sign(m_one,y(i+1)) < m_zero) nnode = nnode + 1
1179 t = h(i) - e * s(i)
1180 g = y(i) / (m_one - t / 12._real64)
1181 gsg = gsg + g * s(i) * g
1182 end do
1183
1184
1185!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1186! !
1187! The kink radius is defined as the point where !
1188! psi first turns downward. this usually means at the outermost !
1189! maximum !
1190! !
1191! the inward integration is now complete !
1192! !
1193!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1194
1195 end subroutine numin
1196
1197
1198 subroutine numout(e, h, s, y, ncor, knk, nnode, y2, g, gsg, x)
1199 real(real64), intent(in) :: e, h(:), s(:)
1200 real(real64), intent(inout) :: y(:)
1201 integer, intent(in) :: ncor
1202 integer, intent(inout) :: knk
1203 integer, intent(out) :: nnode
1204 real(real64), intent(in) :: y2
1205 real(real64), intent(out) :: g, gsg, x
1206
1207 integer :: i, nm4
1208 real(real64) :: t, xl
1209
1210 ! no PUSH SUB, called too often
1211
1212 y(1)=m_zero
1213 y(2)=y2
1214 t=h(2)-e*s(2)
1215 g=y(2)/(m_one-t/12._real64)
1216 gsg=g*s(2)*g
1217 y(3)=m_one
1218 t=h(3)-e*s(3)
1219 g=y(3)/(m_one-t/12._real64)
1220 gsg=gsg+g*s(3)*g
1221 x=y(3)-y(2)
1222 nnode=0
1223
1224!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1225! !
1226! begin the outward integration by the Numerov method !
1227! !
1228!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1229
1230 nm4 = knk - 4
1231 knk = nm4
1232 do i = 4, nm4
1233 xl = x
1234 x = x + t * g
1235 y(i) = y(i-1) + x
1236 ! y(i) * y(i-1) can lead to overflows. This is now prevented
1237 ! by comparing only the signs
1238 if (sign(m_one,y(i)) * sign(m_one,y(i-1)) < m_zero) nnode = nnode + 1
1239 t = h(i) - e * s(i)
1240 g = y(i) / (m_one - t / 12._real64)
1241 gsg = gsg + g * s(i) * g
1242 if (.not. (nnode < ncor .or. xl*x > m_zero)) then
1243 knk = i
1244 exit
1245 end if
1246 end do
1247
1248!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1249! !
1250! the outward integration is now complete !
1251! !
1252!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1253 end subroutine numout
1254
1255end module atomic_oct_m
1256
1257!! Local Variables:
1258!! mode: f90
1259!! coding: utf-8
1260!! End:
accel_oct_m::accel_create_buffer
Definition: accel.F90:301
accel_oct_m::accel_get_device_pointer
Definition: accel.F90:364
accel_oct_m::accel_read_buffer
Definition: accel.F90:338
accel_oct_m::accel_write_buffer
Definition: accel.F90:313
sort_oct_m::sort
This is the common interface to a sorting routine. It performs the shell algorithm,...
Definition: sort.F90:156
sqrt
double sqrt(double __x) __attribute__((__nothrow__
accel_oct_m
Definition: accel.F90:120
accel_oct_m::accel_free_buffer
subroutine, public accel_free_buffer(this, async)
Definition: accel.F90:1006
accel_oct_m::accel_finish
subroutine, public accel_finish()
Definition: accel.F90:1124
accel_oct_m::accel_mem_read_write
integer, parameter, public accel_mem_read_write
Definition: accel.F90:186
atomic_oct_m
Definition: atomic.F90:116
atomic_oct_m::bcrmax
subroutine bcrmax(e, dr, rmax, h, s, n, yn, a, b)
Definition: atomic.F90:1206
atomic_oct_m::atomhxc
subroutine, public atomhxc(namespace, functl, g, nspin, dens, v, extra)
Definition: atomic.F90:233
atomic_oct_m::nrmlzg
subroutine nrmlzg(namespace, g, s, n)
Definition: atomic.F90:1071
atomic_oct_m::yofe
subroutine yofe(e, de, dr, rmax, h, s, y, nmax, l, ncor, nnode, z, a, b)
Definition: atomic.F90:959
atomic_oct_m::atomxc
subroutine, public atomxc(namespace, FUNCTL, AUTHOR, NR, MAXR, RMESH, NSPIN, DENS, EX, EC, DX, DC, V_XC)
Finds total exchange-correlation energy and potential for a ! spherical electron density distribution...
Definition: atomic.F90:342
atomic_oct_m::valconf_unpolarized_to_polarized
subroutine, public valconf_unpolarized_to_polarized(conf)
Definition: atomic.F90:213
atomic_oct_m::valconf_string_length
integer, parameter, public valconf_string_length
Following stuff is for the "val_conf_t" data type, made to handle atomic configurations.
Definition: atomic.F90:150
atomic_oct_m::valconf_copy
subroutine, public valconf_copy(cout, cin)
Definition: atomic.F90:169
atomic_oct_m::valconf_sort_by_l
subroutine, public valconf_sort_by_l(c, lmax)
Definition: atomic.F90:187
atomic_oct_m::numin
subroutine numin(e, h, s, y, n, nnode, yn, g, gsg, x, knk)
Definition: atomic.F90:1234
atomic_oct_m::numout
subroutine numout(e, h, s, y, ncor, knk, nnode, y2, g, gsg, x)
Definition: atomic.F90:1294
atomic_oct_m::bcorgn
subroutine bcorgn(e, h, s, l, zdr, y2)
Definition: atomic.F90:1138
atomic_oct_m::vhrtre
subroutine, public vhrtre(rho, v, r, drdi, srdrdi, nr, a)
VHRTRE CONSTRUCTS THE ELECTROSTATIC POTENTIAL DUE TO A SUPPLIED ! ELECTRON DENSITY....
Definition: atomic.F90:681
atomic_oct_m::egofv
subroutine, public egofv(namespace, h, s, n, e, g, y, l, z, a, b, rmax, nprin, nnode, dr, ierr)
egofv determines the eigenenergy and wavefunction corresponding ! to a particular l,...
Definition: atomic.F90:819
debug_oct_m
Definition: debug.F90:116
global_oct_m
Definition: global.F90:116
global_oct_m::m_two
real(real64), parameter, public m_two
Definition: global.F90:202
global_oct_m::m_zero
real(real64), parameter, public m_zero
Definition: global.F90:200
global_oct_m::m_four
real(real64), parameter, public m_four
Definition: global.F90:204
global_oct_m::m_pi
real(real64), parameter, public m_pi
some mathematical constants
Definition: global.F90:198
global_oct_m::m_epsilon
real(real64), parameter, public m_epsilon
Definition: global.F90:216
global_oct_m::m_half
real(real64), parameter, public m_half
Definition: global.F90:206
global_oct_m::m_one
real(real64), parameter, public m_one
Definition: global.F90:201
global_oct_m::m_three
real(real64), parameter, public m_three
Definition: global.F90:203
logrid_oct_m
Definition: logrid.F90:116
messages_oct_m
Definition: messages.F90:117
messages_oct_m::messages_warning
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:525
messages_oct_m::message
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:162
messages_oct_m::messages_fatal
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:410
namespace_oct_m
Definition: namespace.F90:105
profiling_oct_m
Definition: profiling.F90:118
profiling_oct_m::profiling_out
subroutine, public profiling_out(label)
Increment out counter and sum up difference between entry and exit time.
Definition: profiling.F90:631
profiling_oct_m::profiling_in
subroutine, public profiling_in(label, exclude)
Increment in counter and save entry time.
Definition: profiling.F90:554
sort_oct_m
This module is intended to contain "only mathematical" functions and procedures.
Definition: sort.F90:119
types_oct_m
Definition: types.F90:116
types_oct_m::type_float
type(type_t), parameter, public type_float
Definition: types.F90:135
xc_functional_oct_m
Definition: xc_functional.F90:116
xc_functional_oct_m::xc_functional_init_func
subroutine, public xc_functional_init_func(conf, func_id, nspin, on_cpu)
Definition: xc_functional.F90:1488
xc_oct_m
Definition: xc.F90:120
xc_oct_m::xc_is_using_device
logical function, public xc_is_using_device(namespace)
Check if XC will run on device (parses config and checks availability)
Definition: xc.F90:328
accel_oct_m::accel_mem_t
Definition: accel.F90:229
atomic_oct_m::valconf_t
Definition: atomic.F90:152
namespace_oct_m::namespace_t
Definition: namespace.F90:117
true
int true(void)
Definition: symmetries_finite.c:3150