PROGRAM HELISIMP C C*********************************************************************** C C * Càlcul variacional senzill dels estats excitats 2P i 3D dels C àtoms (ions) de dos electrons C C * El programa minimitza l'energia E(x,y) calculada a classe C respecte els paràmetres x,y amb un mètode de tipus 'SIMPLEX' C (x=alpha, y=beta en la notació de classe) C C M. Centelles C mario@ecm.ub.es C C*********************************************************************** C IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(NP=2) COMMON/VAR/ Z,ILEVEL,ISIM COMMON/ENERGS/ EI1,EI2,EJ,EK DIMENSION PARAM(NP),D1(NP,NP+1),F1(NP+1),D0(NP) EXTERNAL FUNC c HART= 27.2113845D0 c c Factor de precisió i número màxim d'iteracions per a la c subrutina de minimització TOL= 1.D-16 NITER= 1000 c c Lectura de dades c WRITE(6,*) '-----------------------------------------------------' WRITE(6,*) ' ATOMS (IONS) DE DOS ELECTRONS' WRITE(6,*) ' Estats excitats 2P i 3D: càlcul variacional senzill' WRITE(6,*) '-----------------------------------------------------' 10000 WRITE(6,*) 'Càrrega nuclear Z ...' READ(5,*) Z WRITE(6,*) 'Estat 2P (1s,2p) o 3D (1s,3d) [entra 0 o 1] ...' READ(5,*) ILEVEL WRITE(6,*) 'Para o orto [entra 0 o 1] ...' READ(5,*) ISIM c c Valors inicials per començar a buscar el mínim c PARAM(1)= 1.D0 PARAM(2)= (Z-1.D0)/Z CALL FUNC(PARAM,ENER) FSTART= ENER c c Minimització c CALL DSIMPL(TOL,0,NITER,0,NP,NP+1,PARAM,D1,F1,D0,FUNC,FSTART) c c Valors finals i presentació c CALL FUNC(PARAM,ENER) x= PARAM(1) y= PARAM(2) WRITE(6,*) '-----------------------------------------------------' WRITE(6,*) 'Resultat del càlcul variacional:' WRITE(6,*) WRITE(6,*) 'ENERGIA ...' WRITE(6,*) ' E[1s,nl]=', ENER,' u.a. =', ENER*HART,' eV' WRITE(6,*) WRITE(6,*) 'Càrregues efectives ...' WRITE(6,*) ' Z(1s)= x*Z=',x*Z, ', Z(nl)= y*Z=',y*Z WRITE(6,*) WRITE(6,*)'Integrals I[1s], I[nl], J[1s;nl], K[1s;nl] (u.a.) ...' WRITE(6,*) EI1, EI2, EJ, EK WRITE(6,*) WRITE(6,*)'Integrals I[1s], I[nl], J[1s;nl], K[1s;nl] (eV) ...' WRITE(6,*) EI1*HART, EI2*HART, EJ*HART, EK*HART WRITE(6,*) '-----------------------------------------------------' WRITE(6,*) c c Per comparar, també estimarem el valor de E[1s,nl] en una c aproximació més esquemàtica d'apantallament total: c càrrega efectiva Z per a l'electró intern i c càrrega efectiva Z-1 per a l'electró extern c IF (ILEVEL.EQ.0) THEN c c cas 2P ENERX= -0.5d0*( Z**2 + (Z-1)**2/2.d0**2) ELSE c c cas 3D ENERX= -0.5d0*( Z**2 + (Z-1)**2/3.d0**2) ENDIF c WRITE(6,*) 'Resultat d''aproximació d''apantallament total:' WRITE(6,*) '(no distingeix entre nivells orto i para)' WRITE(6,*) WRITE(6,*) ' E[1s,nl]= - 0.5*Z^2/1^2 - 0.5*(Z-1)^2/n^2' WRITE(6,*) ' =', ENERX,' u.a. =', ENERX*HART,' eV' WRITE(6,*) '-----------------------------------------------------' WRITE(6,*) WRITE(6,*) 'Entra 0 per a acabar, 1 per a un altre càlcul ...' READ(5,*) MORE WRITE(6,*) GOTO (10000,10002),MORE 10002 STOP END SUBROUTINE FUNC(PARAM,ENER) c c Càlcul de l'energia ENER(x,y;Z) c IMPLICIT REAL*8 (A-H,O-Z) COMMON/VAR/ Z,ILEVEL,ISIM COMMON/ENERGS/ EI1,EI2,EJ,EK DIMENSION PARAM(*) C x= PARAM(1) y= PARAM(2) IF (ILEVEL.EQ.0) THEN c c Integrals I,J,K pel cas 2P EI1= (x**2/2.d0 - x)*Z**2 EI2= (y**2/2.d0 - y)*Z**2/4.d0 EJ= y*Z/4.d0 - (x*y**5*Z)/(2.d0*x+y)**5 - & (y**5*Z)/(4.d0*(2.d0*x+y)**4) EK= (112.d0*x**3*y**5*Z)/(3.d0*(2.d0*x+y)**7) ELSE c c Integrals I,J,K pel cas 3D EI1= (x**2/2.d0 - x)*Z**2 EI2= (y**2/2.d0 - y)*Z**2/9.d0 EJ= y*Z/9.d0 - (x*y**7*Z)/(3.d0*x+y)**7 - & (y**7*Z)/(9.d0*(3.d0*x+y)**6) EK= (324.d0*x**3*y**7*Z)/(5.d0*(3.d0*x+y)**9) ENDIF c c Signe de K segons sigui para o orto EK= float(1-2*ISIM)*EK c c Energia ENER= EI1 + EI2 + EJ + EK RETURN END SUBROUTINE DSIMPL(TOL,ND,NI,NS,N3,N4,D2,D1,F1,D0,FUNC,FSTART) c c Subrutina de minimització d'una funció FTOTAL que depèn de c D2(N3) variables respecte d'aquestes. c El càlcul de FTOTAL per a cada conjunt de valors de D2(N3) es c fa cridant la subrutina externa FUNC(D2,FTOTAL). c DSIMPL assaja diferents valors de D2(N3) i evoluciona de forma c que el corresponent valor de FTOTAL va fent-se menor. c c Per a detalls del mètode SIMPLEX de minimització/maximització: c http://www.nrbook.com/b/bookfpdf/f10-4.pdf c (versió online del llibre "Numerical Recipes, The Art of c Scientific Computing") c c c * Si vols veure l'evolució de la minimització descomenta c els WRITE(6,101), WRITE(6,103) i WRITE(6,109) de baix. c c IMPLICIT REAL*8 (A-H,O-Z) DIMENSION D1(N3,N4),F1(N4),D0(N3),D2(N3) 101 FORMAT(' INIT',1P,E20.7) 103 FORMAT( ' ',1P,7E16.6) C INITIAL LOOP C RN0=59.D0 SF=1.*N4 NM=(N4-1)**2 NC=0 NIR=5 DO 20 J=1,N4 IF(J.GT.1.AND.ND.EQ.0) GO TO 8 C READ(NIR,103)(D2(I),I=1,N3),F1(J) F1(1)=FSTART ccc WRITE(6,103)(D2(I),I=1,N3) 8 CONTINUE DO 10 I=1,N3 10 D1(I,J)=D2(I) IF(ND.NE.2) CALL FUNC(D2,FTOTAL) IF(ND.NE.2) F1(J)=FTOTAL ccc WRITE(6,101) F1(J) ccc WRITE(6,103) D2 IF(ND.NE.0.OR.J.EQ.N4) GO TO 20 DO 15 I=1,N3 15 D2(I)=D1(I,1) IF((J+NS).GT.N3) GO TO 20 RN0=DMOD(1.2207031250D9*RN0,3.43597383680D10) RND=SNGL(RN0*2.910383045673370D-11) D2(J+NS)=(1.+.1*(2.*RND-1.))*D1(J+NS,1) IF(D2(J+NS).EQ.0.)D2(J+NS)=.1 20 CONTINUE 40 CONTINUE C SORT LOOP I9=1 I1=I9 I2=2 S=F1(1) B1=S B2=F1(2) NC=NC+1 DO 70 I=2,N4 IF(S.LT.F1(I)) GO TO 50 S=F1(I) I9=I 50 IF(B1.GT.F1(I)) GO TO 60 B2=B1 I2=I1 B1=F1(I) I1=I GO TO 70 60 IF(B2.GT.F1(I)) GO TO 70 B2=F1(I) I2=I 70 CONTINUE IF(DABS((B1-S)/S).LT.TOL ) GO TO 800 IF(NC.EQ.NI) GO TO 800 IF(DABS((B1-S)/S).LT.5.*TOL) GO TO 300 C FORMING D0 DO 100 J=1,N3 D0(J)=0. DO 90 I=1,N4 IF(I.EQ.I1) GO TO 90 D0(J)=D0(J)+D1(J,I) 90 CONTINUE D0(J)=D0(J)/(N4-1) 100 CONTINUE C REFLECTION DO 110 I=1,N3 D2(I)=2.*D0(I)-D1(I,I1) 110 CONTINUE CALL FUNC(D2,FTOTAL) FT=FTOTAL IF(FT.GE.B2) GO TO 150 120 F1(I1)=FT DO 130 I=1,N3 130 D1(I,I1)=D2(I) c ccc IF (FT.LT.F1(I9)) WRITE(6,109) FT,(D2(I),I=1,N3) c 109 FORMAT(' NEW MIN',1P,E17.8,' CONST ARE', (1P,8E15.6)) IF(FT.LT.S) GO TO 200 GO TO 40 C CONTRACTION 150 FS=FT IF(FT.GT.B1) GO TO 165 FS=B1 B1=FT DO 160 I=1,N3 160 D1(I,I1)=D2(I) 165 CONTINUE AM=.75 168 CONTINUE DO 170 I=1,N3 D2(I)=AM*D1(I,I1)+(1.-AM)*D0(I) 170 CONTINUE CALL FUNC(D2,FTOTAL) FT=FTOTAL S=-1.E35 IF(FT.LT.B2)GO TO 120 IF(AM.LT..3) GO TO 300 AM=AM-.25 GO TO 168 C EXPANSION 200 S=FT DO 210 I=1,N3 D2(I)=2*D2(I)-D0(I) 210 CONTINUE CALL FUNC(D2,FTOTAL) FT=FTOTAL IF(FT.GT.S) GO TO 40 S=-1.E35 GO TO 120 C SHRINKAGE 300 CONTINUE C FORMING THE AVERAGE OF THE SQUARES OF THE DISTANCES FROM THE MIN POINT DO 280 J=1,N3 D0(J)=0. DO 270 I=1,N4 IF(I.EQ.I9)GO TO 270 D0(J)=D0(J)+(D1(J,I)-D1(J,I9))**2 270 CONTINUE 280 D0(J)=DSQRT(D0(J)/(N4-1)) SF=0.75D0*SF ccc WRITE(6,103) SF C CALL RSEED(75579321+1000*NC,179321195+15*NC) MM=NS J=0 310 J=J+1 IF(J.EQ.I9) GO TO 310 MM=MM+1 IF(MM.GT.N3)GO TO 40 DO 340 I=1,N3 D2(I)=D1(I,I9) 340 CONTINUE RN0=DMOD(1.2207031250D9*RN0,3.43597383680D10) RND=SNGL(RN0*2.910383045673370D-11) D2(MM)=SF*(2.*RND-1.)*(D0(MM))+D1(MM,I9) CALL FUNC(D2,FTOTAL) F1(J)=FTOTAL ccc WRITE(6,103) D2 DO 345 I=1,N3 345 D1(I,J)=D2(I) I1=J FT=F1(J) IF(FT.LT.F1(I9)) GO TO 120 GO TO 310 800 CONTINUE DO 810 J=1,N4 ccc WRITE(6,103) (D1(I,J),I=1,N3),F1(J) 810 CONTINUE DO 820 J=1,N3 820 D2(J)=D1(J,I9) RETURN END