      program elioscf                       !inizio programma principale: nome
      parameter (nmax=500)                !serve nelle dichiarazioni successive
      dimension r(nmax),chi(nmax),q(nmax),aux(nmax),vh(nmax)  !dimensioni array
      dimension vtot(nmax),up(nmax),upp(nmax),cf(nmax)
      character*40 messag           !stringa alfanumerica per messaggi d'errore

      open(unit=70,status='unknown',file='funzion.txt')  !70 = file griglia.out
      open(unit=80,status='unknown',file='energie.txt')  !80 = file energia.out
      z=2                                             !numero atomico dell'elio
      r1=0.01/z                                !parametri della griglia radiale
      grigl=1.05
      al=alog(grigl)
      nx=alog(z*75./0.005)/al
      r0=r1/sqrt(grigl)                    !raggio per primo pezzetto integrali
      if(nx.gt.nmax) then                   !controllo dimensioni: ho sfondato?
         write(70,3000) nx,nmax                                  !scrivo avviso
 3000    format(/1x,'nx = ',i5,3x,'nmax = ',i5)                !formato avviso
         messag='griglia con troppi punti, aumentare nmax'
         call err(messag)
         stop                                              !stop se ho sfondato
      endif                                                      !fine ciclo if
      call griglia(r1,grigl,nx,nmax,r)         !costruisco griglia radiale r(i)

c scegli Z* ottimale = Z-5/16 e calcola risultati variazionali
         zstar=z-(float(5)/16.)

            do i=1,nx                         !definisco chi(i), di simmetria s
            chi(i)=2.*(zstar**1.5)*r(i)*exp(-zstar*r(i))
            end do
         coef=2.*(zstar**1.5) !coef: derivata 1a radiale (analit.) di chi a r=0
      call poisson
     &(r,al,z,r0,coef,chi,q,aux,auxx,vh,vtot,hartree,nx,nmax)
      etotal=-(zstar**2)+hartree+2.*(zstar-z)*auxx    !energia totale (variaz.)
      write(80,1000) 0, etotal                             !scrivo su unita' 80
 1000 format(1x,i5,5x,1p4e15.4)                           !formato di scrittura

c comincia ciclo autoconsistente partendo da  fz. variazionale (full feedback)
      n=1                                           !numero quantico principale
      l=0                                             !numero quantico angolare
      eps=(-zstar**2)/2         !valore iniziale dell'autovalore da cui partire

         do i=1,100                 !ciclo iterativo, max numero iterazioni=100

         relerr=1.e-5   !precisione relativa (criterio accuratezza) per schroed
         epsold=eps
         call schroed(n,l,eps,relerr,mch,z,r,vtot,chi,up,upp,cf,nx,nmax)

c coef = derivata prima radiale (approx. numerica) di chi in r=0
         coef=((chi(1)*r(2)/r(1))-(chi(2)*r(1)/r(2)))/(r(2)-r(1))
         deltae=abs(eps-epsold)        !cambio energia da un'iteraz. all'altra
         if(deltae.lt.1e-8) go to 100 !esci da ciclo se deltae < 0.00000001 au

         call poisson
     &(r,al,z,r0,coef,chi,q,aux,auxx,vh,vtot,hartree,nx,nmax)

      etotal=2*eps-hartree
      write(80,1000) i, etotal                            !scrivo su unita' 80

         end do                                          !fine ciclo iterativo
