*---------------------------------------------------------------
      program hydrogen_gauss
*---------------------------------------------------------------
*
*	Atomo di idrogeno mediante espansione su base gaussiana
*	e diagonalizzazione con metodi di algebra lineare.
*
*	Compilazione:
*	g77 hydrogen_gauss.f -lblas -llapack

      implicit none

      integer nalpx
      real*8  e2,pi,tpi
      parameter (nalpx=100, e2=2.d0, pi=3.14159265358979d0,tpi=2.d0*pi)
      integer ia,ib,n_alpha
      real*8 zeta, lambda_min, lambda_max, lambda
      real*8 alpha(nalpx), aa, e(nalpx)
      real*8 s(nalpx,nalpx), h(nalpx,nalpx), v(nalpx,nalpx)
      real*8 r,f,prob,q,dr
      integer i,j, nrx
*
*	Dati di input
*
      zeta = 1.d0
      write (*,'(a,f10.6)') ' Carica atomica = ', zeta
      write (*,'($,a)') ' n_alpha, lambda_min, lambda_max >> '
      read(*,*) n_alpha, lambda_min, lambda_max
      if (n_alpha.gt.nalpx) stop 'n_alpha.gt.nalpx '
      if (n_alpha.lt.1) stop 'n_alpha.lt.1 '

*
*	Definizione parametri delle gaussiane
*	(lambda=lunghezze caratteristiche in u.a.)
*

      if (n_alpha.gt.1) then
         do ia=1,n_alpha
             lambda = (lambda_min*(ia-1) + lambda_max*(n_alpha-ia))/
     $                (n_alpha-1)
             alpha(ia) =  (zeta/lambda)**2
         end do
      else 
         alpha(1) = (zeta/lambda_min)**2
      end if

********************************************************************
*   Stati S
********************************************************************
      write(*,'(/,A)') "Problema idrogenoide: 1 elettrone, stati S"

*
*	Assegnazione valori integrali di overlap e elementi
*	di matrice dell'hamiltoniano sulla base S:
*
      do ia=1,n_alpha
         do ib=1,n_alpha
            aa = alpha(ia)+alpha(ib)
            s(ia,ib) = (pi/aa)**1.5d0
            h(ia,ib) = s(ia,ib)*6.d0*alpha(ia)*alpha(ib)/aa -
     +                 e2*zeta*tpi/aa 
         end do
      end do

*
*	Risoluzione:
*
      call diaghg(n_alpha,h,s,nalpx,e,v)

      write (*,*) 'Primi tre autov.:', (e(ia),ia=1,min(3,n_alpha) )
*
*	I coefficienti dell'espansione sono immagazzinati
*	in v(j,i)  [j=componente j-ima funzione di base,
*	, i=indice autovalore]
*
*
*	Output coefficienti stato a energia piu' bassa:
*
      open (7,file='s-coeff.out',status='unknown',form='formatted')
      do j=1,n_alpha
         write(7,*) j, zeta/sqrt(alpha(j)), v(j,1)
      end do
      close(7)
*
*	Output stato a energia piu' bassa
*
      open (7,file='s-wfc.out',status='unknown',form='formatted')
      dr =0.1d0/zeta
      nrx = 200
      q = 0.d0
      do i=0, nrx
         r = dr*i
         f = 0.d0
         do j=1,n_alpha
            f = f + exp(-alpha(j)*r*r)*v(j,1)
         end do
         prob = 4.d0*pi*(r*f)**2
         q   = q + prob*dr
         write(7,*) r, f, r*f, prob
      end do 
*	verifica normalizzazione (se desiderata):
c      write (*,*) q
      close(7)

********************************************************************
*   Stati P
********************************************************************
      write(*,'(/,A)') "Problema idrogenoide: 1 elettrone, stati P"

*
*	Assegnazione valori integrali di overlap e elementi
*	di matrice dell'hamiltoniano sulla base P:
*
      do ia=1,n_alpha
         do ib=1,n_alpha
            aa = alpha(ia)+alpha(ib)
            s(ia,ib) = 0.5d0/aa*(pi/aa)**1.5d0
            h(ia,ib) = s(ia,ib)*10.d0*alpha(ia)*alpha(ib)/aa -
     +                 e2*zeta*tpi/aa/aa/3.d0
         end do
      end do

