program harmonic1 * * Risolve l'oscillatore armonico quantistico. * Integrazione in avanti e all'indietro con metodo di Numerov. * Ricerca degli autovalori mediante "shooting method". * Matching delle soluzioni a un punto di inversione. * implicit none * Numero massimo di punti nella griglia integer mshx parameter (mshx = 2000) real*8 pi parameter (pi=3.14159265358979323846d0) integer mesh, i, icl integer nodes, hnodes, ncross, kkk, iterate real*8 xmax, dx, ddx12, xmcl, norm, yicl, djump real*8 eup, elw, e real*8 x(0:mshx), y(0:mshx), p(0:mshx), vpot(0:mshx), f(0:mshx) character*80 fileout * * Input e inizializzazione della griglia * print '($,a)','Massimo x (valore tipico: 10) ? ' read*, xmax print '($,a,i5,a)', $ 'Numero punti griglia (massimo:',mshx,') ? ' read*,mesh if (mesh.gt.mshx) stop dx = xmax/mesh ddx12 = dx*dx/12.d0 * * Definizione del potenziale (simmetrico rispetto a x=0) * do i = 0, mesh x(i) = dfloat(i) * dx vpot(i) = 0.5d0 * x(i)*x(i) enddo print '($,a)','Nome file output = ' read '(a)',fileout open(1,file=fileout,form="formatted",status="new") * * Entrata loop di ricerca nuovo autovalore * 999 continue * * Input numero nodi (<0 per terminare) * print '($,a)', 'Numero nodi (-1=esci) ? ' read*,nodes if (nodes.lt.0) then close(1) stop endif * * Definizione limiti all'autovalore * eup = vpot(mesh) elw = eup do i=0,mesh elw = min(elw,vpot(i)) eup = max(eup,vpot(i)) enddo * * Definizione energia per singola valutazione * print '($,a)','Energia da provare (0=ricerca con bisezione) ? ' read*,e if ( e .eq. 0.d0 ) then * Modalita' ricerca automatica autovalore e = 0.5d0 * (elw + eup) iterate = 1 else * Modalita' test di un singolo autovalore iterate = 0 endif kkk = 0 * * Punto di entrata per calcolo soluzione a energia e * 1 continue kkk = kkk + 1 djump = 0.d0 * * Definizione funzione g usata nell'algoritmo di Numerov * f(0) = ddx12*(2.d0*(vpot(0)-e)) do i=1,mesh * f <= 0 regione permessa classicamente * f > 0 regione proibita classicamente f(i)=ddx12*2.d0*(vpot(i)-e) * Nessun f(i) puo' essere zero altrimenti * il cambio di segno non viene rilevato. * La linea seguente e' un trucco per evitare * di imbattersi in f(i) nulli: if( f(i) .eq. 0.d0 ) f(i)=1.d-20 * Registra l'indice a cui e' avvenuto * l'ultimo cambio di segno: if( f(i) .ne. sign(f(i),f(i-1)) ) icl=i enddo if(icl.ge.mesh-2) then print*,'Ultimo cambio segno troppo lontano.' stop endif * * f richiesta dal metodo di Numerov: * do i=0,mesh f(i) = 1.0d0 - f(i) enddo do i=0,mesh y(i) = 0.d0 enddo * * Determinazione funzione d'onda nei primi due punti * hnodes = nodes/2 if (2*hnodes.eq.nodes) then y(0) = 1.0d0 y(1) = 0.5d0*(12.d0-10.d0*f(0))*y(0)/f(1) else y(0) = 0.d0 y(1) = dx endif * * Integrazione verso l'esterno fermandosi all'ultimo punto * di inversione, e conteggio dei cambi di segno * ncross=0 do i=1,icl-1 y(i+1) = ((12.d0-10.d0*f(i))*y(i)-f(i-1)*y(i-1))/f(i+1) if ( y(i) .ne. sign(y(i),y(i+1)) ) ncross=ncross+1 enddo yicl = y(icl) * * Verifica il numero di cambi di segno * if ( iterate .ne. 0 ) then if (ncross.ne.hnodes) then print '(2i4,e16.8,f14.8)',kkk, ncross, djump, e * Numero errato di cambi di segno if (ncross.gt.hnodes) then * Avvenuti piu' cambi di segno del previsto. * L'energia e' troppo alta e va abbassata * (--> semiintervallo inferiore). eup = e else * Avvenuti meno cambi di segno del previsto. * L'energia e' troppo bassa e va alzata * (--> semiintervallo superiore). elw = e endif * Nuovo valore da provare: e = 0.5d0 * (eup+elw) * Riiteriamo comunque (il numero di nodi era errato) go to 1 endif endif * * Il numero di nodi e' giusto. Completiamo l'integrazione * e utilizziamo il matching della derivata al punto di * giunzione come criterio per la ricerca dell'autovalore. * y(mesh) = dx y(mesh-1) = (12.d0-10.d0*f(mesh))*y(mesh)/f(mesh-1) * * Integrazione dall'esterno verso l'interno * do i=mesh-1,icl+1,-1 y(i-1)=((12.d0-10.d0*f(i))*y(i)-f(i+1)*y(i+1))/f(i-1) enddo * * Riscala per ottenere matching a icl: * yicl = yicl/y(icl) do i=icl,mesh y(i) = y(i)*yicl enddo * * Rinormalizzazione della funzione d'onda in modo che * int |psi|^2 dx = 1: * norm = 0.d0 do i=1,mesh norm = norm + y(i)*y(i) enddo norm = dx * (2.d0 * norm + y(0)*y(0)) norm = sqrt(norm) do i=0,mesh y(i) = y(i) / norm enddo * * Calcolo discontinuita' della derivata prima al punto di * raccordo * if ( iterate .ne. 0 ) then i = icl * Calcola y'(i;RIGHT) - y'(i;LEFT) djump = ( y(i+1) + y(i-1) - (14.d0-12.d0*f(i))*y(i) ) / dx print '(2i4,e16.8,f14.8)',kkk, ncross, djump, e if (djump*y(i).gt.0.d0) then * Energia troppo alta --> semiintervallo inferiore eup = e else * Energia troppo bassa --> semiintervallo superiore elw = e endif e = 0.5d0 * (eup+elw) if (eup-elw .gt. 1.d-10) go to 1 endif print '(2i4,e16.8,f14.8)',kkk, ncross, djump, e * * La convergenza e' stata raggiunta, o non erano richieste * iterazioni. * * * Calcolo della densita' classica relativa all'energia e * (per confronto): * xmcl = sqrt(2.d0*e) do i=icl,mesh p(i) = 0.d0 enddo do i=0,icl-1 p(i) = 1.d0/sqrt(xmcl**2 - x(i)**2)/pi enddo * Semiasse negativo: do i=mesh,1,-1 write (1,'(f7.3,3d16.8,f12.6)') $ -x(i),(-1)**nodes*y(i),y(i)*y(i),p(i),vpot(i) enddo * Semiasse positivo: do i=0,mesh write (1,'(f7.3,3d16.8,f12.6)') % x(i),y(i),y(i)*y(i),p(i),vpot(i) enddo go to 999 end