1!
2! Copyright (C) 2004-2007 Quantum ESPRESSO group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!=----------------------------------------------------------------------------=!
9  MODULE upf_to_internal
10!=----------------------------------------------------------------------------=!
11  !! Contains two subroutines:
12  !! add_upf_grid  generates the radial grid from data contained in upf
13  !! set_upf_q     builds the Q(r) functions if not explicitly present
14  !!
15  USE pseudo_types
16  IMPLICIT NONE
17  PRIVATE
18  PUBLIC :: add_upf_grid, set_upf_q
19  SAVE
20
21!=----------------------------------------------------------------------------=!
22  CONTAINS
23!=----------------------------------------------------------------------------=!
24!
25!---------------------------------------------------------------------
26SUBROUTINE add_upf_grid (upf, grid)
27  !---------------------------------------------------------------------
28  !
29  USE radial_grids, ONLY: radial_grid_type, allocate_radial_grid
30  !
31  IMPLICIT NONE
32  !
33  TYPE (pseudo_upf), INTENT(in) :: upf
34  TYPE (radial_grid_type), INTENT(out) :: grid
35  !
36  CALL allocate_radial_grid(grid,upf%mesh)
37  !
38  grid%dx   = upf%dx
39  grid%xmin = upf%xmin
40  grid%zmesh= upf%zmesh
41  grid%mesh = upf%mesh
42  grid%r  (1:upf%mesh) = upf%r  (1:upf%mesh)
43  grid%rab(1:upf%mesh) = upf%rab(1:upf%mesh)
44  !
45  grid%r2 = upf%r**2
46  grid%sqr= sqrt(upf%r)
47  ! Prevent FP error if r(1) = 0
48  IF ( upf%r(1) > 1.0D-16 ) THEN
49     grid%rm1 = upf%r**(-1)
50     grid%rm2 = upf%r**(-2)
51     grid%rm3 = upf%r**(-3)
52  ELSE
53     grid%rm1(1) =0.0_dp
54     grid%rm2(1) =0.0_dp
55     grid%rm3(1) =0.0_dp
56     grid%rm1(2:)= upf%r(2:)**(-1)
57     grid%rm2(2:)= upf%r(2:)**(-2)
58     grid%rm3(2:)= upf%r(2:)**(-3)
59  END IF
60  !
61END SUBROUTINE add_upf_grid
62!
63!---------------------------------------------------------------------
64SUBROUTINE set_upf_q (upf)
65  !---------------------------------------------------------------------
66  !
67  ! For USPP we set the augmentation charge as an l-dependent array in all
68  ! cases. This is already the case when upf%tpawp or upf%q_with_l are .true.
69  ! For vanderbilt US pseudos, where nqf and rinner are non zero, we do here
70  ! what otherwise would be done multiple times in many parts of the code
71  ! (such as in init_us_1, addusforce_r, bp_calc_btq, compute_qdipol)
72  ! whenever the q_l(r) were to be constructed.
73  ! For simple rrkj3 pseudos we duplicate the information contained in q(r)
74  ! for all q_l(r).
75  !
76  ! This requires a little extra memory but unifies the treatment of q_l(r)
77  ! and allows further weaking with the augmentation charge.
78  !
79  IMPLICIT NONE
80  !
81  TYPE (pseudo_upf) :: upf
82  !
83  !     Local variables
84  !
85  INTEGER :: nb, mb, ijv, ir, ilast, l, l1, l2
86  !
87  IF ( upf%tvanp .and. .not.upf%q_with_l ) THEN
88     ALLOCATE( upf%qfuncl ( upf%mesh, upf%nbeta*(upf%nbeta+1)/2, 0:upf%nqlc-1 ) )
89     upf%qfuncl  = 0.0_DP
90
91     DO nb = 1, upf%nbeta
92        DO mb = nb, upf%nbeta
93           ! ijv is the combined (nb,mb) index
94           ijv = mb * (mb-1) / 2 + nb
95           l1=upf%lll(nb) ; l2=upf%lll(mb)
96           ! copy q(r) to the l-dependent grid
97           DO l=abs(l1-l2),l1+l2,2
98              upf%qfuncl(1:upf%mesh,ijv,l) = upf%qfunc(1:upf%mesh,ijv)
99           END DO
100! adjust the inner values on the l-dependent grid if nqf and rinner are defined
101           IF ( upf%nqf > 0 ) THEN
102              DO l = abs(l1-l2),l1+l2, 2
103                 IF ( upf%rinner (l+1) > 0.0_dp) THEN
104                    DO ir = 1, upf%kkbeta
105                       if (upf%r(ir) <upf%rinner (l+1) ) ilast = ir
106                    END DO
107                    CALL setqfnew( upf%nqf,upf%qfcoef(:,l+1,nb,mb), ilast, &
108                         upf%r, l, 2, upf%qfuncl(:,ijv,l) )
109                 END IF
110              END DO
111           END IF
112        END DO
113     END DO
114  END IF
115
116END SUBROUTINE set_upf_q
117!------------------------------------------------------------------------
118SUBROUTINE setqfnew( nqf, qfcoef, mesh, r, l, n, rho )
119  !-----------------------------------------------------------------------
120  !
121  ! ... Computes the Q function from its polynomial expansion (r < rinner)
122  ! ... On input: nqf = number of polynomial coefficients
123  ! ...    qfcoef(nqf)= the coefficients defining Q
124  ! ...          mesh = number of mesh point
125  ! ...        r(mesh)= the radial mesh
126  ! ...             l = angular momentum
127  ! ...             n = additional exponent, result is multiplied by r^n
128  ! ... On output:
129  ! ...      rho(mesh)= r^n * Q(r)
130  !
131  USE upf_kinds, ONLY: dp
132  !
133  IMPLICIT NONE
134  !
135  INTEGER,  INTENT(in):: nqf, l, mesh, n
136  REAL(dp), INTENT(in) :: r(mesh), qfcoef(nqf)
137  REAL(dp), INTENT(out) :: rho(mesh)
138  !
139  INTEGER  :: ir, i
140  REAL(dp) :: rr
141  !
142  DO ir = 1, mesh
143     rr = r(ir)**2
144     rho(ir) = qfcoef(1)
145     DO i = 2, nqf
146        rho(ir) = rho(ir) + qfcoef(i)*rr**(i-1)
147     ENDDO
148     rho(ir) = rho(ir)*r(ir)**(l+n)
149  ENDDO
150  !
151  RETURN
152  !
153END SUBROUTINE setqfnew
154!
155
156!=----------------------------------------------------------------------------=!
157  END MODULE upf_to_internal
158!=----------------------------------------------------------------------------=!
159