1!
2! Copyright (C) 2004-2012 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!---------------------------------------------------------------
10subroutine ascheqps ( nam, lam, jam, e0, mesh, ndm, grid, vpot, thresh,&
11                      y, beta, ddd, qq, nbeta, nwfx, lls, jjs, ikk, nstop )
12  !---------------------------------------------------------------
13  !
14  !  numerical integration of a generalized radial schroedinger equation,
15  !  using Numerov with outward and inward integration and matching.
16  !  Works for both norm-conserving nonlocal and US pseudopotentials
17  !  Requires in input a good estimate "e0" of the energy
18  !
19  use io_global, only : stdout
20  use kinds, only : DP
21  use radial_grids, only: radial_grid_type, series
22
23  implicit none
24
25  type(radial_grid_type), intent(in) :: grid
26  integer, intent(in) :: &
27       nam, &
28       lam, &      ! l angular momentum
29       mesh,&      ! size of radial mesh
30       ndm, &      ! maximum radial mesh
31       nbeta,&     ! number of beta function
32       nwfx, &     ! maximum number of beta functions
33       ikk(nbeta),&! for each beta the point where it become zero
34       lls(nbeta)  ! for each beta the angular momentum
35
36  real(DP), intent(in) :: &
37       jam,       & ! j angular momentum
38       vpot(mesh),& ! the local potential
39       thresh,    & ! precision of eigenvalue
40       jjs(nwfx), & ! the j angular momentum
41       beta(ndm,nwfx), &            ! the beta functions
42       ddd(nwfx,nwfx),qq(nwfx,nwfx) ! parameters for computing B_ij
43
44  real(DP), intent(inout) :: &
45       e0,      &  ! output eigenvalue
46       y(mesh)     ! the output solution
47
48  integer, intent(out) :: &
49       nstop       ! error code, used to check the behavior of the routine
50  !
51  !    the local variables
52  !
53  integer :: &
54       ndcr,  &    ! number of required nodes
55       n1, n2, &   ! counters
56       ikl         ! auxiliary variables
57  real(DP) :: &
58       work(nbeta),& ! auxiliary space
59       e,          &  ! energy
60       ddx12,      &  ! dx**2/12 used for Numerov integration
61       sqlhf,      &  ! the term for angular momentum in equation
62       ze2,        &  ! possible coulomb term aroun the origin (set 0)
63       b(0:3),     &  ! coefficients of taylor expansion of potential
64       eup,elw,    & ! actual energy interval
65       ymx,        & ! the maximum value of the function
66       fe,integ,dfe,de, &! auxiliary for numerov computation of e
67       eps,        & ! the epsilon of the delta e
68       yln, xp, expn,& ! used to compute the tail of the solution
69       int_0_inf_dr  ! integral function
70
71  real(DP), allocatable :: &
72       fun(:),  &   ! integrand function
73       f(:),    &   ! the f function
74       el(:),c(:) ! auxiliary for inward integration
75
76  integer, parameter :: &
77       maxter=100    ! maximum number of iterations
78
79  integer :: &
80       n,  &    ! counter on mesh points
81       iter,&   ! counter on iteration
82       ik,  &   ! matching point
83       ns,  &   ! counter on beta functions
84       l1,  &   ! lam+1
85       nst, &   ! used in the integration routine
86       ierr, &
87       ncross,& ! actual number of nodes
88       nstart  ! starting point for inward integration
89
90  logical, save :: first(0:10,0:10) = .true.
91
92
93   if (mesh.ne.grid%mesh) &
94        call errore('compute_solution','mesh dimension is not as expected',1)
95  !
96  !  set up constants and allocate variables the
97  !
98  allocate(fun(mesh), stat=ierr)
99  allocate(f(mesh), stat=ierr)
100  allocate(el(mesh), stat=ierr)
101  allocate(c(mesh), stat=ierr)
102  nstop=0
103  nstart=0
104  e=e0
105!  write(6,*) 'entering ', nam,lam, e
106  eup=0.3_DP*e
107  elw=1.3_dp*e
108  ndcr=nam-lam-1
109!  write(6,*) 'entering ascheqps', vpot(mesh-20)*grid%r(mesh-20)
110
111  ddx12=grid%dx*grid%dx/12.0_dp
112  l1=lam+1
113  nst=l1*2
114  sqlhf=(DBLE(lam)+0.5_dp)**2
115  !
116  !  series developement of the potential near the origin
117  !
118  do n=1,4
119     y(n)=vpot(n)
120  enddo
121  call series(y,grid%r,grid%r2,b)
122
123  !      write(stdout,*) 'enter lam,eup,elw,e',lam,nbeta,eup,elw,e
124  !
125  !  set up the f-function and determine the position of its last
126  !  change of sign
127  !  f < 0 (approximatively) means classically allowed   region
128  !  f > 0         "           "        "      forbidden   "
129  !
130  do iter=1,maxter
131!  write(6,*) 'starting iter', iter, elw, e, eup
132  ik=1
133  f(1)=ddx12*(grid%r2(1)*(vpot(1)-e)+sqlhf)
134  do n=2,mesh
135     f(n)=ddx12*(grid%r2(n)*(vpot(n)-e)+sqlhf)
136     if( f(n).ne.sign(f(n),f(n-1)).and.n.lt.mesh-5 ) ik=n
137  enddo
138!  if (ik.eq.0.and.nbeta.eq.0) ik=mesh*3/4
139  if (ik.eq.1.or.grid%r(ik)>4.0_DP) ik=mesh*3/4
140
141  if(ik.ge.mesh-2) then
142     do n=1,mesh
143        write(stdout,*) grid%r(n), vpot(n), f(n)
144     enddo
145     call errore('compute_solution', 'No point found for matching',1)
146  endif
147  !
148  !     determine if ik is sufficiently large
149  !
150  do ns=1,nbeta
151     if (lls(ns).eq.lam .and. jjs(ns).eq.jam .and. ikk(ns).gt.ik) ik=ikk(ns)+3
152  enddo
153  !
154  !     if everything is ok continue the integration and define f
155  !
156  do n=1,mesh
157     f(n)=1.0_dp-f(n)
158  enddo
159  !
160  !  determination of the wave-function in the first two points by
161  !  series developement
162  !
163  ! no coulomb divergence in the origin for a pseudopotential
164  ze2=0.d0
165  call start_scheq( lam, e, b, grid, ze2, y )
166  !
167  !    outward integration before ik
168  !
169  call integrate_outward (lam,jam,e,mesh,ndm,grid,f,b,y,beta,ddd,qq,&
170       nbeta,nwfx,lls,jjs,ikk,ik)
171
172  ncross=0
173  ymx=0.0_dp
174  do n=2,ik-1
175     if ( y(n) .ne. sign(y(n),y(n+1)) ) ncross=ncross+1
176     ymx=max(ymx,abs(y(n+1)))
177  end do
178!
179!  If at this point the number of nodes is wrong it means that something
180!  is probably wrong in the calling routines. A ghost might be present
181!  in the pseudopotential. With a nonlocal pseudopotential there is no
182!  node theorem so strictly speaking the following instructions are
183!  wrong but sometimes they help so we keep them here.
184!
185  if( ndcr /= ncross .and. first(nam,lam) ) then
186      write(stdout,"(/,7x,'Warning: n=',i1,', l=',i1,' expected ',i1,' nodes,',&
187        & ' found ',i1)") nam,lam,ndcr,ncross
188      write(stdout,'(7x,a,/,7x,a,/)') 'Setting wfc to zero for this iteration',&
189                    '(This warning will only be printed once per wavefunction)'
190      first(nam,lam) = .false.
191  endif
192
193  if(ndcr < ncross) then
194     !
195     !  too many crossings. e is an upper bound to the true eigenvalue.
196     !  increase abs(e)
197     !
198     eup=e
199     e=0.9_dp*elw+0.1_dp*eup
200!     write(6,*) 'too many crossing', ncross, ndcr
201!     call errore('aschqps','wrong number of nodes. Probably a Ghost?',1)
202     y=0.0_DP
203     ymx=0.0_dp
204     go to 300
205  else if (ndcr > ncross) then
206     !
207     !  too few crossings. e is a lower bound to the true eigenvalue.
208     !  decrease abs(e)
209     !
210     elw=e
211     e=0.9_dp*eup+0.1_dp*elw
212!     write(6,*) 'too few crossing', ncross, ndcr
213!     call errore('aschqps','wrong number of nodes. Probably a Ghost?',1)
214     y=0.0_DP
215     ymx=0.0_dp
216     go to 300
217  end if
218  !
219  !    inward integration up to ik
220  !
221  call integrate_inward(e,mesh,ndm,grid,f,y,c,el,ik,nstart)
222  !
223  !  if necessary, improve the trial eigenvalue by the cooley's
224  !  procedure. jw cooley math of comp 15,363(1961)
225  !
226  fe=(12.0_dp-10.0_dp*f(ik))*y(ik)-f(ik-1)*y(ik-1)-f(ik+1)*y(ik+1)
227  !
228  ! audjust the normalization if needed
229  !
230  if(ymx.ge.1.0e10_dp) y=y/ymx
231  !
232  !  calculate the normalization
233  !
234  do n1=1,nbeta
235     if (lam.eq.lls(n1).and.abs(jam-jjs(n1)).lt.1.e-7_dp) then
236        ikl=ikk(n1)
237        do n=1,ikl
238           fun(n)=beta(n,n1)*y(n)*grid%sqr(n)
239        enddo
240        work(n1)=int_0_inf_dr(fun,grid,ikl,nst)
241     else
242        work(n1)=0.0_dp
243     endif
244  enddo
245  do n=1,nstart
246     fun(n)= y(n)*y(n)*grid%r(n)
247  enddo
248  integ=int_0_inf_dr(fun,grid,nstart,nst)
249  do n1=1,nbeta
250     do n2=1,nbeta
251        integ = integ + qq(n1,n2)*work(n1)*work(n2)
252     enddo
253  enddo
254  dfe=-y(ik)*f(ik)/grid%dx/integ
255  de=-fe*dfe
256  eps=abs(de/e)
257!  write(6,'(i5, 3f20.12)') iter, e, de
258  if(abs(de).lt.thresh) go to 600
259  if(eps.gt.0.25_dp) de=0.25_dp*de/eps
260  if(de.gt.0.0_dp) elw=e
261  if(de.lt.0.0_dp) eup=e
262  e=e+de
263  if(e.gt.eup) e=0.9_dp*eup+0.1_dp*elw
264  if(e.lt.elw) e=0.9_dp*elw+0.1_dp*eup
265300 continue
266  enddo
267  nstop=1
268  if (nstart==0) goto 900
269600 continue
270  !
271  !   exponential tail of the solution if it was not computed
272  !
273  if (nstart.lt.mesh) then
274     do n=nstart,mesh-1
275        if (y(n) == 0.0_dp) then
276           y(n+1)=0.0_dp
277        else
278           yln=log(abs(y(n)))
279           xp=-sqrt(12.0_dp*abs(1.0_dp-f(n)))
280           expn=yln+xp
281           if (expn.lt.-80.0_dp) then
282              y(n+1)=0.0_dp
283           else
284              y(n+1)=sign(exp(expn),y(n))
285           endif
286        endif
287     enddo
288  endif
289  !
290  !  normalize the eigenfunction as if they were norm conserving.
291  !  If this is a US PP the correct normalization is done outside this
292  !  routine.
293  !
294  do n=1,mesh
295     el(n)=grid%r(n)*y(n)*y(n)
296  enddo
297  integ=int_0_inf_dr(el,grid,mesh,nst)
298  if (integ>0.0_DP) then
299     integ=sqrt(integ)
300     do n=1,mesh
301        y(n)=grid%sqr(n)*y(n)/integ
302     enddo
303     e0=e
304  else
305     nstop=1
306  endif
307
308900 continue
309
310  deallocate(el)
311  deallocate(f )
312  deallocate(c )
313  deallocate(fun )
314  return
315
316end subroutine ascheqps
317