C -*- FORTRAN -*- C..FILE 'common_const.h' C C DELCX debug C..DELCM NEUTRAL ATOMS C DELCN negative ion (not all correct) C DELC1 ion +1 C DELC2 ion +2 C DELC3 ion +3 C DELC4 ion +4 C DELC5 ion +5 C C C Calculates total energies of atoms after the data given by C J.C.Slater in his Phys.Rev. paper 1930: C C J.C.Slater: "Atomic Shielding Constants", Phys.Rev. 36 (1930) 57 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C P.Reinhardt, Paris, 1/2000 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PROGRAM SLATER IMPLICIT DOUBLE PRECISION (A-H,O-Z) C the screening constants 1-6s, 2-6p, 3-5d, 4-5f COMMON /VALS/ SCREEN(17,17) CHARACTER*2 SYMBAT,SHSYM CHARACTER*1 CST,C1 COMMON /SYMTAB/ XNQEFF(7),SYMBAT(0:92),SHSYM(17),CST(4),NQUANT(17) $ ,IMAP(17),NGAS(7),NSUBS(17),IOCCM(17),IOCC(17) C.. INCLUDE 'common_const.h' DIMENSION IZEFF(92),ETOT(92),EION(92) DATA SYMBAT/ 'XX','H ','He','Li','Be','B ','C ','N ','O ' - ,'F ','Ne','Na','Mg','Al','Si','P ','S ','Cl','Ar','K ' - ,'Ca','Sc','Ti','V ','Cr','Mn','Fe','Co','Ni','Cu','Zn' - ,'Ga','Ge','As','Se','Br','Kr','Rb','Sr','Y ','Zr','Nb' - ,'Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te' - ,'I ','Xe','Cs','Ba','La','Ce','Pr','Nd','Pm','Sm','Ee' - ,'Gd','Tb','Dy','Ho','Er','Tm','Y ','Lu','Hf','Ta','W ' - ,'Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At' - ,'Rn','Fr','Ra','Ac','Th','Pa','U '/ DATA CST /'S','P','D','F'/ DATA NGAS /2,10,18,36,54,86,200/ DATA NSUBS /2,2,6,2,6,10,2,6,10,14,2,6,10,14,2,6,10/ DATA SHSYM /'1s','2s','2p','3s','3p','3d','4s','4p','4d','4f','5s' $ ,'5p','5d','5f','6s','6p','7s'/ DATA IMAP /1,2,3,4,5,7,6,8,11,9,12,15,10,13,16,17,14/ DATA NQUANT /1,2,2,3,3,3,4,4,4,4,5,5,5,5,6,6,7/ C DATA XNQEFF /1.D0,2.D0,3.D0,3.7D0,4.0D0,4.2D0,4.3D0/ DATA XNQEFF /1.D0,2.D0,3.D0,4.0D0,5.0D0,6.0D0,7.0D0/ DIMENSION ZEFF(17),XNEFF(17) X=CPTIME(3) CALL DATING('SLATER ',1) WRITE(6,*) ' --------------------------------------------' WRITE(6,*) WRITE(6,*) ' S L A T E R ' WRITE(6,*) ' - calculation of atomic data' WRITE(6,*) ' after J.C.Slater Phys.Rev. 36 (1930) 57 ' WRITE(6,*) ' P. Reinhardt (Paris) 2001 ' WRITE(6,*) WRITE(6,*) ' --------------------------------------------' WRITE(6,*) CALL FILL WRITE(6,*) ' effective n: ' WRITE(6,9821) (XNQEFF(I),I=1,7) 9821 FORMAT(' n = 1 2 3 4 5 6 7',/ - ,' n_eff= ',7F6.2) WRITE(6,*) WRITE(6,*) C WRITE(6,*) ' PLEASE GIVE THE ATOMIC NUMBER OF THE ELEMENT' READ(5,*) IZ IF (IZ.GT.92) THEN WRITE(6,*) ' WE CAN TREAT ELEMENTS UP TO 92' STOP ' IZ TOO HIGH ' ELSE IF (IZ.LT.0) THEN WRITE(6,*) ' ATOMIC NUMBERS SHOULD BE POSITIVE ' STOP ' IZ NEGATIVE ' END IF WRITE(6,*) WRITE(6,*) ' ATOMIC SYMBOL OF THE ELEMENT: ',SYMBAT(IZ) C C which next rare gas? I=1 100 CONTINUE IF (IZ.LT.NGAS(I)) GO TO 200 I=I+1 GOTO 100 C 200 CONTINUE I=I-1 IF (NGAS(I).EQ.IZ) THEN WRITE(6,*) ' the atom is a rare gaz ' ELSE WRITE(6,*) ' the core corresponds to the rare gas ' $ ,SYMBAT(NGAS(I)) END IF WRITE(6,*) C C electronic configuration C WRITE(6,*) ' Electronic Configuration according', - ' to the Aufbau Principle: ' IDUM=IZ DO I=1,17 IOCCM(IMAP(I))=MIN(NSUBS(IMAP(I)),IDUM) IDUM=MAX(IDUM-NSUBS(IMAP(I)),0) IOCC(IMAP(I))=IOCCM(IMAP(I)) END DO C DO I=1,17 WRITE(6,9901) I,SHSYM(I),IOCCM(I) END DO C C positive IONS C C//SKIPCM cM/SKIPCN cMhighest shell: cM IHIGH=0 cM DO I=1,17 cM IF (IOCC(I).NE.0) IHIGH=I cM END DO cM WRITE(6,*) ' FIRST SHELL TO BE IONIZIED: ',SHSYM(IHIGH) cM IOCC(IHIGH)=IOCC(IHIGH)-1 cM/SKIPC1 cM IF (IZ.GT.1) THEN cMhighest shell: cM IHIGH=0 cM DO I=1,17 cM IF (IOCC(I).NE.0) IHIGH=I cM END DO cM WRITE(6,*) ' SECOND SHELL TO BE IONIZIED: ',SHSYM(IHIGH) cM IOCC(IHIGH)=IOCC(IHIGH)-1 cMhighest shell: cM/SKIPC2 cM IF (IZ.GT.2) THEN cM IHIGH=0 cM DO I=1,17 cM IF (IOCC(I).NE.0) IHIGH=I cM END DO cM WRITE(6,*) ' THIRD SHELL TO BE IONIZIED: ',SHSYM(IHIGH) cM IOCC(IHIGH)=IOCC(IHIGH)-1 cMhighest shell: cM/SKIPC3 cM IF (IZ.GT.3) THEN cM IHIGH=0 cM DO I=1,17 cM IF (IOCC(I).NE.0) IHIGH=I cM END DO cM WRITE(6,*) ' FOURTH SHELL TO BE IONIZIED: ',SHSYM(IHIGH) cM IOCC(IHIGH)=IOCC(IHIGH)-1 cMhighest shell: cM/SKIPC4 cM IF (IZ.GT.4) THEN cM IHIGH=0 cM DO I=1,17 cM IF (IOCC(I).NE.0) IHIGH=I cM END DO cM WRITE(6,*) ' FIFTH SHELL TO BE IONIZIED: ',SHSYM(IHIGH) cM IOCC(IHIGH)=IOCC(IHIGH)-1 cM END IF cM/ENDC4 cM END IF cM/ENDC3 cM END IF cM/ENDC2 cM END IF cM/ENDC1 cM/ENDCN cM cM cMnegative IONS cM cM/SKIPC1 cM/SKIPC2 cM/SKIPC3 cM/SKIPC4 cM/SKIPC5 cMhighest shell: cM IHIGH=0 cM DO I=1,17 cM IF (IOCC(I).NE.0) THEN cM IF (IOCC(I).EQ.NSUBS(I)) THEN cM IHIGH=I+1 cM ELSE cM IHIGH=I cM END IF cM END IF cM END DO cM WRITE(6,*) ' WE ADD ONE ELECTRON TO SHELL: ',SHSYM(IHIGH) cM IOCC(IHIGH)=IOCC(IHIGH)+1 cM/ENDC5 cM/ENDC4 cM/ENDC3 cM/ENDC2 cM/ENDC1 cM/ENDCM C 9901 FORMAT(5X,' Shell ',I2,' (',A2,'): ',I2,' electrons') WRITE(6,*) C C this is the occupation according to the Aufbau principle C now we reorder according to the main quantum number n C 115 CONTINUE WRITE(6,*) ' Do you want to modify the configuration? (y/n)' READ(5,'(A1)') C1 IF (C1.EQ.'Y'.OR.C1.EQ.'y') THEN 117 CONTINUE WRITE(6,*) ' Which shell occupation should be modified? ' READ(5,*) IACT 118 CONTINUE WRITE(6,*) ' How many electrons should be there? ' READ(5,*) INEW IF (INEW.GT.NSUBS(IACT)) THEN WRITE(6,*) ' THIS SHELL CAN ONLY HOLD ',NSUBS(IACT) $ ,' ELECTRONS! ' GO TO 118 END IF IOCC(IACT)=INEW WRITE(6,*) ' Is that all? (y/n)' READ(5,'(A1)') C1 IF (C1.EQ.'N'.OR.C1.EQ.'n') GO TO 117 C print actual configuration WRITE(6,*) WRITE(6,*) ' Actual electronic configuration:' IELEC=0 DO I=1,17 IF (IOCC(I).NE.0) THEN IELEC=IELEC+IOCC(I) WRITE(6,9901) I,SHSYM(I),IOCC(I) END IF END DO WRITE(6,*) WRITE(6,*) ' total charge: ',IZ-IELEC WRITE(6,*) goto 115 END IF C WRITE(6,*) WRITE(6,*) ' chosen electronic configuration:' IELEC=0 DO I=1,17 IF (IOCC(I).NE.0) THEN WRITE(6,9901) I,SHSYM(I),IOCC(I) IELEC=IELEC+IOCC(I) END IF END DO WRITE(6,*) WRITE(6,*) ' total charge: ',IZ-IELEC WRITE(6,*) WRITE(6,*) WRITE(6,*) ' Calculation of the total energy: ' WRITE(6,*) C C calculate the Zeff C DO I=1,17 ZEFF(I)=0.D0 IF (IOCC(I).NE.0) THEN CX WRITE(6,*) ' TREATING SHELL ',I,': (',SHSYM(I),')' Z=DBLE(IZ) DO J=1,I IF (I.NE.J.AND.IOCC(J).NE.0) THEN Z=Z-DBLE(IOCC(J))*SCREEN(I,J) CX WRITE(6,*) ' TREATING COUPLE ',I,J,': (' CX $ ,SHSYM(I),',',SHSYM(J),') :',SCREEN(I CX $ ,J),IOCC(J),Z END IF IF (I.EQ.J) THEN Z=Z-DBLE(IOCC(J)-1)*SCREEN(I,I) CX WRITE(6,*) ' TREATING COUPLE ',I,J,': (' CX $ ,SHSYM(I),',',SHSYM(J),') :',SCREEN(I CX $ ,J),IOCC(J)-1,Z END IF END DO ZEFF(I)=Z END IF END DO C C correct for sp shells: C if p, then Zeff(s)=Zeff(p) C if no p, then no correction C C this concerns shells 2/3 4/5 7/8 11/12 15/16 IF (IOCC(3) .NE.0) ZEFF(2) =ZEFF(3) IF (IOCC(5) .NE.0) ZEFF(4) =ZEFF(5) IF (IOCC(8) .NE.0) ZEFF(7) =ZEFF(8) IF (IOCC(12).NE.0) ZEFF(11)=ZEFF(12) IF (IOCC(16).NE.0) ZEFF(15)=ZEFF(16) C BOHR=0.529177249D0 C WRITE(6,9903) HARTREE=-13.6D0 HARTREE=-1.D0 DO I=1,17 XNEFF(I)=XNQEFF(NQUANT(I)) IF (IOCC(I).NE.0) WRITE(6,9902) SHSYM(I),IOCC(I),ZEFF(I) $ ,XNEFF(I),BOHR*XNEFF(I)*XNEFF(I)/ZEFF(I), - 0.5d0*HARTREE*ZEFF(I)*ZEFF(I)/XNEFF(I)/XNEFF(I), - -13.6d0*ZEFF(I)*ZEFF(I)/XNEFF(I)/XNEFF(I) END DO 9902 FORMAT(4X,A2,8X,I2,6X,F12.6,F12.2,6X,F12.6,6X,F12.6,F14.2) 9903 FORMAT(2X,'Shell',4X,'Occupation',3X,'effective Z',5X $ ,'effective n',4X,'radius (A)',4X,'energy (a.u.)' - ,4X,' (eV)') C C calculate the total energy of a configuration C ETOTAL=0.D0 DO I=1,17 IF (IOCC(I).NE.0) ETOTAL=ETOTAL+DBLE(IOCC(I))*(ZEFF(I)/XNEFF(I)) $ *(ZEFF(I)/XNEFF(I)) END DO EEV=-ETOTAL*13.6D0 EAU=-ETOTAL*0.5D0 ECM=-ETOTAL*109700D0 EKJ=-ETOTAL*1311.04D0 WRITE(6,*) WRITE(6,*) ' TOTAL ENERGY (in 1/2 a.u.) : ',-ETOTAL WRITE(6,*) WRITE(6,*) ' TOTAL ENERGY (in eV): ',EEV WRITE(6,*) ' TOTAL ENERGY (in a.u.): ',EAU WRITE(6,*) ' TOTAL ENERGY (in cm-1): ',ECM WRITE(6,*) ' TOTAL ENERGY (in kJ/mol): ',EKJ C C calculate the atomic radius C which is the outermost shell? IOUTER=17 107 CONTINUE IF (IOCC(IOUTER).EQ.0) THEN IOUTER=IOUTER-1 GO TO 107 END IF C RADIUS=BOHR*XNEFF(IOUTER)*XNEFF(IOUTER)/ZEFF(IOUTER) WRITE(6,9928) RADIUS 9928 FORMAT(/,' RADIUS OF THE ATOM (in Angstrom): ',F12.8,/) X=CPTIME(4) CALL DATING('SLATER ',2) END C SUBROUTINE FILL IMPLICIT DOUBLE PRECISION (A-H,O-Z) C the screening constants 1-6s, 2-6p, 3-5d, 4-5f COMMON /VALS/ SCREEN(17,17) CHARACTER*2 SYMBAT,SHSYM CHARACTER*1 CST,C1 COMMON /SYMTAB/ XNQEFF(7),SYMBAT(0:92),SHSYM(17),CST(4),NQUANT(17) $ ,IMAP(17),NGAS(7),NSUBS(17),IOCCM(17),IOCC(17) C.. INCLUDE 'common_const.h' C C parametrization Z30=0.31D0 Z35=0.35D0 Z85=0.85D0 C WRITE(6,*) $ ' the Slater model is used with the following constants ' WRITE(6,*) WRITE(6,*) ' screening constants:' WRITE(6,*) ' 1s - 1s: ',Z30 WRITE(6,*) ' ns - ns: ',Z35 WRITE(6,*) ' ns - (n-1)s: ',Z85 WRITE(6,*) C C 1s, 2s, 2p, 3s, 3p, 3d, 4s, 4p, 4d, 4f, 5s, 5p, 5d, 5f, 6s, 6p, 7s C 1s .31 C 2s .85 .35 C 2p .85 .35 .35 C 3s 1 .85 .85 .35 C 3p 1 .85 .85 .35 .35 C 3d 1 1 1 1 1 .35 C 4s 1 1 1 .85 .85 .85 .35 C 4p 1 1 1 .85 .85 .85 .35 .35 C 4d 1 1 1 1 1 1 1 1 .35 C 4f 1 1 1 1 1 1 1 1 1 .35 C 5s 1 1 1 1 1 1 .85 .85 .85 .85 .35 C 5p 1 1 1 1 1 1 .85 .85 .85 .85 .35 .35 C 5d 1 1 1 1 1 1 1 1 1 1 1 1 .35 C 5f 1 1 1 1 1 1 1 1 1 1 1 1 1 .35 C 6s 1 1 1 1 1 1 1 1 1 1 .85 .85 .85 .35 .35 C 6p 1 1 1 1 1 1 1 1 1 1 .85 .85 .85 .85 .35 .35 C 7s 1 1 1 1 1 1 1 1 1 1 1 1 1 1 .85 .85 .35 C DO I=1,17 DO J=I+1,17 SCREEN(I,J)=0.D0 END DO DO J=1,I SCREEN(I,J)=1.D0 END DO END DO C C 1s - 1s SCREEN (1,1)=Z30 C the remainder of the diagonal DO I=2,17 SCREEN(I,I)=Z35 END DO C 2s,26 - 1s SCREEN (2,1)=Z85 SCREEN (3,1)=Z85 C 2p,3s,3p - 2s SCREEN (3,2)=Z35 SCREEN (4,2)=Z85 SCREEN (5,2)=Z85 C 3s,3p - 2p SCREEN (4,3)=Z85 SCREEN (5,3)=Z85 C 3p,4s,4p - 3s SCREEN (5,4) =Z35 SCREEN (7,4) =Z85 SCREEN (8,4) =Z85 C 4s,4p - 3p SCREEN (7,5) =Z85 SCREEN (8,5) =Z85 C 4s,4p - 3d SCREEN (7,6) =Z85 SCREEN (8,6) =Z85 C 4p,5s,5p - 4s SCREEN (8,7) =Z35 SCREEN (11,7)=Z85 SCREEN (12,7)=Z85 C 5s,5p - 4p SCREEN (11,8)=Z85 SCREEN (12,8)=Z85 C 5s,5p - 4d SCREEN (11,9)=Z85 SCREEN (12,9)=Z85 C 5s,5p - 4f SCREEN (11,10)=Z85 SCREEN (12,10)=Z85 C 5p,6s,6p - 5s SCREEN (12,11)=Z35 SCREEN (15,11)=Z85 SCREEN (16,11)=Z85 C 6s,6p - 5p SCREEN (15,12)=Z85 SCREEN (16,12)=Z85 C 6s,6p - 5d SCREEN (15,13)=Z85 SCREEN (16,13)=Z85 C 6s,6p - 5f SCREEN (15,14)=Z85 SCREEN (16,14)=Z85 C 6p,7s - 6s SCREEN (16,15)=Z35 SCREEN (17,15)=Z85 C 7s - 6p SCREEN (17,16)=Z85 C DO I=1,17 DO J=I+1,17 SCREEN(I,J)=SCREEN(J,I) END DO END DO RETURN END C C..DELCF WE HAVE A ROUTINE FLU SUBROUTINE CPUTIM(TIM) IMPLICIT REAL*8 (A-H,O-Z) C C UNIX SYSTEM C REAL T(2),SEC1,ETIME SEC1=ETIME(T) TIM=DBLE(T(1)+T(2)) C C dummy C C TIM=0. C C AIX/PS2 C C INTEGER CLOCK C I=CLOCK(I) C TIM=DBLE(I/1000000.) C C RS6000 C C REAL T C T=MCLOCK(TIMES) C TIM=T/100. RETURN END C REAL*8 FUNCTION CPTIME(I) IMPLICIT REAL*8 (A-H,O-Z) INTEGER I C IF (I.EQ.1) THEN CALL CPUTIM(X) CPTIME=X ELSEIF (I.EQ.4.OR.I.EQ.3) THEN CALL CPUTIM(X) IX1=INT(X/3600.) IX2=INT((X-REAL(IX1*3600))/60.) X2=X-REAL(IX1*3600 + IX2*60) IF (I.EQ.4) THEN CALL TIMING(' ALL ') WRITE(6,*) WRITE (6,9991) IX1,IX2,X2 END IF 9991 FORMAT(' TOTAL CPUTIME: ',I5,' hrs ',I3,' mins and ',f7.4,' secs' $) C folded 1 (fixf $Revision: 1.3 $) CPTIME=X ENDIF RETURN END C SUBROUTINE SECNND(X) REAL*8 X,CPTIME X=CPTIME(1) RETURN END C SUBROUTINE TIMG(NOMZ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*4 NOMZ CALL SECNND(FT) WRITE(6,100)NOMZ,FT 100 FORMAT(1X,53('T'),1X,A4,' CPU',F11.3) RETURN END C SUBROUTINE TIMING(NOMZ) IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*4 NOMZ CALL SECNND(FT) WRITE(6,100)NOMZ,FT 100 FORMAT(1X,53('T'),1X,A4,' CPU',F11.3) CALL FLUSH(6) RETURN END C C -*- Fortran -*- SUBROUTINE DATIM(A,B) CHARACTER*8 A,B CHARACTER*24 DATE CHARACTER*3 MON(12) CHARACTER*2 M(13) DATA MON /'Jan','Feb','Mar','Apr','May','Jun','Jul', * 'Aug','Sep','Oct','Nov','Dec'/ DATA M /'01','02','03','04','05','06','07','08', * '09','10','11','12','**'/ DO 3 I=1,24 3 DATE(I:I)='-' CALL FDATE(DATE) DO 1 I=1,12 IF (DATE(5:7).EQ.MON(I)) GOTO 2 1 CONTINUE I=13 2 CONTINUE A(1:2)=DATE(9:10) A(3:3)='.' A(4:5)=M(I) A(6:6)='.' A(7:8)=DATE(23:24) B(1:8)=DATE(12:19) RETURN END C SUBROUTINE DATING(PNAME,IACT) CHARACTER*8 PNAME,A,B CALL DATIM(A,B) C remove leading or trailing blanks from pname IFST=1 100 CONTINUE IF (PNAME(IFST:IFST).EQ.' ') THEN IFST=IFST+1 IF (IFST.LT.8) GO TO 100 END IF ILST=8 IF (IFST.EQ.8) GO TO 300 200 CONTINUE IF (PNAME(ILST:ILST).EQ.' ') THEN ILST=ILST-1 IF (ILST.GT.1) GO TO 200 END IF 300 CONTINUE IF (IACT.EQ.1) THEN WRITE(6,9901) PNAME(IFST:ILST),A,B ELSE IF (IACT.EQ.2) THEN WRITE(6,9902) PNAME(IFST:ILST),A,B ELSE WRITE(6,9903) PNAME(IFST:ILST),A,B END IF 9901 FORMAT(//,' Started ',A,'; date: ',A8,' time: ',A8,//) 9902 FORMAT(//,' Finished ',A,'; date: ',A8,' time: ',A8,//) 9903 FORMAT(//,1X,A,'; date: ',A8,' time: ',A8,//) RETURN END C