*---------------------------------------------------------------
      program helium_gauss
*---------------------------------------------------------------
*
*	Atomo di elio mediante espansione su base gaussiana
*	e diagonalizzazione con metodi di algebra lineare.
*
*	Compilazione:
*	g77 helium_gauss.f ~/lib/lapack.a ~/lib/blas.a

      implicit none

      integer nalpx, nnalpx 
      real*8  e2,pi,tpi
      parameter (nalpx=40, nnalpx=nalpx*(nalpx+1),e2=2.d0, 
     +           pi=3.14159265358979d0,tpi =2.d0*pi)
      integer ia,ib,ic,id,n_alpha,nn
      real*8 zeta, lambda_min, lambda_max, lambda, vee
      real*8 alpha(nalpx), aa, e(nalpx), ee(nalpx*nalpx)
      real*8 s(nalpx,nalpx,0:1), h(nalpx,nalpx,0:1), v(nalpx,nalpx)
      real*8 ss(nnalpx,nnalpx), hh(nnalpx,nnalpx), vv(nnalpx,nnalpx)
      real*8 r,r1,r2,f,q,dr, prob(0:200), prob2(0:200,0:200)
      integer i,j, i1, i2, nrx
*
*	Dati di input
*
      zeta = 2.0d0
      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 '

*
*	Definizione parametri delle gaussiane
*	(lambda=lunghezze caratteristiche in u.a.)
*
      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

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

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

********************************************************************
*	2 elettroni, interazioni S-S
********************************************************************
      write(*,'(/,A)') 
     $    "Problema elio: 2 elettroni, interazioni solo S-S"

*
*	Vi sono nn funzioni nella base di funzioni a due elettroni
*	(usata per applicare il metodo variazionale).
*	Ciascuna di esse e' ottenuta combinando in modo simmetrico
*	una funzione della base a un elettrone con un'altra 
*	funzione della base a un elettrone.
*	Nel loop sotto gli indici i e j varieranno tra 1 e nn.
*
      nn = n_alpha*(n_alpha+1)/2
      i = 0
      do ia=1,n_alpha
      do ib=ia,n_alpha
      i = i+1
      j = 0
      do ic=1,n_alpha
      do id=ic,n_alpha
      j = j+1
         ss(i,j) = (s(ia,ic,0)*s(ib,id,0)+s(ia,id,0)*s(ib,ic,0))
         hh(i,j) = (h(ia,ic,0)*s(ib,id,0)+s(ia,id,0)*h(ib,ic,0)+
     +              s(ia,ic,0)*h(ib,id,0)+h(ia,id,0)*s(ib,ic,0))
         hh(i,j) = hh(i,j) + vee(alpha,n_alpha,ia,ic,ib,id,"ss") +
     +                       vee(alpha,n_alpha,ia,id,ib,ic,"ss") 
      end do
      end do
      end do
      end do
*
*	Risoluzione:
*
      call diaghg(nn,hh,ss,nnalpx,ee,vv)
      write (*,*) 'Primi tre autovalori:', (ee(ia),ia=1,3)
*
*	Calcolo densita' di probabilita' prob(r) e prob2(r1,r2)
*
      dr =0.1d0/zeta
      nrx = 60
      q = 0.d0
      do i1=0, nrx
         r1 = dr*i1
         prob(i1) = 0.d0
         do i2 =0,nrx
            r2 = dr*i2
            f = 0.d0
            i = 0
            do ia=1,n_alpha
            do ib=ia,n_alpha
               i=i+1
               f = f + exp(-alpha(ia)*r1*r1-alpha(ib)*r2*r2)*vv(i,1)
               f = f + exp(-alpha(ib)*r1*r1-alpha(ia)*r2*r2)*vv(i,1)
            end do
            end do
            f = f/dsqrt(2.d0)
            prob2(i1,i2) = (4.d0*pi*r1*r2*f)**2
            prob(i1) = prob(i1) + prob2(i1,i2)*dr
         end do
         q   = q + prob(i1)*dr
      end do
c      write (*,*) q
*
*	Output su file
*
      open (7,file='gs-wfc.out',status='unknown',form='formatted')
      do i1=0, nrx
         r1 = dr*i1
         do i2 =0,nrx
            r2 = dr*i2
            write(7,*) r1, r2, prob2(i1,i2),  prob(i1)*prob(i2), 
     +                         prob2(i1,i2) - prob(i1)*prob(i2)
         end do
         write(7,*) 
      end do
      close(7)

********************************************************************
*	2 elettroni, interazioni S-S, S-P, P-P
********************************************************************
      write(*,'(/,A)') 
     $    "Problema elio: 2 elettroni, interazioni S-S,S-P,P-P"

*
*	Aggiunta dei nuovi integrali
*
      i = 0
      do ia=1,n_alpha
      do ib=ia,n_alpha
      i = i+1
      j = 0
      do ic=1,n_alpha
      do id=ic,n_alpha
      j = j+1
         ss(i+nn,j+nn) =3*(s(ia,ic,1)*s(ib,id,1)+s(ia,id,1)*s(ib,ic,1))
         hh(i+nn,j+nn) =3*(h(ia,ic,1)*s(ib,id,1)+s(ia,id,1)*h(ib,ic,1)+
     +                     s(ia,ic,1)*h(ib,id,1)+h(ia,id,1)*s(ib,ic,1))
         hh(i+nn,j+nn) = hh(i+nn,j+nn) + 
     +                   (vee(alpha,n_alpha,ia,ic,ib,id,"pp") +
     +                    vee(alpha,n_alpha,ia,id,ib,ic,"pp") ) 
         ss(j+nn,i) = 0.d0
         hh(i,j+nn) = (vee(alpha,n_alpha,ia,ic,ib,id,"sp") +
     +                 vee(alpha,n_alpha,ia,id,ib,ic,"sp") ) 
         ss(j+nn,i) = ss(i,j+nn)
         hh(j+nn,i) = hh(i,j+nn)
      end do
      end do
      end do
      end do
      nn = 2 * nn
*
*	Risoluzione:
*
      call diaghg(nn,hh,ss,nnalpx,ee,vv)
      write (*,*) 'Primi tre autovalori:', (ee(ia),ia=1,3)

      end


*---------------------------------------------------------------------
      function vee(alpha,n_alpha,ia,ib,ic,id,type)
*---------------------------------------------------------------------
*
*	Calcolo termine di repulsione coulombiana fra le distribuzioni
*	di carica dei due elettroni negli integrali H_ij
*
      implicit none 
      real*8 pi, e2
      parameter (pi=3.14159265358979d0, e2=2.d0)
      integer ia,ib,ic,id,n_alpha
      real*8 vee, alpha(n_alpha), s, r
      character*2 type

      s = alpha(ia)+alpha(ib)
      r = alpha(ic)+alpha(id)

      if (type.eq."ss") then
*		caso S-S
         vee= 2.d0*e2*pi**2.5d0/(s*r)/sqrt(s+r)
      else if (type.eq."sp") then
*		caso S-P
         vee= e2*pi**2.5d0/(s*r)/(s+r)**1.5d0
      else if (type.eq."pp") then
*		caso P-P
         vee= e2*pi**2.5d0/(s*r)**2*(1.5d0*s*r+(s+r)**2)/(s+r)**2.5d0
      else 
         write (*,*) "Integrale di tipo sconosciuto: ",type
         stop
      end if

      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
