CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PROGRAM HF3 C Hartree-Fock solutions of small atomic systems in the C filling approximation CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C output units INTEGER MUNIT !output unit INTEGER OUNIT !density output unit MUNIT=12 OUNIT=13 OPEN (MUNIT,FILE='P3.OUT',STATUS='UNKNOWN') OPEN (OUNIT,FILE='P32.OUT',STATUS='UNKNOWN') C CALL INIT !display header screen, setup parameters 5 CONTINUE !main loop/ execute once for each set of param CALL ARCHON !find the Hartree-Fock wavefunctions GOTO 5 END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ARCHON C find the Hartree-Fock wavefunctions for the specified atom CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Local variables: REAL E(MAXSTT+1,8) !all energies of all states REAL FOCK(0:MAXSTP,MAXSTT) !Fock terms REAL RHO(0:MAXSTP) !density REAL PSTOR(0:MAXSTP,MAXSTT) !radial wavefunction REAL PHI(0:MAXSTP) !electron potential REAL ESP !single particle energy of state INTEGER ITER !iteration index INTEGER STATE !single particle state index REAL ZSTAR !optimal effective nuclear charge INTEGER ISTOP,ISTART !current limits on iteration INTEGER MUNIT !output unit MUNIT=12 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C begin iterations with a good guess MIX=1. !no old density to mix with new ZSTAR=Z CALL HYDRGN(ZSTAR,PSTOR) !find hydrogenic wavefunctions CALL ENERGY(E,FOCK,RHO,PHI,PSTOR) !and energy C optimal ZSTAR using Virial theorem ZSTAR=-Z*(E(NSTATE+1,IVTOT)/(2*E(NSTATE+1,IKTOT)) ) CALL HYDRGN(ZSTAR,PSTOR) !find new hydrogenic wavefunctions CALL ENERGY(E,FOCK,RHO,PHI,PSTOR) ! and energies C C output summary of parameters CALL PRMOUT(MUNIT,ZSTAR) C MIX=.5 !mix old and new density for stability ITER=0 !zero iteration counters ISTART=0 ISTOP=0 C C output initial energies CALL TXTOUT(MUNIT,ITER,E,RHO,PSTOR) C 10 CONTINUE !loop over iterations ISTART=ISTOP+1 !update iteration counters ISTOP=ISTOP+NITER DO 100 ITER=ISTART,ISTOP DO 50 STATE=1,NSTATE IF (NOCC(STATE) .NE. 0) THEN !loop over all occpd states ESP=E(STATE,ITOT) !single particle energy IF (ESP .GT. 0) ESP=-10 !keep particle bound C find single part wavefunc CALL SNGWFN(STATE,ESP,PSTOR,FOCK,PHI) END IF 50 CONTINUE C calculate new energies and output CALL ENERGY(E,FOCK,RHO,PHI,PSTOR) CALL TXTOUT(MUNIT,ITER,E,RHO,PSTOR) 100 CONTINUE C C allow for more iterations PRINT*,' Increase number of iterations?' PRINT*,' Then add [ ] more steps [enter 0 to stop]:' READ(*,*) NITER IF (NITER .NE. 0 ) THEN GOTO 10 ELSE PRINT*,' Output in P3.out: energies of each orbital' PRINT*,' Output in P32.out: densities' STOP ENDIF C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE HYDRGN(ZSTAR,PSTOR) C creates hydrogenic wavefunctions PSTOR with nuclear charge=ZSTAR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Input variables: REAL ZSTAR !effective nuclear charge REAL PSTOR(0:MAXSTP,MAXSTT) !radial wavefunction C Local variables: INTEGER IR !radial index REAL RSTAR !scaled radius REAL ERSTAR !useful exponential INTEGER STATE !state index REAL NORM !norm of wavefunction CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C store radial parts of hydrogenic wavefunctions DO 20 IR=0,NR RSTAR=IR*DR*ZSTAR/ABOHR !scaled radius ERSTAR=EXP(-RSTAR/2.) !useful exponential IF (NOCC(1) .NE. 0) PSTOR(IR,1)=RSTAR*ERSTAR**2 IF (NOCC(2) .NE. 0) PSTOR(IR,2)=(2-RSTAR)*RSTAR*ERSTAR IF (NOCC(3) .NE. 0) PSTOR(IR,3)=(RSTAR**2)*ERSTAR 20 CONTINUE C C normalize wavefunctions DO 40 STATE=1,NSTATE IF (NOCC(STATE) .NE. 0) THEN NORM=0. DO 30 IR=0,NR NORM=NORM+PSTOR(IR,STATE)**2 30 CONTINUE NORM=1./(SQRT(NORM*DR)) DO 35 IR=0,NR PSTOR(IR,STATE)=PSTOR(IR,STATE)*NORM 35 CONTINUE END IF 40 CONTINUE C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE ENERGY(E,FOCK,RHO,PHI,PSTOR) C subroutine to calculate the energies of a normalized C set of single-particle wavefunctions (PSTOR); C also calculates Fock terms, density, and electron potential CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Passed variables: REAL PSTOR(0:MAXSTP,MAXSTT) !radial wavefunction (input) REAL FOCK(0:MAXSTP,MAXSTT) !Fock terms (output) REAL PHI(0:MAXSTP) !electron potential (output) REAL RHO(0:MAXSTP) !density (output) REAL E(MAXSTT+1,8) !all energies of all states (output) C Local variables: INTEGER STATE !state index INTEGER IR !radial index INTEGER IE !energy index REAL R !current radius REAL LL1 !square of angular momentum REAL PM,PZ,PZ2 !values to calc d(rho)/dr CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL SOURCE(PSTOR,FOCK,RHO) !calculate Fock terms and density CALL POISSN(PHI,RHO) !calc potntl due to electron charge C DO 20 IE=IKEN,ITOT !zero total energies E(NSTATE+1,IE)=0 20 CONTINUE C DO 38 STATE=1,NSTATE IF (NOCC(STATE) .NE. 0) THEN C DO 10 IE=IKEN,ITOT !zero energy for this state E(STATE,IE)=0 10 CONTINUE C LL1=ANGMOM(STATE)*(ANGMOM(STATE)+1) PM=0 DO 48 IR=1,NR !integrate the energy densities R=IR*DR PZ=PSTOR(IR,STATE) PZ2=PZ**2 C E(STATE,IKEN)=E(STATE,IKEN)+(PZ-PM)**2 !kinetic E(STATE,ICEN)=E(STATE,ICEN)+PZ2*LL1/R**2 !centrifugal E(STATE,IVEN)=E(STATE,IVEN)-PZ2/R !elec-nucl E(STATE,IVEE)=E(STATE,IVEE)+PHI(IR)*PZ2 !elec-elec E(STATE,IVEX)=E(STATE,IVEX)+FOCK(IR,STATE)*PZ !exchange C PM=PZ !roll values for derivative 48 CONTINUE C C put in constant factors E(STATE,IKEN)=E(STATE,IKEN)*HBARM/(2*DR) E(STATE,ICEN)=E(STATE,ICEN)*DR*HBARM/2 E(STATE,IVEN)=E(STATE,IVEN)*ZCHARG*DR E(STATE,IVEE)=E(STATE,IVEE)*DR E(STATE,IVEX)=E(STATE,IVEX)*DR C C calculate totals for this state E(STATE,IKTOT)=E(STATE,IKEN)+E(STATE,ICEN) E(STATE,IVTOT)=E(STATE,IVEN)+E(STATE,IVEE)+E(STATE,IVEX) E(STATE,ITOT)=E(STATE,IVTOT)+E(STATE,IKTOT) C C add this state's contribution to total energy DO 30 IE=IKEN,IVEX E(NSTATE+1,IE)=E(NSTATE+1,IE)+E(STATE,IE)*NOCC(STATE) 30 CONTINUE C END IF 38 CONTINUE C STATE=NSTATE+1 !calculate total energies C don't double count electron-electron and exchange energies E(STATE,IVEE)=E(STATE,IVEE)/2 E(STATE,IVEX)=E(STATE,IVEX)/2 C C find total kinetic, potential and total energies E(STATE,IKTOT)=E(STATE,IKEN)+E(STATE,ICEN) E(STATE,IVTOT)=E(STATE,IVEN)+E(STATE,IVEE)+E(STATE,IVEX) E(STATE,ITOT)=E(STATE,IVTOT)+E(STATE,IKTOT) C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SOURCE(PSTOR,FOCK,RHO) C subroutine to compute the density and the Fock terms given PSTOR CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Passed variables: REAL PSTOR(0:MAXSTP,MAXSTT)!radial wavefunction (input) REAL RHO(0:MAXSTP) !density (output) REAL FOCK(0:MAXSTP,MAXSTT) !Fock terms (output) C Local variables: REAL DF !increment in Fock term REAL FAC !constant factor in Fock term INTEGER IR !radial index REAL R !current radius REAL RLAM !r**lambda REAL RLAM1 !r**(lambda+1) INTEGER STATE !indexes state INTEGER STATE2 !indexes second state in Fock term REAL SUM !sum for Fock integrals REAL TERM !temporary term for sum INTEGER L1,L2 !angular momentum for two states INTEGER LSTART,LSTOP !limits on sum of L1+L2 INTEGER LAM !current value of ang mom sum REAL THREEJ !three j coefficient CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO 20 IR=1,NR !include fraction of old density RHO(IR)=(1.-MIX)*RHO(IR) 20 CONTINUE C DO 30 STATE=1,NSTATE IF (NOCC(STATE) .NE. 0) THEN !loop over occupied states C DO 40 IR=1,NR !contribution of this state to density RHO(IR)=RHO(IR)+MIX*NOCC(STATE)*(PSTOR(IR,STATE)**2) FOCK(IR,STATE)=0 !zero FOCK term 40 CONTINUE C C begin calculation of Fock term DO 50 STATE2=1,NSTATE IF (NOCC(STATE2) .NE. 0) THEN !loop over occupied states C L1=ANGMOM(STATE) L2=ANGMOM(STATE2) LSTART=IABS(L1-L2) !limits on ang mom sum LSTOP=L1+L2 C DO 60 LAM=LSTART,LSTOP,2 !loop over ang mom values C CALL SQR3J(L1,L2,LAM,THREEJ) FAC=-CHARGE/2*NOCC(STATE2)*THREEJ C SUM=0 DO 80 IR=1,NR !outward integral for Fock R=IR*DR RLAM=R**LAM TERM=PSTOR(IR,STATE2)*PSTOR(IR,STATE)*RLAM/2 SUM=SUM+TERM DF=PSTOR(IR,STATE2)*FAC*SUM*DR/(RLAM*R) FOCK(IR,STATE)=FOCK(IR,STATE)+DF SUM=SUM+TERM 80 CONTINUE C SUM=0 DO 90 IR=NR,1,-1 !inward integral for Fock R=IR*DR RLAM1=R**(LAM+1) TERM=PSTOR(IR,STATE2)*PSTOR(IR,STATE)/RLAM1/2 SUM=SUM+TERM DF=PSTOR(IR,STATE2)*FAC*SUM*DR*RLAM1/R FOCK(IR,STATE)=FOCK(IR,STATE)+DF SUM=SUM+TERM 90 CONTINUE C 60 CONTINUE !end loop over lam END IF 50 CONTINUE !end loop over state2 END IF 30 CONTINUE !end loop over state C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE POISSN(PHI,RHO) C subroutine to solve Poisson's equation for the direct potential C given the electron density RHO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Passed variables: REAL PHI(0:MAXSTP) !electron potential (output) REAL RHO(0:MAXSTP) !density (input) C Local variables: REAL CON !useful constant INTEGER IR !radial index REAL M !linear behavior to subtract REAL R !radius REAL SM,SP,SZ !source terms REAL SUM !total charge CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUM = 0. !quadrature of density to get DO 19 IR = 1,NR ! initial value for PHI(1) SUM = SUM+RHO(IR)/REAL(IR) 19 CONTINUE C CON = DR**2/12 !initial values for outward integ SM = 0 SZ = -CHARGE*RHO(1)/DR PHI(0) = 0 PHI(1) = SUM*CHARGE*DR C DO 29 IR = 1,NR-1 !Numerov algorithm SP=-CHARGE*RHO(IR+1)/((IR+1)*DR) PHI(IR+1)=2*PHI(IR)-PHI(IR-1)+CON*(10*SZ+SP+SM) SM=SZ SZ=SP 29 CONTINUE C M=(PHI(NR)-PHI(NR-10))/(10*DR) !subtract off linear behavior DO 39 IR=1,NR R=IR*DR PHI(IR)=PHI(IR)/R-M !factor of 1/r for true potl 39 CONTINUE C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SNGWFN(STATE,E,PSTOR,FOCK,PHI) C subroutine to solve the single-particle wavefunction as an C inhomogeneous boundary-value problem for a given state, energy E, C source term (FOCK), and potential (PHI) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Passed variables: INTEGER STATE !which single particle state (input) REAL E !single particle energy (input) REAL FOCK(0:MAXSTP,MAXSTT) !Fock terms (input) REAL PHI(0:MAXSTP) !electron potential (input) REAL PSTOR(0:MAXSTP,MAXSTT) !radial wavefunction (output) C Local variables: REAL PSIIN(0:MAXSTP),PSIOUT(0:MAXSTP) !homogeneous solutions INTEGER NR2 !midpoint for the lattice REAL LL1 !angular momentum REAL DRHBM !useful constant REAL K2M,K2Z,K2P !local wave numbers REAL NORM !normalization REAL R !current radius REAL SUM !sum for integration REAL TERM !temp term in integration REAL WRON !Wronskian at middle of lattice INTEGER IR !radial index CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DRHBM=DR**2/HBARM/6 !useful constants LL1=ANGMOM(STATE)*(ANGMOM(STATE)+1)*HBARM/2 C K2M=0 !integrate outward homogeneous soln K2Z=DRHBM*(E-PHI(1)+(ZCHARG-LL1/DR)/DR) PSIOUT(0)=0 PSIOUT(1)=1.0E-10 DO 49 IR=2,NR R=DR*IR !Numerov algorithm K2P=DRHBM*(E-PHI(IR)+(ZCHARG-LL1/R)/R) PSIOUT(IR)=(PSIOUT(IR-1)*(2-10*K2Z)- + PSIOUT(IR-2)*(1+K2M))/(1+K2P) K2M=K2Z !roll values K2Z=K2P 49 CONTINUE C K2P=0 !integrate inward homogeneous soln R=(NR-1)*DR K2Z=DRHBM*(E-PHI(NR-1)+(ZCHARG-LL1/R)/R) PSIIN(NR)=0 PSIIN(NR-1)=1.0E-10 DO 59 IR=NR-2,1,-1 R=DR*IR !Numerov algorithm K2M=DRHBM*(E-PHI(IR)+(ZCHARG-LL1/R)/R) PSIIN(IR)=(PSIIN(IR+1)*(2-10*K2Z)- + PSIIN(IR+2)*(1+K2P))/(1+K2M) K2P=K2Z !roll values K2Z=K2M 59 CONTINUE C NR2=NR/2 !Wronskian at middle of mesh WRON=(PSIIN(NR2+1)-PSIIN(NR2-1))/(2*DR)*PSIOUT(NR2) WRON=WRON-(PSIOUT(NR2+1)-PSIOUT(NR2-1))/(2*DR)*PSIIN(NR2) C SUM=0 !outward integral in Green's soln DO 69 IR=1,NR TERM=-PSIOUT(IR)*FOCK(IR,STATE)/2 SUM=SUM+TERM PSTOR(IR,STATE)=PSIIN(IR)*SUM*DR SUM=SUM+TERM 69 CONTINUE C SUM=0 !inward integral in Green's soln DO 79 IR=NR,1,-1 TERM=-PSIIN(IR)*FOCK(IR,STATE)/2 SUM=SUM+TERM PSTOR(IR,STATE)=(PSTOR(IR,STATE)+PSIOUT(IR)*SUM*DR)/WRON SUM=SUM+TERM 79 CONTINUE C IF (STATE .EQ. 2) THEN !keep 1s and 2s states orthogonal SUM=0 DO 89 IR=1,NR SUM=SUM+PSTOR(IR,1)*PSTOR(IR,2) 89 CONTINUE SUM=SUM*DR DO 99 IR=1,NR PSTOR(IR,2)=PSTOR(IR,2)-SUM*PSTOR(IR,1) 99 CONTINUE END IF C NORM=0 !normalize the soln DO 18 IR=1,NR NORM=NORM+PSTOR(IR,STATE)**2 18 CONTINUE NORM=1/SQRT(NORM*DR) DO 28 IR=1,NR PSTOR(IR,STATE)=PSTOR(IR,STATE)*NORM 28 CONTINUE C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE SQR3J(L1,L2,LAM,THREEJ) C subroutine to calculate square of the 3-j coefficient appearing C in the exchange energy CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Passed variables: INTEGER L1,L2,LAM !angular momentum (input) REAL THREEJ !three j coefficient (output) C Local variables: REAL DELTA !intermediate term for calc of 3J INTEGER IMAX,P !useful combinations of ang mom CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMAX=L1+L2+LAM+1 !useful combinations P=(L1+L2+LAM)/2 C DELTA=FACTRL(L1+L2-LAM)*FACTRL(-L1+L2+LAM) !temp values DELTA=DELTA*FACTRL(L1-L2+LAM)/FACTRL(IMAX) THREEJ=DELTA*(FACTRL(P)**2) !calculate 3J squared THREEJ=THREEJ/(FACTRL(P-L1)**2) THREEJ=THREEJ/(FACTRL(P-L2)**2) THREEJ=THREEJ/(FACTRL(P-LAM)**2) C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE INIT C intializes constants, displays header screen, C gets parameters from screen C ends program on request C calculates all derivative parameters CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Local parameters: INTEGER I !index for factorial loop INTEGER MUNIT !output unit INTEGER IOPT !run option CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MUNIT=12 C display header screen WRITE(MUNIT,*)'PROJECT 3' WRITE(MUNIT,*)'Hartree-Fock solutions of small atomic systems' WRITE(MUNIT,*)'in the filling approximation' C C calculate constants c atomic constants HBARM=7.6359 !hbar*2/(mass) CHARGE=14.409 !charge of the electron ABOHR=HBARM/CHARGE !Bohr radius C NSTATE=MAXSTT !descriptions of states ANGMOM(1)=0 ANGMOM(2)=0 ANGMOM(3)=1 ID(1)=' 1S ' ID(2)=' 2S ' ID(3)=' 2P ' ID(4)='TOTAL' C FACTRL(0)=1 !factorials DO 10 I=1,10 FACTRL(I)=I*FACTRL(I-1) 10 CONTINUE C C get input from terminal C C default parameters for a test run Z=6. NOCC(1)=2 NOCC(2)=2 NOCC(3)=2 DR=0.1 RMAX=2. NITER=4 PRINT*,' enter 0 [1] for default [your] parameters' READ(*,*)IOPT IF(IOPT .EQ. 0)GOTO 111 C PRINT*,'Enter nuclear charge - Min = 1, Max = 10' READ(*,*) Z C PRINT*,'Enter number of electrons in the 1s state (0, 1, or 2)' READ(*,*) NOCC(1) C PRINT*,'Enter number of electrons in the 2s state (0, 1, or 2)' READ(*,*) NOCC(2) C PRINT*,'Enter number of electrons in the 2p state (0 to 6)' READ(*,*) NOCC(3) C 123 PRINT*,'Enter radial step size (Angstroms) (.0001 to .5)' READ(*,*) DR C PRINT*,'Enter outer radius of the lattice (Angstroms)' PRINT*,' (.01 to 10.)' READ(*,*) RMAX C PRINT*,'Enter number of iterations (1 to 50)' READ(*,*) NITER C 111 ZCHARG=Z*CHARGE !nuclear charge NR=INT(RMAX/DR) !number of radial steps C ensure that the number of radial steps is not greater than the C size of the arrays IF ( NR .GT. MAXSTP) THEN WRITE(*,15) REAL(NR),MAXSTP PRINT*,' Enter a larger step' GOTO 123 END IF C 15 FORMAT (' Total number of radial steps (=',1PE9.2, + ') is larger than maxstp (=',i5,')') C NOCC(NSTATE+1)=0 !total number of electrons DO 20 I=1,NSTATE NOCC(NSTATE+1)=NOCC(NSTATE+1)+NOCC(I) 20 CONTINUE C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE PRMOUT(MUNIT,ZSTAR) C outputs text results to the specified unit CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables INCLUDE 'PAR.P3' C Passed variables INTEGER MUNIT !current unit number (input) REAL ZSTAR !optimal charge (input) INTEGER IS !indexes states (input) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C WRITE (MUNIT,19) WRITE (MUNIT,21) WRITE (MUNIT,27) Z,ZSTAR WRITE (MUNIT,23) RMAX,DR WRITE (MUNIT,30) (NOCC(IS),IS=1,NSTATE) WRITE (MUNIT,33) WRITE (MUNIT,19) C 19 FORMAT (' ') 21 FORMAT (' Output from project 3:', + ' Hartree-Fock solutions for small atomic systems') 27 FORMAT (' nuclear charge=',F6.3,5X,' zstar =',F6.3) 23 FORMAT (' Rmax (Angstroms)=',F6.3,5X, + ' radial step (Angstroms)=', + 1PE12.5) 30 FORMAT (' Occupation of the states are:',3(2X,I2)) 33 FORMAT (' All energies are in eV') C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE TXTOUT(MUNIT,ITER,E,RHO,PSTOR) C writes energies for one iteration to the requested unit CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C Global variables: INCLUDE 'PAR.P3' C Passed variables: REAL PSTOR(0:MAXSTP,MAXSTT) !radial wavefunction REAL RHO(0:MAXSTP) !density REAL E(MAXSTT+1,8) !all energies of all states (input) INTEGER MUNIT !current unit (input) INTEGER ITER !current iteration (input) C Local variables: REAL R !radius for OUNIT REAL DEN !density for gunit INTEGER IS !state index INTEGER IR !indexes radius INTEGER STATE !current state DIMENSION DEN(MAXSTT+1) INTEGER OUNIT !density output unit CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C OUNIT=13 WRITE (MUNIT,20) ITER WRITE (MUNIT,30) WRITE (*,20) ITER WRITE (*,30) C DO 50 IS=1,NSTATE+1 IF (NOCC(IS) .NE. 0) THEN !occupied states only WRITE (MUNIT,35) ID(IS),NOCC(IS),E(IS,IKTOT),E(IS,IVEN), + E(IS,IVEE),E(IS,IVEX),E(IS,IVTOT),E(IS,ITOT) WRITE (*,35) ID(IS),NOCC(IS),E(IS,IKTOT),E(IS,IVEN), + E(IS,IVEE),E(IS,IVEX),E(IS,IVTOT),E(IS,ITOT) END IF 50 CONTINUE C 20 FORMAT (26X,'------- Iteration ',I2,' -------') 30 FORMAT (2X,'State',2X,'Nocc',4X,'Ktot',7X,'Ven',7X,'Vee',7X, + 'Vex',6X,'Vtot',6X,'Etot') 35 FORMAT (2X,A5,3X,I2,1X,6(1X,F9.3)) C IF(ITER .EQ. NITER) THEN WRITE(OUNIT,40) WRITE(OUNIT,60) ENDIF 40 FORMAT (4X,' radius ',4X,'1s density',5X,'2s density', + 5X,'2p density',4X,'total density') 60 FORMAT (4X,'-----------',4X,'----------',5X,'----------', + 5X,'----------',4X,'-------------') c C output of densities IF(ITER.EQ.NITER) THEN DO 120 IR=0,NR R=IR*DR DO 110 STATE=1,NSTATE+1 IF (STATE .EQ. NSTATE+1) THEN DEN(STATE)=RHO(IR) ELSE DEN(STATE)=NOCC(STATE)*PSTOR(IR,STATE)**2 END IF 110 CONTINUE WRITE (OUNIT,70) R,(DEN(STATE),STATE=1,NSTATE+1) 120 CONTINUE ENDIF C 70 FORMAT (5(4X,1PE11.3)) RETURN END