    	program hidrogeno		!!      18 May 2017
	implicit real*8(a-h,o-z)

c...... v.1.1

c-------hydrogen wavefunctions

c.......with Lapack

	logical existjell
	common/readjellbk/existjell
	
c....... input files
	open(unit=9,file='hidrogeno.inp',status='old')
c....... output files
	open(unit=10,file='hidrogeno.out',status='unknown')
	
c.......read input parameters
	call readinput(llq)
	
	do 100 lq=llq,llq

c...... construction of Hamiltonian
	call constrH(lq)

c....... diagonalization
       	call diagH(lq)

c....... output 
	call printH(lq)

100	continue
	
	stop
	end

c
c**************************************************************************
c
	subroutine readinput(lq)
	implicit real*8(a-h,o-z)

c.......read input parameters
	
c.......Parameters:
c	mxpts: max number of rmesh points
      	parameter (mxpts=10000)
   
	common/bkmesh/rmin,rmax,dx,nxpts
	common/bkener/ener0,enerf,nenerpt
	common/bkcharge/iZ
	common/bkpifunc/nfunctionsp

	namelist /meshlist/nxpts,dx
	namelist /wavelist/nfunctionsp,lq,iZ
	
	data autoev/27.2113845/

!       meshlist:
!     nxpts  -- number of total spatial mesh points
!     dx     -- space integration interval (a.u.)
	read(9,meshlist)
	if (nxpts.gt.mxpts) then
		print*,'nxpts .gt. mxpts=:',mxpts,'  (reset)'
		nxpts=mxpts
	endif
      	rmin=0.0
      	rmax=nxpts*dx
      	write(10,30) rmin,rmax
30   	format(/5x,'rmin=:',f10.4,5x,'rmax=:',f10.4)

!       wavelist
!       nfunctionsp
	read(9,wavelist)
	write(10,50) nfunctionsp,lq,iZ
50	format(/5x,'number of np functions:',i4,
     +         /5x,'lq=:',i2,5x,'Z=:',i2)

	return
	end

c
c**************************************************************************
c
	subroutine constrH(lq)
	implicit real*8(a-h,o-z)
	
c-------Generates the Hamiltonian.

	parameter (mxpts=10000)
	common/bkmesh/rmin,rmax,dx,nxpts
	common/bkhamilt/d(mxpts),e(mxpts)
	common/bkpot/potential(mxpts)
	
	data half,one,two/0.5d0,1.0d0,2.0d0/
 	
c....... reading or generating potential 
	call potentialgen(lq)

c-----------building hamiltonian on the grid
	    x0 = (rmax + rmin)/two
            do 40 i=1,nxpts
	    	r = i*dx
c...............diagonal terms
                d(i)= one/dx**2 + potential(i)
c...............non-diagonal terms
                e(i)=-half/dx**2
40          continue

	return
	end
c
c**************************************************************************
c
	subroutine diagH(lq)
	implicit real*8(a-h,o-z)

c------- Matrix diagonalization.
c------- Also stores the Eigenvalues and Eigenvectors.

c........ uses DSTEVX from LAPACK

	parameter (mxpts=10000,mxfunc=500)
	common/bkmesh/rmin,rmax,dx,nxpts
	common/bkeigen/ener(mxpts),chi(mxpts,mxfunc)
	common/bk1s/psi1s(mxpts),ener1s
	common/bkpifunc/nfunctionsp
	common/bkhamilt/d(mxpts),e(mxpts)
c	d(mxpts)			! diagonal part of the Hamiltonian
c	e(mxpts)			! non-diagonal part of H
	common/bkener/ener0,enerf,nenerpt

c........ dtqli
c	dimension  wksp(mxpts),iwksp(mxpts) ! working arrays for sorting
c	dimension  v(mxpts,mxpts)	! eigenvectors

 	data rzero,one/0.0d0,1.0d0/

c.......LAPACK
	parameter (nwork = 5*mxpts)
	dimension work(nwork)
	integer iwork(nwork),ifail(mxpts)
      	character*1 JOBZ, RANGE
	
c----------matrix diagonalization
	jobz = 'V'
	il = 1
	iu = nfunctionsp
	range = 'I'

	abstol = 1.e-8

	print*,' calling dstevx'
	call DSTEVX( JOBZ, RANGE, nxpts, d, e, VL, VU, IL, IU, ABSTOL,
     $               M, ener, chi, mxpts, WORK, IWORK, IFAIL, INFO )

	if (M.gt.mxfunc) stop 'increase mxfunc parameter'
	nenerpt = M
	if (info.ne.0) then
		print*,' diagonalization failed: ',info
		print*,' ifail=:',ifail
	endif

c....... normalization
	do 200 nw=1,M
	  rsum = rzero
	  do 100 i=1,nxpts
	   rsum = rsum + dabs(chi(i,nw)**2)
100	  continue
	  rfact = one/dsqrt(rsum * dx)
	  do 200 i=1,nxpts
	   chi(i,nw) = chi(i,nw)*rfact
200	  continue
	
	if (lq.eq.0) then
	   ener1s = ener(1)
	   do 300 i=1,nxpts
		psi1s(i)=chi(i,1)
300	   continue
	endif
		  	  			
	return
     	end
c
c**************************************************************************
c
	subroutine printH(lq)
	implicit real*8(a-h,o-z)

c-------Prints the energies and selected wavefunctions.

	parameter (mxpts=10000,mxfunc=500)
	common/bkmesh/rmin,rmax,dx,nxpts
	common/bkeigen/ener(mxpts),chi(mxpts,mxfunc)
	common/bkener/ener0,enerf,nenerpt
	character*30  outputfile
	
        character*1 filec,filed,fileu
	character*1 filel

 	logical existwave

c....... open files
	ich0 = ichar('0')
	il = lq + ich0
	filel = char(il)
	inquire (file = 'waves'//filel//'/.', exist = existwave)
	if (.not.existwave) CALL system('mkdir waves'//filel)

	do 500 ie=1,nenerpt
45	   format(5x,'n =:',i5,10x,'energy =: ',f15.4,' a.u.')
 	   write(10,45) ie,ener(ie)
c 	   write(25,*) ener(ie)
	
c....... open files
	ich0 = ichar('0')
        ic = ie/100
        id = (ie - 100*ic)/10
        iu = (ie - 100*ic - 10*id)
        ic = ic + ich0
        id = id + ich0
        iu = iu + ich0
        filec = char(ic)
        filed = char(id)
        fileu = char(iu)
           
	open(unit=45,file='waves'//filel//'/dbwl.'//filec//filed//fileu,
     +       form='formatted',status='unknown')
 
	do 75 irow=1,nxpts
	    r= irow*dx
  	    write(45,333) r,chi(irow,ie)
75	continue
333	format(1pg20.8,5x,1pg26.16)

	close(45)
500	continue
	
	return
     	end

c
c**************************************************************************
c
	subroutine potentialgen(lq)
	implicit real*8(a-h,o-z)
	
c....... generates potential

	parameter(mxpts=10000)
	common/bkwell/a01,a02,x1,x2,w1,w2
	data rzero,fourth,one,onepfive,two
     +      /0.0d0,0.25d0,1.0d0,1.50d0,2.0d0/
	common/bkmesh/rmin,rmax,dx,nxpts
	common/bkpot/potential(mxpts)
	common/bkcharge/iZ

	  do 100 i=1,nxpts
	    r = i*dx
	    potential(i) = -iZ/r + lq*(lq+one)/(two * r**2 )
c	  write(11,*) r,' ',potential(i)
100	 continue

500	return
	end	
c
c**************************************************************************
c

