1!
2! Copyright (C) 2010 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!---------------------------------------------------------------
9subroutine start_scheq(lam,e,b,grid,ze2,solution)
10  !---------------------------------------------------------------
11  !
12  !  determines the wave-function in the first two points by
13  !  series developement. It receives as input:
14  !  lam the angular momentum
15  !  e   the energy in Ry
16  !  b(0:3) the coefficients of a polynomial that interpolates the
17  !         potential in the first three points
18  !  grid  the mesh
19  !  ze2   the zed of the mesh
20  !  in output solution(1:2) contains the solution in the first two points
21  !
22USE kinds, ONLY : DP
23USE radial_grids, ONLY : radial_grid_type
24IMPLICIT NONE
25TYPE(radial_grid_type), INTENT(IN) :: grid
26INTEGER, INTENT(IN) :: lam
27REAL(DP), INTENT(IN) :: b(0:3), e
28REAL(DP) :: ze2, xl1, x4l6, x6l12, x8l20, b0e, c1, c2, c3, c4, rr1, rr2
29REAL(DP) :: solution(1:2)
30INTEGER :: l1
31!
32!  set up constants and initialize
33!
34l1=lam+1
35xl1=lam+1.0_DP
36x4l6=4.0_dp*lam+6.0_dp
37x6l12=6.0_dp*lam+12.0_dp
38x8l20=8.0_dp*lam+20.0_dp
39!
40!
41b0e=b(0)-e
42c1=0.5_dp*ze2/xl1
43c2=(c1*ze2+b0e)/x4l6
44c3=(c2*ze2+c1*b0e+b(1))/x6l12
45c4=(c3*ze2+c2*b0e+c1*b(1)+b(2))/x8l20
46rr1=(1.0_dp+grid%r(1)*(c1+grid%r(1)*(c2+grid%r(1)*(c3+grid%r(1)*c4))))*grid%r(1)**l1
47rr2=(1.0_dp+grid%r(2)*(c1+grid%r(2)*(c2+grid%r(2)*(c3+grid%r(2)*c4))))*grid%r(2)**l1
48solution(1)=rr1/grid%sqr(1)
49solution(2)=rr2/grid%sqr(2)
50
51return
52end subroutine start_scheq
53