1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
|
!
! Copyright (C) 2007 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!--------------------------------------------------------------------------
SUBROUTINE compute_potps_new(ik,v_in,v_out,xc)
!--------------------------------------------------------------------------
!
! This routine computes the local pseudopotential by smoothing the
! all_electron potential. In input it receives, the point
! ik where the cut is done.
!
! The smooth potential has the following form for r < rc
!
! V(r) = A + B * f( q * r)
!
! where f(x) = [a*a*sen(x)/x - sen(a*x)/(a*x)]/[a*a-1]
!
! V,V' and V'' are matched to the ae value in rc and
! both V'(0) and V''(0) vanish.
!
! According to Troulier and Martins this is an optimal choice for
! rapidly converging and transferable potentials.
!
use kinds, only: dp
use radial_grids, only: ndmx
use ld1inc, only: grid
IMPLICIT NONE
REAL(DP), PARAMETER :: a = .707106781186547 !=sqrt(0.5_dp) ! the parameter defining f(x)
REAL(DP) :: &
v_in(ndmx), & ! input: the potential to pseudize
v_out(ndmx),& ! output: the pseudized potential
xc(8) ! output: the coefficients of the fit
INTEGER :: &
ik ! input: the point corresponding to rc
REAL(DP) :: &
fae, & ! the value of the all-electron function
f1ae, & ! its first derivative
f2ae ! the second derivative
REAL(DP) :: &
AA, & ! the matching value
deriv, & ! auxilairy quantities
j1(ndmx,2) ! auxiliary functions
REAL(DP) :: &
deriv_7pts, deriv2_7pts
INTEGER :: &
n, & ! counter on mesh points
nc ! counter on bessel
!
! compute first and second derivative
!
fae=v_in(ik)
f1ae=deriv_7pts(v_in,ik,grid%r(ik),grid%dx)
f2ae=deriv2_7pts(v_in,ik,grid%r(ik),grid%dx)
!
! set the matching value
!
AA = f2ae*grid%r(ik)/f1ae
! find the solution
call find_matching_x(AA,xc(5))
! rescale the solution
xc(5) = xc(5)/grid%r(ik)
xc(4) = sqrt(0.5_dp)*xc(5)
! build the auxiliary functions const and f(x)
DO nc=1,2
call sph_bes(ik+1,grid%r,xc(3+nc),0,j1(1,nc))
end do
j1(1:ik+1,2) = 0.5_dp * j1(1:ik+1,2) - j1(1:ik+1,1)
j1(1:ik+1,1) = 1.0_dp
! determine A and B
deriv = 0.5_dp *(j1(ik+1,2)-j1(ik,2))/(grid%r(ik+1)-grid%r(ik)) + &
0.5_dp *(j1(ik,2)-j1(ik-1,2))/(grid%r(ik)-grid%r(ik-1))
xc(2) = f1ae / deriv
xc(1) = fae - xc(2) * j1(ik,2)
!
! define the v_out function
!
DO n=1,ik
v_out(n)=xc(1)*j1(n,1)+xc(2)*j1(n,2)
ENDDO
DO n=ik+1,grid%mesh
v_out(n)=v_in(n)
ENDDO
RETURN
CONTAINS
SUBROUTINE find_matching_x(AA,x)
real(dp) :: x, AA
real(dP) :: aa_min, aa_max, aa_test, xmin, xmax, xtest, dx
!
dx = 0.1_dp
!
xmin = 0._dp
aa_min = 3._dp - AA
IF (aa_min < 0._dp) CALL errore('compute_potps_new', &
'unable to find a solution.. try lloc=-1',1)
xmax = dx
call fz(xmax, AA, aa_max)
! bracket the first solution
do while (aa_min * aa_max >0)
xmin = xmax
aa_min = aa_max
xmax = xmin + dx
call fz(xmax, AA, aa_max)
end do
! refine by bisection
do while (abs(aa_max-aa_min) > 1.d-8)
xtest= (xmax+xmin)/2.d0
call fz(xtest, AA, aa_test)
if (aa_test * aa_max > 0 ) then
aa_max = aa_test
xmax = xtest
else
aa_min = aa_test
xmin = xtest
endif
end do
x = (xmax+xmin)/2.d0
return
END SUBROUTINE find_matching_x
!
! the auxiliary target function x*f''(x)/f'(x)-AA
SUBROUTINE fz(x,AA,aa_f)
use kinds, only: dp
real(dp) :: x, AA, aa_f
real(dp) :: fs, f1s, f2s, fac
real(dp) :: gs, g1s, g2s, xx
fac = 1._dp/(a*a-1._dp)
call ffz(x,fs,f1s,f2s)
xx = a*x
call ffz(xx,gs,g1s,g2s)
fs = fac*(a*a*fs-gs)
f1s= fac*(a*a*f1s-a*g1s)
f2s= fac*(a*a*f2s-a*a*g2s)
aa_f = x * f2s / f1s - AA
return
END SUBROUTINE
!
! the first bessel function and its derivatives
SUBROUTINE ffz(x,fs,f1s,f2s)
use kinds, only: dp
implicit none
real(dp) :: x, fs, f1s, f2s
fs = (sin(x))/x
f1s = (x*cos(x)-sin(x))/x**2
f2s = (-x*x*sin(x) - 2*x*cos(x) + 2*sin(x))/x**3
return
END SUBROUTINE ffz
END SUBROUTINE compute_potps_new
|
