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, 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 hartold=hartree 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 etovar=2*eps-hartold write(80,1000) i, etotal, etovar !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