c se arriva qui vuol dire che deltar non e' minore di 0.00000001 au: avvertire
      write(70,4000) i,deltae
 4000 format(/1x,'completate ', i3,' iterazioni; 
     &             delta autovalore = ',1pe12.5,' a.u.'/)

  100 continue
      iteraz=i

      write(6,5000) etotal,iteraz
 5000 format(/1x,'energia totale = ',1pe15.4,2x,' iterazioni = ',i3/)
      write(70,1000) (i,r(i),chi(i),q(i),vh(i),i=1,nx)        !scrivo unita' 70
      stop                                           !stop quando si arriva qui
      end                                                !fine modulo programma

      subroutine poisson
     &(r,al,z,r0,coef,chi,q,aux,auxx,vh,vtot,hartree,nx,nmax)
c input: funzione d'onda radiale chi e suo coefficiente lineare in r=0
c output: potenziale di hartree, potenziale totale ed energia di hartree

      dimension r(nmax),chi(nmax),q(nmax),aux(nmax),vh(nmax),vtot(nmax)

         an=(coef**2)*(r0**3)/3.     !pezzetto da r=0 al primo punto di griglia
            do i=1,nx                                    !ciclo: carica radiale
            an=an+(chi(i)*chi(i)*al*r(i))
            q(i)=an
            end do                                   !fine ciclo carica radiale

         an=(coef**2)*(r0**2)/2.     !pezzetto da r=0 al primo punto di griglia
            do i=1,nx              !ciclo: aux = integrale in dr di (chi*chi/r)
            an=an+(chi(i)*chi(i)*al)
            aux(i)=an 
            end do                                              !fine ciclo aux

         auxx=aux(nx)
            do i=1,nx                !ciclo: sottraggo auxx (cost. d'integraz.)
            aux(i)=aux(i)-auxx                  !per ottenere vh=1/r a r grande
            end do                                              !fine ciclo aux
            do i=1,nx                          !ciclo: potenziale di hartree vh
            vh(i)=(q(i)/r(i))-aux(i)
            vtot(i)=-(z/r(i))+vh(i)
            end do                               !fine ciclo potenziale hartree

         hartree=0.0             !qui per il primo pezzetto va bene hartree = 0
            do i=1,nx                                !ciclo: energia di hartree
            hartree=hartree+vh(i)*chi(i)*chi(i)*al*r(i)      !basate su fz. chi
            end do                                       !fine hartree e totale

      return
      end

      subroutine schroed(n,l,eps,relerr,mch,z,r,v,u,up,upp,cf,mmax,nmax)
c integra schroedinger radiale su griglia logaritmica
c utilizzando formule d'estrapolazione ed interpolazione di adams 
c per integrazione outward e inward (abramowitz & stegun p. 896)

      dimension r(nmax),v(nmax),u(nmax),up(nmax)
      dimension upp(nmax),cf(nmax)
      character*40 messag           !stringa alfanumerica per messaggi d'errore
      real kap

c*** FASE A: preliminari************************************************
c***********************************************************************

c parametri griglia, errore relativo energia, l(l+1), energia max e min
      amesh=r(2)/r(1)
      al=alog(amesh)
      als=al**2
      sls=float(l*(l+1))
      emax=v(mmax)+0.5*sls/r(mmax)**2
      emin=-0.0
         do i=1,mmax
         emin=amin1(emin,v(i)+0.5*sls/r(i)**2)
         end do

c riporta l'energia (negativa) dell'input nella finestra emin-emax
      if(eps .gt. emax) eps=1.25*emax
      if(eps .lt. emin) eps=0.75*emin
      if(eps .gt. emax) eps=0.5*(emax+emin)

c azzera il contatore delle iterazioni
      iteraz=0

c***********************************************************************
c*** FINE FASE A *******************************************************

c*** FASE B: iterazioni successive che migliorano eps rispetto al valore
c*** di input fino alla precisione voluta*******************************
c***********************************************************************

   10 iteraz=iteraz+1
c iteraz conta iterazioni (era nint: ribattezzato in onore di Marchetti)

c stop (con messaggio) se si superano 100 iterazioni
      if(iteraz .gt. 100) then
      messag='schroed, iteraz .gt. 100                  '
      call err(messag)
      endif

c definisci array coefficienti nell'eq. differ. per u (cioe' per chi)
      do i=1,mmax
      cf(i)=als*sls + 2.0*als*(v(i)-eps)*r(i)**2
      end do

c trova il punto d'inversione classica sapendo l'array cf e l'energia eps,
c e' il punto di raccordo della fz. d'onda radiale (cioe' l'indice mch)
      do ii=1,mmax-2
      i=mmax-ii
      if(cf(i-1) .gt. 0.0) go to 30
      if(cf(i) .le. 0.0) go to 30
      mch=i
      go to 40
   30 continue
      end do

c stampa messaggio d'errore se non si trova il punto
      messag='schroed: non si trova punto inversione classica   '
      call err(messag)
   40 continue

c>>> FASE B1: integrazione da r=0 fino al punto di raccordo r(mch)>>>>>>

c funzione d'onda a piccoli r (i=1,4): serie di taylor
      c2=-z/float(l+1)
      c3=z**2/((2*l+3)*(l+1)) - eps/(2*l+3)
         do i=1,4
         u(i)=r(i)**(l+1) + c2*r(i)**(l+2) + c3*r(i)**(l+3)
         up(i)=al*((l+1)*r(i)**(l+1) + c2*(l+2)*r(i)**(l+2)
     &    + c3*(l+3)*r(i)**(l+3))
         upp(i)=al*up(i)+cf(i)*u(i)
         end do

c integrazione verso l'esterno, cioe' da r=r(5) fino al raggio r(mch)
c usando predictor una volta, corrector due volte
      node=0
         do i=4,mch-1
         u(i+1)=u(i)+aeo(up,i,nmax)
         up(i+1)=up(i)+aeo(upp,i,nmax)
            do it=1,2
            upp(i+1)=al*up(i+1)+cf(i+1)*u(i+1)
            up(i+1)=up(i)+aio(upp,i,nmax)
            u(i+1)=u(i)+aio(up,i,nmax)
            end do
         if(u(i+1)*u(i) .le. 0.0) node=node+1
         end do
      if(node-n+l+1) 80,100,90

c troppo pochi nodi: inutile andare avanti, provare subito nuova energia
   80 emin=eps
      eps=0.75*eps
      if(eps .gt. emax) eps=0.5*(emin+emax)
      go to 10

c troppi nodi: inutile andare avanti, provare subito nuova energia
   90 emax=eps
      eps=1.25*eps
      if(eps .lt. emin) eps=0.5*(emin+emax)
      go to 10

c numero di nodi giusto: ok, procedere all'integrazione verso l'interno
  100 uout=u(mch)
      upout=up(mch)

c>>> FINE FASE B1 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

c>>> FASE B2: integrazione da "infinito" al punto di raccordo r(mch)>>>>

c ultimo punto r(nin): 10 volte raggio del punto d'inversione classico
      nin=mch+2.3/al
      if(nin+4 .gt. mmax) nin=mmax-4

c funzione d'onda a grandi r (ultimi 4 punti): esponenziale semplice
      kap=sqrt((sls/r(nin)**2) + 2.0*(v(nin)-eps))
         do i=nin,nin+4
         u(i)=exp(-kap*(r(i)-r(nin)))
         up(i)=-r(i)*al*kap*u(i)
         upp(i)=al*up(i)+cf(i)*u(i)
         end do

c integrazione verso l'interno, cioe' da un raggio grande r(nin) verso
c raggi piu' piccoli, fino al raggio di raccordo r(mch)
      i=nin+1
  120 i=i-1
      if(i .le. mch) go to 140
      u(i-1)=u(i)+aei(up,i,nmax)
      up(i-1)=up(i)+aei(upp,i,nmax)
         do it=1,2
         upp(i-1)=al*up(i-1)+cf(i-1)*u(i-1)
         up(i-1)=up(i)+aii(upp,i,nmax)
         u(i-1)=u(i)+aii(up,i,nmax)
         end do
      go to 120
  140 continue

c>>> FINE FASE B2 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

c>>> FASE B3: raccordo la funzione in r(mch) e ne calcolo la norma>>>>>>

c raccordo con continuita': riscalo u(i>mch)
      sc=uout/u(mch)
         do i=mch,nin
         up(i)=sc*up(i)
         u(i)=sc*u(i)
         end do
c calcolo norma
      ro=r(1)/sqrt(amesh)
      sn=ro**(2*l+3)/(2*l+3) + 2.0*c2*ro**(2*l+4)/(2*l+4)
     & +(2.0*c3+c2**2)*ro**(2*l+5)/(2*l+5)
         do i=1,nin
         sn=sn+al*r(i)*u(i)**2
         end do

c>>> FINE FASE B3 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

c>>> FASE B4: valutazione errore energia e decisione se andare avanti>>>

c teoria perturbativa per energia della prossima iterazione*********
c basata su relazione fra mismatch derivata e errore energia********
      deps=0.5*uout*(upout-up(mch))/(sn*al*r(mch))

c se l'errore relativo fra questa e la precedente energia e' minore di
c relerr (vedi input), normalizza u ed esci dalla subroutine (FASE C)
      if(abs(deps/eps) .lt. relerr) go to 170

c se invece l'errore relativo rispetto alla precedente energia e' 
c maggiore o uguale di relerr (vedi input), correggi l'energia e fai 
c un'altra iterazione (cioe' torna all'inizio della FASE B)
      if(abs(deps/eps) .gt. 0.25) deps=0.25*deps*abs(eps/deps)
      if(deps .gt. 0.0) emin=eps
      if(deps .lt. 0.0) emax=eps
      eps=eps+deps
      if(eps .gt. emax .or. eps .lt. emin) eps=0.5*(emax+emin)
      go to 10

c>>> FINE FASE B3 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c***********************************************************************
c*** FINE FASE B *******************************************************

c*** FASE C: normalizza funzione, pulisci array, metti nella variabile
c*** relerr l'effettivo errore relativo ed esci dalla subroutine********
c***********************************************************************

c normalizza la funzione d'onda
  170 cn=1.0/sqrt(sn)
         do i=1,nin
         up(i)=cn*up(i)
         u(i)=cn*u(i)
         end do

c pulizia: azzera la funzione d'onda al di la' di r(nin)
         do i=nin+1,mmax
         u(i)=0.0
         end do

c metti in relerr (output) l'effettivo errore relativo 
      relerr=abs(deps/eps)

c esci dalla subroutine
      return                                  !ritorna nel programma principale
c***********************************************************************
c*** FINE FASE c *******************************************************
      end                                       !fine modulo subroutine schroed

c 4 funzioni che realizzano le formule di adams (abramowitz & stegun p. 896)

      real function aio(y,j,nmax)                  !adams interpolation outward
      dimension y(nmax)
      aio=(4.1666667e-2)*(9.0*y(j+1)+19.0*y(j)
     & -5.0*y(j-1)+y(j-2))
      return
      end

      real function aei(y,j,nmax)                   !adams extrapolation inward
      dimension y(nmax)
      aei=-(4.1666667e-2)*(55.0*y(j)-59.0*y(j+1)
     & +37.0*y(j+2)-9.0*y(j+3))
      return
      end

      real function aii(y,j,nmax)                   !adams interpolation inward
      dimension y(nmax)
      aii=-(4.1666667e-2)*(9.0*y(j-1)+19.0*y(j)
     & -5.0*y(j+1)+y(j+2))
      return
      end

      real function aeo(y,j,nmax)                  !adams extrapolation outward
      dimension y(nmax)
      aeo=(4.1666667e-2)*(55.0*y(j)-59.0*y(j-1)
     & +37.0*y(j-2)-9.0*y(j-3))
      return
      end

c subroutine che stampa il messaggio contenuto in "messag" e ferma (stop)

      subroutine err(messag)
      character*40 messag
      write(70,1000) messag
 1000 format(/1x,' stop : ',a40/)
      stop
      end

c subroutine che costruisce la griglia

      subroutine griglia(r1,grigl,nx,nmax,r) !inizio subroutine: nome argomenti
      dimension r(nmax)                                     !dimensione array r
      r(1)=r1                                        !primo punto della griglia
         do i=2,nx                                  !ciclo: costruzione griglia
         r(i)=grigl*r(i-1)
         end do                                 !fine ciclo costruzione griglia
      return                                  !ritorna nel programma principale
      end                                       !fine modulo subroutine griglia