*
*	Risoluzione:
*
      call diaghg(n_alpha,h,s,nalpx,e,v)
      write (*,*) 'Primi tre autov.:', (e(ia),ia=1,min(3,n_alpha) )

*
*	Output coefficienti stato a energia piu' bassa:
*
      open (7,file='p-coeff.out',status='unknown',form='formatted')
      do j=1,n_alpha
         write(7,*) j, zeta/sqrt(alpha(j)), v(j,1)
      end do
      close(7)
*
*	Output stato a energia piu' bassa:
*
      open (7,file='p-wfc.out',status='unknown',form='formatted')
      dr =0.1d0/zeta
      nrx = 200
      q = 0.d0
      do i=0, nrx
         r = dr*i
         f = 0.d0
         do j=1,n_alpha
            f = f + r*exp(-alpha(j)*r*r)*v(j,1)
         end do
         prob = 4.d0/3.0d0*pi*(r*f)**2
         q   = q + prob*dr
         write(7,*) r, f, r*f, prob
      end do 
*	verifica normalizzazione (se desiderata):
c     write (*,*) q
      close(7)

      end


*-----------------------------------------------------------------------
      subroutine diaghg (n,h,s,ldh,e,v)
*     =================
*
*	Calcolo autovalori e autovettori del problema generalizzato
*	Hv=eSv, con H matrice hermitiana, S matrice di overlap.
*	Lascia le matrici S e H invariate dopo l'esecuzione.
*
*-----------------------------------------------------------------------

      implicit none

* on INPUT
      integer    n,          ! dimension of the matrix to be diagonalized
     *           ldh         ! leading dimension of h,s and v, as declared
                             ! in the calling pgm unit
      real*8 h(ldh,n)        ! matrix to be diagonalized
      real*8 s(ldh,n)        ! overlap matrix

* on OUTPUT
      real*8     e(n)        ! eigenvalues
      real*8 v(ldh,n)        ! eigenvectors (column-wise)

* LOCAL variables
      integer 
     +          lwork,       ! dimension of lapack routine ZHEGV workspace
     +          info,        ! error flag
     +          i,j,nn       ! local counter

      real*8 small
      parameter (small=1.d-10)
      external dcopy, ddot
      real*8 ddot, rwork
      real*8 work(3*n), aux(n,n), uno(n,n)

      info = 0
      lwork=3*n
*
*	Copia S in matrice ausiliaria perche' dsyev distrugge
*	la matrice
*
      do i=1,n
         do j=1,n
            aux(i,j)=s(i,j)
         end do
      enddo
*
*	Diagonalizza S
*
      call dsyev('V','U',n,aux,n,e,work,lwork,info)
      if (info.ne.0) stop 'S-matrix diagonalization failed '
*
*	Tiene solo le combinazioni linearmente indipendenti
*	(dal punto di vista numerico)
*
      nn = 0
      do i=1,n
         if (e(i).gt.small) then
            nn = nn + 1
            call dcopy(n,aux(1,i),1,aux(1,nn),1)
            call dscal(n,1.d0/dsqrt(e(i)),aux(1,nn),1)
         end if
      end do
      write(*,*) " # di vettori linearmente indipendenti =", nn, n

*
*	Trasforma H
*
      do i=1,n
      do j=1,nn
         v(i,j) = ddot(n,h(i,1),ldh,aux(1,j),1)
      end do
      end do
      do i=1,nn
      do j=1,nn
         uno(i,j) = ddot(n,aux(1,i),1,v(1,j),1)
      end do
      end do

*
*	Diagonalizza H
*
      info = 0
      call dsyev('V','U',nn,uno,n,e,work,lwork,info)
      do i=1,n
      do j=1,nn
         v(i,j) = ddot(nn,aux(i,1),n,uno(1,j),1)
      end do
      end do
      
      if (info.ne.0) stop 'H-matrix diagonalization failed '
        
      end

