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 ps_in_grid_oct_m
22  use atomic_oct_m
23  use global_oct_m
24  use logrid_oct_m
25  use messages_oct_m
26  use namespace_oct_m
27  use profiling_oct_m
28
29  implicit none
30
31  private
32  public ::                   &
33    ps_in_grid_t,             &
34    ps_in_grid_init,          &
35    ps_in_grid_end,           &
36    ps_in_grid_vlocal,        &
37    ps_in_grid_kb_projectors, &
38    ps_in_grid_kb_cosines,    &
39    ps_in_grid_cutoff_radii,  &
40    ps_in_grid_check_rphi,    &
41    first_point_extrapolate
42
43  type ps_in_grid_t
44    ! Components are public by default
45    type(logrid_t) :: g                !< log grid where the pseudos are defined
46
47    FLOAT          :: zval             !< valence charge
48
49    integer        :: no_l_channels    !< number of l channels to consider
50    FLOAT, pointer :: vps(:, :)        !< the pseudopotential (l=0 .. no_l_channels-1)
51    FLOAT, pointer :: KB(:,:)          !< Kleinman-Bylander projectors
52    FLOAT, pointer :: dkbcos(:)        !< Kleinman-Bylander cosine
53    FLOAT, pointer, private :: dknorm(:)        !< Kleinman-Bylander norm
54    FLOAT, pointer :: kb_radius(:)     !< radius of KB projectors
55
56    integer        :: so_no_l_channels
57    FLOAT, pointer :: so_vps(:,:)      !< spin-orbit components (l=1 .. so_no_l_channels)
58    FLOAT, pointer, private :: so_KB(:,:)       !< Kleinman-Bylander projectors
59    FLOAT, pointer, private :: so_dkbcos(:)     !< Kleinman-Bylander cosine
60    FLOAT, pointer, private :: so_dknorm(:)     !< Kleinman-Bylander norm
61    FLOAT, pointer, private :: so_kb_radius(:)  !< radius of KB projectors
62
63    FLOAT, pointer :: vlocal(:)        !< local part of the pseudopotential
64    FLOAT, pointer :: rphi(:, :,:)     !< pseudo wavefunctions
65
66    logical        :: core_corrections
67    FLOAT, pointer :: chcore(:)        !< core charge density
68
69  end type ps_in_grid_t
70
71contains
72
73  subroutine ps_in_grid_init(ps, flavor, a, b, nrval, no_l, so_no_l)
74    type(ps_in_grid_t), intent(out) :: ps
75    integer,            intent(in)  :: flavor, nrval
76    FLOAT,              intent(in)  :: a, b
77    integer,            intent(in)  :: no_l, so_no_l
78
79    PUSH_SUB(ps_in_grid_init)
80
81    ! initialize logarithmic grid
82    call logrid_init(ps%g, flavor, a, b, nrval)
83
84    ! copy data
85    ps%   no_l_channels =    no_l
86    ps%so_no_l_channels = so_no_l
87
88    ! Allocate some stuff
89    SAFE_ALLOCATE(ps%vps    (1:ps%g%nrval, 1:no_l))
90    SAFE_ALLOCATE(ps%chcore (1:ps%g%nrval))
91    SAFE_ALLOCATE(ps%vlocal (1:ps%g%nrval))
92
93    SAFE_ALLOCATE(ps%rphi  (1:ps%g%nrval, 1:no_l, 1:3))
94    SAFE_ALLOCATE(ps%KB    (1:ps%g%nrval, 1:no_l))
95    SAFE_ALLOCATE(ps%dkbcos(1:no_l))
96    SAFE_ALLOCATE(ps%dknorm(1:no_l))
97    SAFE_ALLOCATE(ps%kb_radius(1:no_l+1))
98
99    if(so_no_l > 0) then
100      SAFE_ALLOCATE(ps%so_vps(1:ps%g%nrval, 1:so_no_l))
101      SAFE_ALLOCATE(ps%so_KB (1:ps%g%nrval, 1:so_no_l))
102      SAFE_ALLOCATE(ps%so_dkbcos(1:so_no_l))
103      SAFE_ALLOCATE(ps%so_dknorm(1:so_no_l))
104      SAFE_ALLOCATE(ps%so_kb_radius(1:so_no_l))
105    end if
106
107    POP_SUB(ps_in_grid_init)
108  end subroutine ps_in_grid_init
109
110
111  ! ---------------------------------------------------------
112  subroutine ps_in_grid_end(ps)
113    type(ps_in_grid_t), intent(inout) :: ps
114
115    PUSH_SUB(ps_in_grid_end)
116
117    SAFE_DEALLOCATE_P(ps%vps)
118    SAFE_DEALLOCATE_P(ps%chcore)
119    SAFE_DEALLOCATE_P(ps%vlocal)
120
121    SAFE_DEALLOCATE_P(ps%rphi)
122    SAFE_DEALLOCATE_P(ps%KB)
123    SAFE_DEALLOCATE_P(ps%dkbcos)
124    SAFE_DEALLOCATE_P(ps%dknorm)
125    SAFE_DEALLOCATE_P(ps%kb_radius)
126
127    if(ps%so_no_l_channels > 0) then
128      SAFE_DEALLOCATE_P(ps%so_vps)
129      SAFE_DEALLOCATE_P(ps%so_KB)
130      SAFE_DEALLOCATE_P(ps%so_dkbcos)
131      SAFE_DEALLOCATE_P(ps%so_dknorm)
132      SAFE_DEALLOCATE_P(ps%so_kb_radius)
133    end if
134
135    call logrid_end(ps%g)
136
137    POP_SUB(ps_in_grid_end)
138  end subroutine ps_in_grid_end
139
140
141  ! ---------------------------------------------------------
142  subroutine ps_in_grid_vlocal(ps, l_loc, rcore, namespace)
143    type(ps_in_grid_t), intent(inout) :: ps
144    integer,            intent(in)    :: l_loc
145    FLOAT,              intent(in)    :: rcore
146    type(namespace_t),  intent(in)    :: namespace
147
148    integer :: ir
149    FLOAT :: a, b, qtot
150    FLOAT, allocatable :: rho(:)
151
152    PUSH_SUB(ps_in_grid_vlocal)
153
154    if(l_loc >= 0) then
155      ps%vlocal(:) = ps%vps(:, l_loc+1)
156
157    else if(l_loc == -1) then
158      if(ps%g%flavor /= LOGRID_PSF) then
159        message(1) = "For the moment, Vanderbilt local potentials are only possible with tm grids."
160        call messages_fatal(1, namespace=namespace)
161      end if
162
163      a = CNST(1.82) / rcore
164      b = M_ONE
165      SAFE_ALLOCATE(rho(1:ps%g%nrval))
166
167      do ir = 1, ps%g%nrval
168        rho(ir) = exp( -( sinh(a*b*ps%g%rofi(ir)) / sinh(b) )**2 )
169        rho(ir) = M_FOUR * M_PI * rho(ir) * ps%g%rofi(ir)**2
170      end do
171!      do ir =1,2
172      qtot = sum(rho(:)*ps%g%drdi(:))
173      rho(:) = - rho(:)*(ps%zval/qtot)
174
175      call vhrtre(rho, ps%vlocal, ps%g%rofi, ps%g%drdi, ps%g%s, ps%g%nrval, ps%g%a)
176      ps%vlocal(1) = ps%vlocal(2)
177
178      SAFE_DEALLOCATE_A(rho)
179    end if
180
181    POP_SUB(ps_in_grid_vlocal)
182  end subroutine ps_in_grid_vlocal
183
184
185  ! ---------------------------------------------------------
186  !> KB-projectors
187  !!   kb = (vps - vlocal) |phi> * dknorm
188  subroutine ps_in_grid_kb_projectors(ps)
189    type(ps_in_grid_t), intent(inout) :: ps
190
191    integer :: l
192
193    PUSH_SUB(ps_in_grid_kb_projectors)
194
195    do l = 1, ps%no_l_channels
196      ps%KB(2:, l) = (ps%vps(2:, l) - ps%vlocal(2:))*(ps%rphi(2:, l, 1)/ps%g%rofi(2:))*ps%dknorm(l)
197
198      ps%KB(1, l) = first_point_extrapolate(ps%g%rofi, ps%KB(:, l))
199    end do
200
201    do l = 1, ps%so_no_l_channels
202      ps%so_KB(2:, l) = ps%so_vps(2:, l)*(ps%rphi(2:, l, 1)/ps%g%rofi(2:))*ps%dknorm(l)
203
204      ps%so_KB(1, l) = first_point_extrapolate(ps%g%rofi, ps%so_KB(:, l))
205    end do
206
207    POP_SUB(ps_in_grid_kb_projectors)
208  end subroutine ps_in_grid_kb_projectors
209
210
211  ! ---------------------------------------------------------
212  !> KB-cosines and KB-norms:
213  !!       dkbcos stores the KB "cosines:"
214  !!               || (v_l - v_local) phi_l ||^2 / < (v_l - v_local)phi_l | phi_l >  [Rydberg]
215  !!       dknorm stores the KB "norms:"
216  !!               1 / || (v_l - v_local) phi_l || [1/Rydberg]
217  subroutine ps_in_grid_kb_cosines(ps, lloc)
218    type(ps_in_grid_t), intent(inout) :: ps
219    integer,            intent(in)    :: lloc
220
221    integer :: ir, l
222    FLOAT :: dnrm, avgv, vphi
223
224    PUSH_SUB(ps_in_grid_kb_cosines)
225
226    do l = 1, ps%no_l_channels
227      if(l-1 == lloc) then
228        ps%dkbcos(l) = M_ZERO
229        ps%dknorm(l) = M_ZERO
230        cycle
231      end if
232
233      dnrm = M_ZERO
234      avgv = M_ZERO
235      do ir = 1, ps%g%nrval
236        vphi = (ps%vps(ir, l) - ps%vlocal(ir))*ps%rphi(ir, l, 1)
237        dnrm = dnrm + vphi*vphi*ps%g%drdi(ir)
238        avgv = avgv + vphi*ps%rphi(ir, l, 1)*ps%g%drdi(ir)
239      end do
240      ps%dkbcos(l) = dnrm/(avgv + CNST(1.0e-20))
241      ps%dknorm(l) = M_ONE/(sqrt(dnrm) + CNST(1.0e-20))
242    end do
243
244    do l = 1, ps%so_no_l_channels
245      dnrm = M_ZERO
246      avgv = M_ZERO
247      do ir = 1, ps%g%nrval
248        vphi = ps%so_vps(ir, l)*ps%rphi(ir, l, 1)
249        dnrm = dnrm + vphi*vphi*ps%g%drdi(ir)
250        avgv = avgv + vphi*ps%rphi(ir, l, 1)*ps%g%drdi(ir)
251      end do
252      ps%so_dkbcos(l) = dnrm/(avgv + CNST(1.0e-20))
253      ps%so_dknorm(l) = M_ONE/(sqrt(dnrm) + CNST(1.0e-20))
254    end do
255
256    POP_SUB(ps_in_grid_kb_cosines)
257  end subroutine ps_in_grid_kb_cosines
258
259
260  ! ---------------------------------------------------------
261  subroutine ps_in_grid_cutoff_radii(ps, lloc)
262    type(ps_in_grid_t), intent(inout) :: ps
263    integer,            intent(in)    :: lloc
264
265    integer          :: l, ir
266    FLOAT            :: dincv, phi
267    FLOAT, parameter :: threshold = CNST(1.0e-6)
268
269    PUSH_SUB(ps_in_grid_cutoff_radii)
270
271    ! local part ....
272    do ir = ps%g%nrval-1, 2, -1
273      dincv = abs(ps%vlocal(ir)*ps%g%rofi(ir) + M_TWO*ps%zval)
274      if(dincv > threshold) exit
275    end do
276    ps%kb_radius(ps%no_l_channels+1) = ps%g%rofi(ir + 1)
277
278    ! non-local part....
279    do l = 1, ps%no_l_channels
280      if(l-1 == lloc) then
281        ps%kb_radius(l) = M_ZERO
282        cycle
283      end if
284
285      do ir = ps%g%nrval-1, 2, -1
286        phi = (ps%rphi(ir, l, 1)/ps%g%rofi(ir))*ps%dknorm(l)
287        dincv = abs((ps%vps(ir, l) - ps%vlocal(ir))*phi)
288
289        if(dincv > threshold) exit
290
291        phi = (ps%rphi(ir, l, 1)/ps%g%rofi(ir))*ps%dknorm(l)
292      end do
293
294      ps%kb_radius(l) = ps%g%rofi(ir + 1)
295    end do
296
297    ! now the SO part
298    do l = 1, ps%so_no_l_channels
299      do ir = ps%g%nrval-1, 2, -1
300        phi = (ps%rphi(ir, l, 1)/ps%g%rofi(ir))*ps%so_dknorm(l)
301        dincv = abs((ps%so_vps(ir, l))*phi)
302
303        if(dincv > threshold) exit
304      end do
305
306      ps%so_kb_radius(l) = ps%g%rofi(ir + 1)
307    end do
308
309    POP_SUB(ps_in_grid_cutoff_radii)
310  end subroutine ps_in_grid_cutoff_radii
311
312
313  ! ---------------------------------------------------------
314  !> checks normalization of the pseudo wavefunctions
315  subroutine ps_in_grid_check_rphi(ps, namespace)
316    type(ps_in_grid_t), intent(in) :: ps
317    type(namespace_t),  intent(in) :: namespace
318
319    integer :: l
320    FLOAT   :: nrm
321
322    PUSH_SUB(ps_in_grid_check_rphi)
323
324    !  checking normalization of the wavefunctions
325    do l = 1, ps%no_l_channels
326      nrm = sqrt(sum(ps%g%drdi(:)*ps%rphi(:, l, 1)**2))
327      nrm = abs(nrm - M_ONE)
328      if (nrm > CNST(1.0e-5)) then
329        write(message(1), '(a,i2,a)') "Eigenstate for l = ", l-1, ' is not normalized.'
330        write(message(2), '(a, f12.6,a)') '(abs(1 - norm) = ', nrm, ')'
331        call messages_warning(2, namespace=namespace)
332      end if
333    end do
334
335    POP_SUB(ps_in_grid_check_rphi)
336  end subroutine ps_in_grid_check_rphi
337
338  ! ---------------------------------------------------------
339
340  FLOAT function first_point_extrapolate(x, y, high_order) result(y0)
341    FLOAT,             intent(in)    :: x(:)
342    FLOAT,             intent(in)    :: y(:)
343    logical, optional, intent(in)    :: high_order
344
345    FLOAT :: x1, x2, x3
346    FLOAT :: y1, y2, y3
347
348    x1 = x(2) - x(1)
349    x2 = x(3) - x(1)
350    x3 = x(4) - x(1)
351    y1 = y(2)
352    y2 = y(3)
353    y3 = y(4)
354
355    if(optional_default(high_order, .false.)) then
356
357      y0 = y1*x2*x3*(x2 - x3) + y2*x1*x3*(x3 - x1) + y3*x1*x2*(x1 - x2);
358      y0 = y0/((x1 - x2)*(x1 - x3)*(x2 - x3));
359
360    else
361
362      y0 = y1 - (y2 - y1)*x1/(x2 - x1)
363
364    end if
365
366  end function first_point_extrapolate
367
368end module ps_in_grid_oct_m
369
370!! Local Variables:
371!! mode: f90
372!! coding: utf-8
373!! End:
374