C -*- FORTRAN -*-
C..FILE 'common_const.h'
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
C..FILE 'common_const.h'
 
      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),IZ            
     $                                                                  
C..   INCLUDE 'common_const.h'                                          
      DIMENSION IZEFF(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/
      DATA XNQEFF /1.D0,2.D0,3.D0,3.7D0,4.0D0,4.2D0,4.3D0/
      DIMENSION ZEFF(17),XNEFF(17)
      WRITE(6,*)
      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,*)
      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
 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
      CALL ZEFCAL(ZEFF)
C
      BOHR=0.529177249D0
C
      WRITE(6,9903)
      HARTREE=-0.5D0
      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),XNEFF(I)*XNEFF(I)/ZEFF(I),
     -       HARTREE*ZEFF(I)*ZEFF(I)/XNEFF(I)/XNEFF(I)
      END DO
 9902 FORMAT(4X,A2,8X,I2,6X,F12.6,F12.2,4X,F14.6,4X,F14.6)
 9903 FORMAT(2X,'Shell',4X,'Occupation',3X,'effective Z',5X
     $     ,'effective n',5X,'radius (u.a.)',5X,'energy (u.a.)')
C
      WRITE(6,*)
      WRITE(6,*) ' the same in Angstrom and ElectronVolts '
      WRITE(6,*)
      WRITE(6,9904)
      EV=-13.6D0
      DO I=1,17
       IF (IOCC(I).NE.0) WRITE(6,9902) SHSYM(I),IOCC(I),ZEFF(I)
     $      ,XNEFF(I),BOHR*XNEFF(I)*XNEFF(I)/ZEFF(I),
     -       EV*ZEFF(I)*ZEFF(I)/XNEFF(I)/XNEFF(I)
      END DO
 9904 FORMAT(2X,'Shell',4X,'Occupation',3X,'effective Z',5X
     $     ,'effective n',5X,'radius (A)',5X,'energy (eV)')
C
C calculate the total energy of a configuration
C
      CALL ECALC(ETOTAL,ZEFF,XNEFF)
 
      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
      RADIUS=BOHR*XNEFF(IOUTER)*XNEFF(IOUTER)/ZEFF(IOUTER)
C
      WRITE(6,9928) RADIUS
 9928 FORMAT(/,'  ATOMIC RADIUS OF THE ATOM (in Angstrom): ',F12.8,/)
 
      WRITE(6,*)
      WRITE(6,*) ' ionization energies : '
      WRITE(6,*)
C
C we start with the total energy
C
      ERECAL=0.D0
      NELEC=0
      DO I=1,17
       NELEC=NELEC+IOCC(I)
      END DO
      NELTOT=NELEC
 
      IF (NELEC.GT.1) THEN
       ELAST=ETOTAL
       WRITE(6,*) ' ETOTAL (eV) = ',-ETOTAL*13.6
       ILAST=IOUTER
 101   CONTINUE
C we exclude one electron
       IF (IOCC(IOUTER).GT.1) THEN
        IOCC(IOUTER)=IOCC(IOUTER)-1
        CALL ZEFCAL(ZEFF)
C        WRITE(6,*) ' 1 '
        CALL ECALC(ENEW,ZEFF,XNEFF)
        EION=ENEW-ELAST
C        WRITE(6,*) ' ENEW, ELAST, EION :',NELEC,ENEW,ELAST,EION
       ELSE IF (IOCC(IOUTER).EQ.1) THEN
C here we have to take only the last orbital
        CALL ZEFCAL(ZEFF)
        EION=ZEFF(IOUTER)/XNEFF(IOUTER)
C        WRITE(6,*) ' 2 '
        EION=-EION*EION
C        WRITE(6,*) '  EION direct :',NELEC,EION
        ENEW=ELAST+EION
        IOCC(IOUTER)=IOCC(IOUTER)-1
       ELSE IF (IOCC(IOUTER).EQ.0) THEN
        INEXT=IOUTER-1
 202    CONTINUE
        IF (IOCC(INEXT).NE.0) THEN
         IOUTER=INEXT
        ELSE
         INEXT=INEXT-1
         GO TO 202
        END IF
C we continue with a new IOUTER
        IOCC(IOUTER)=IOCC(IOUTER)-1
        CALL ZEFCAL(ZEFF)
C        WRITE(6,*) ' 3 '
        CALL ECALC(ENEW,ZEFF,XNEFF)
        EION=ENEW-ELAST
C        WRITE(6,*) ' ENEW, ELAST, EION :',NELEC,ENEW,ELAST,EION
       END IF
 
       ELAST=ENEW
       NELEC=NELEC-1
 
C       WRITE(6,*) ' EION = ',NELTOT-NELEC,EION
 
       ERECAL=ERECAL+EION
       WRITE(6,9921) NELTOT-NELEC,-EION*13.6D0,-EION*0.5D0
 
 9921  FORMAT(' Ionization No ',I4,' with energy ',F15.6,' (eV)'
     -      ,F15.6,' (u.a.)' )
 
C if there is at least two electrons left we continue
       IF (NELEC.GT.0) GO TO 101
 
       WRITE(6,*)
       WRITE(6,*) ' recalculated total energy (sum of ionization energie
     $s, eV): '
C folded 1 (fixf $Revision: 1.3 $) 
     -      ,ERECAL*13.6
       WRITE(6,*)
 
      ELSE
       WRITE(6,*) ' no interest for NELEC=1 '
      END IF
 
 
 
      WRITE(6,*)
      WRITE(6,*)' all terminated '
      WRITE(6,*)
 
      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),IZ            
     $                                                                  
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 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
      SUBROUTINE ZEFCAL(ZEFF)
      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),IZ            
     $                                                                  
C..   INCLUDE 'common_const.h'                                          
      DIMENSION ZEFF(17)
 
      DO I=1,17
       ZEFF(I)=0.D0
       IF (IOCC(I).NE.0) THEN
        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)
         END IF
         IF (I.EQ.J) THEN
          Z=Z-DBLE(IOCC(J)-1)*SCREEN(I,I)
         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)
 
      RETURN
      END
 
      SUBROUTINE ECALC(ETOTAL,ZEFF,XNEFF)
      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),IZ            
     $                                                                  
C..   INCLUDE 'common_const.h'                                          
      DIMENSION ZEFF(17),XNEFF(17)
 
      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
 
C      WRITE(6,*) ' ETOTAL = ',ETOTAL
 
      RETURN
      END