atom_input.f 38.4 KB
Newer Older
1 2 3 4 5 6
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

7 8
      MODULE m_atominput
      use m_juDFT
9
      INTEGER, PARAMETER ::  l_buffer=512   ! maximum length of e-config string
10
      INTEGER, PARAMETER  :: dbgfh=6, errfh=6, warnfh=6
11 12 13 14 15 16 17 18
      REAL, PARAMETER     :: eps=0.00000001

      CONTAINS
!***********************************************************************
!     reads in the parameters associated with atom types from input
!     file. Part of inp-generator
!***********************************************************************
      SUBROUTINE atom_input(
19
     >                      infh,xl_buffer,bfh,buffer,
20
     >                      input,idlist,xmlCoreRefOccs,
21
     X                      nline,xmlElectronStates,
22
     X                      xmlPrintCoreStates,xmlCoreOccs,
23
     >                      atomTypeSpecies,numSpecies,
24 25 26 27
     <                      nel,atoms,enpara )

      USE m_types
      USE m_juDFT_init
28 29
      USE m_readrecord
      USE m_setatomcore, ONLY : setatom_bystr, setcore_bystr
30
      USE m_constants
31
    
32 33
      IMPLICIT NONE

34
      TYPE(t_input),INTENT(INOUT)    :: input
35
      TYPE(t_enpara),INTENT(OUT)     :: enpara
36 37
      TYPE(t_atoms),INTENT(INOUT)    :: atoms

38 39
! ... Arguments ...
      INTEGER, INTENT (IN)    :: infh  ! file number of input-file
40
      INTEGER, INTENT (IN)    :: bfh
41
      INTEGER, INTENT (IN)    :: numSpecies
42 43 44
      INTEGER, INTENT (INOUT) :: nline ! current line in this file
      INTEGER, INTENT (INOUT) :: nel   ! number of valence electrons

45
      INTEGER, INTENT (IN)     :: xl_buffer
46
      INTEGER, INTENT (IN)     :: atomTypeSpecies(atoms%ntype)
47
      REAL   , INTENT (IN)     :: idlist(atoms%ntype)
48
      REAL   , INTENT (IN)     :: xmlCoreRefOccs(29)
49
      REAL, INTENT (INOUT)     :: xmlCoreOccs(2,29,atoms%ntype)
50
      INTEGER, INTENT (INOUT)  :: xmlElectronStates(29,atoms%ntype)
51
      LOGICAL, INTENT (INOUT)  :: xmlPrintCoreStates(29,atoms%ntype)
52 53 54 55 56 57 58 59 60 61 62
      CHARACTER(len=xl_buffer) :: buffer

!===> data
      INTEGER, PARAMETER ::  nwdd=2         ! maximum number of windows
      INTEGER, PARAMETER ::  nstd=31        ! maximum number of core states

!===> Local Variables
      INTEGER :: nbuffer,ios,n,i,j,l,d1,d10,aoff,up,dn
      INTEGER :: lmax0_def,lnonsph0_def,jri0_def,ncst0_def
      INTEGER :: lmax0,lnonsph0,jri0,ncst0,nlod0,llod
      INTEGER :: natomst,ncorest,nvalst,z,nlo0
63
      INTEGER :: xmlCoreStateNumber, lmaxdTemp
64
      REAL    :: rmt0_def,dx0_def,bmu0_def, upReal, dnReal
65
      REAL    :: rmt0,dx0,bmu0,zat0,id,electronsOnAtom, electronsLeft
66
      LOGICAL :: fatalerror, h_atom, h_allatoms
67 68 69 70 71
      LOGICAL :: idone(atoms%ntype) 
      INTEGER :: lonqn(atoms%nlod,atoms%ntype),z_int(atoms%ntype)
      INTEGER :: coreqn(2,nstd,atoms%ntype),lval(nstd,atoms%ntype)
      INTEGER :: llo0(atoms%nlod)
      REAL    :: nelec(0:nwdd),coreocc(nstd,atoms%ntype)
72 73 74


      CHARACTER(len= l_buffer) :: econfig0_def,econfig0
75 76
      CHARACTER(len= l_buffer) :: econfig(atoms%ntype) ! verbose electronic config
      CHARACTER(len=80) :: lo(atoms%ntype), lo0
77
      CHARACTER(len=13) :: fname
78
      CHARACTER(LEN=20) :: speciesName0
79 80

      CHARACTER(len=1) :: lotype(0:3)
81

82 83 84 85 86 87 88
      DATA lotype /'s','p','d','f'/

!---> initialize some variables

      fatalerror = .false.
      h_atom=.false.;h_allatoms=.false.

89 90 91 92
      idone(1:atoms%ntype) = .false.
      z_int(1:atoms%ntype) = NINT(atoms%zatom(1:atoms%ntype))
      lo = ' ' ; nlod0 = 0 ; atoms%nlo = 0 ; llod = 0 ; lonqn = 0
      atoms%ncst = 0 ; econfig(1:atoms%ntype) = ' '
93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
!
      lmax0_def    = -9999  
      lnonsph0_def = -9999
      rmt0_def     = -9999.9
      dx0_def      = -9999.9
      jri0_def     = -9999
      ncst0_def    = -9999
      econfig0_def = 'NONE'
      bmu0_def     = -9999.9

      WRITE(6,*)
      WRITE(6,'(a50)') '==============================================='
      WRITE(6,'(a50)') '===  modifying atomic input for &(all)atom  ==='
      WRITE(6,'(a50)') '==============================================='
      WRITE(6,*)

!===> continue reading input
      
      nbuffer = len_trim(buffer)

      IF ((buffer(1:9)=='&allatoms') .OR. 
     &    (buffer(1:5)=='&atom') .OR.
!    resetting nbuffer for &qss or &soc interferes with the lapw
!    namelist, therefore, its contributions are also checked for.
!    might interfere with other namelists, too.
!    Klueppelberg Jul 2012
     &    (buffer(1:5)=='&comp') .OR.
     &    (buffer(1:5)=='&exco') .OR.
     &    (buffer(1:5)=='&film') .OR.
     &    (buffer(1:4)=='&kpt') ) THEN
      ELSE
        nbuffer = 0 ! reset, to read in after &qss or &soc 
      ENDIF
      
      loop: DO

      IF (nbuffer == 0) THEN
        DO
131
          CALL read_record(infh,xl_buffer,bfh,nline,nbuffer,buffer,ios)
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155
          IF (ios==1) GOTO 999
          IF (ios == 2)  CALL juDFT_error
     +         ("end of file while reading a record",
     +         calledby ="atom_input")
          IF (buffer(1:1)=='&') EXIT
          CALL err(0)
          fatalerror = .true.
        ENDDO
      ENDIF

!===> allatoms

      IF (buffer(1:9)=='&allatoms') THEN
        IF (h_allatoms) CALL err(1)
        h_allatoms = .true.
        IF (h_atom) then
          WRITE (errfh,*)
          WRITE (errfh,*) 'atom_input: ERROR',
     &     'namelist &allatoms must appear before namelist(s) &atom.'
          WRITE (errfh,*)
          fatalerror = .true.
        ELSE
!--->     read defaults for atom defaults
          CALL read_allatoms(
156
     >                       bfh,
157 158 159 160 161 162
     <                       rmt0_def,dx0_def,jri0_def,lmax0_def,
     <                       lnonsph0_def,ncst0_def,econfig0_def,
     <                       bmu0_def,ios)

          IF (ios.NE.0) GOTO 912
          IF (rmt0_def > -9999.8) THEN
163
            atoms%rmt     = rmt0_def
164 165 166
            WRITE (6,'(a25,f12.6)') 'globally changed rmt to',rmt0_def
          ENDIF
          IF (dx0_def  > -9999.8)   THEN
167
            atoms%dx      = dx0_def
168 169 170
            WRITE (6,'(a25,f12.6)') 'globally changed dx  to',dx0_def
          ENDIF
          IF (jri0_def > -9998  )   THEN
171
            atoms%jri     = jri0_def
172 173 174
            WRITE (6,'(a25,i12)') 'globally changed jri to',jri0_def
          ENDIF
          IF (lmax0_def > -9998 )   THEN
175
            atoms%lmax    = lmax0_def
176 177 178 179
            WRITE (6,'(a26,i12)') 'globally changed lmax to',
     &                                                       lmax0_def
          ENDIF
          IF (lnonsph0_def > -9998) THEN
180
            atoms%lnonsph = lnonsph0_def
181 182 183 184
            WRITE (6,'(a28,i12)') 'globally changed lnonsph to ',
     &                                                    lnonsph0_def
          ENDIF
          IF (ncst0_def > -9998 )   THEN
185
            atoms%ncst    = ncst0_def
186 187 188 189 190 191 192 193 194
            WRITE (6,'(a26,i12)') 'globally changed ncst to',
     &                                                       ncst0_def
          ENDIF
          IF (econfig0_def.NE.'NONE') THEN
            econfig = econfig0_def
            WRITE (6,'(a26,a80)') 'globally set econfig to ',
     &                                                    econfig0_def
          ENDIF
          IF (bmu0_def > -9999.8)   THEN
195
            atoms%bmu     = bmu0_def
196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213
            WRITE (6,'(a25,f12.6)') 'globally changed bmu to',bmu0_def
          ENDIF
        ENDIF

!===> atom
      ELSEIF (buffer(1:5)=='&atom') THEN
        h_atom=.true.

!--->   set atom defaults
        lmax0    = -9999  
        lnonsph0 = -9999
        rmt0     = -9999.9
        dx0      = -9999.9
        jri0     = -9999
        ncst0    = -9999
        econfig0 = 'NONE'
        bmu0     = -9999.9
        lo0      = ' '
214
        speciesName0 = ''
215 216 217

!--->   read namelist
        CALL read_atom(
218
     >                 bfh,lotype,
219
     <                 id,zat0,rmt0,jri0,dx0,lmax0,lnonsph0,
220 221
     <                 ncst0,econfig0,speciesName0,bmu0,lo0,nlod0,llod,
     <                 ios)
222 223 224 225
        IF (ios.ne.0) THEN
          CALL err(3)
        ELSE
!--->     put the data into the correct place
226
          DO n = 1, atoms%ntype
227 228 229 230 231 232
            IF (abs( id - idlist(n) ) > 0.001) CYCLE
            IF (idone(n)) then
              WRITE (errfh,*) 'atom_input: ERROR. did that one already'
              fatalerror=.true.
              EXIT
            ELSE
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250
              IF (speciesName0.NE.'') THEN
                 atoms%speciesName(atomTypeSpecies(n)) = 
     &              TRIM(ADJUSTL(speciesName0))
                 DO i = 1, numSpecies
                    IF (i.NE.atomTypeSpecies(n)) THEN
                       IF((TRIM(ADJUSTL(speciesName0))).EQ.
     &                    (TRIM(ADJUSTL(atoms%speciesName(i))))) THEN
                          WRITE(*,*) ''
                          WRITE(*,*) 'Error for species name'
                          WRITE(*,*) TRIM(ADJUSTL(speciesName0))
                          WRITE(*,*) ''
                          CALL juDFT_error
     +                       ("Same name for different species",
     +                        calledby ="atom_input")
                       END IF
                    END IF
                 END DO
              END IF
251
              IF (rmt0 > -9999.8) THEN
252
                atoms%rmt(n)  = rmt0
253 254 255 256
                WRITE (6,'(a9,i4,2a2,a16,f12.6)') 'for atom ',n,
     &                ' (',namat_const(z_int(n)),') changed rmt to',rmt0
              ENDIF
              IF (dx0 > -9999.8) THEN
257
                atoms%dx(n)  = dx0
258 259 260 261
                WRITE (6,'(a9,i4,2a2,a16,f12.6)') 'for atom ',n,
     &                ' (',namat_const(z_int(n)),') changed dx  to', dx0
              ENDIF
              IF (jri0 > -9998  ) THEN
262
                atoms%jri(n)  = jri0
263 264 265 266
                WRITE (6,'(a9,i4,2a2,a16,i12)') 'for atom ',n,
     &                ' (',namat_const(z_int(n)),') changed jri to',jri0
              ENDIF
              IF (lmax0 > -9998  ) THEN
267
                atoms%lmax(n)  = lmax0
268 269 270 271
                WRITE (6,'(a9,i4,2a2,a17,i12)') 'for atom ',n,
     &              ' (',namat_const(z_int(n)),') changed lmax to',lmax0
              ENDIF
              IF (lnonsph0 > -9998  ) THEN
272
                atoms%lnonsph(n)  = lnonsph0
273 274 275 276
                WRITE (6,'(a9,i4,2a2,a20,i12)') 'for atom ',n,
     &        ' (',namat_const(z_int(n)),') changed lnonsph to',lnonsph0
              ENDIF
              IF (bmu0 > -9999.8  ) THEN
277
                atoms%bmu(n)  = bmu0
278 279 280 281
                WRITE (6,'(a9,i4,2a2,a16,f12.6)') 'for atom ',n,
     &              ' (',namat_const(z_int(n)),  ') changed bmu to',bmu0
              ENDIF
              IF (ncst0 > -9998  ) THEN
282
                atoms%ncst(n)  = ncst0
283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302
                WRITE (6,'(a9,i4,2a2,a17,i12)') 'for atom ',n,
     &             ' (',namat_const(z_int(n)), ') changed ncst to',ncst0
              ENDIF
! ===> electronic configuration
              IF (econfig0.NE.'NONE') THEN
                 econfig(n) = econfig0
                 WRITE (6,'(a9,i4,2a2,a17,a80)') 'for atom ',n,
     &           ' (',namat_const(z_int(n)),') set econfig to ',econfig0
                 CALL setatom_bystr(
     >                              l_buffer,nwdd,econfig(n),
     <                          natomst,ncorest,nvalst,nelec)
                 WRITE (6,'("   corestates =",i3," with",f6.1,
     &                                 " electrons")')  ncorest,nelec(0)
                 WRITE (6,'("   valence st.=",i3," with",f6.1,
     &                                 " electrons")')   nvalst,nelec(1)
                 IF (nelec(2) /= 0) THEN
                 WRITE (6,'("second window found!")')
                 WRITE (6,'("   valence st.=",i3," with",f6.1,
     &                                 " electrons")')   nvalst,nelec(2)
                 ENDIF
303 304
                 IF (nelec(0)+nelec(1)+nelec(2)-
     &               atoms%zatom(n)>0.01) THEN
305 306 307 308 309 310 311
                    CALL juDFT_error
     +                   ("econfig does not fit to this atom type!"
     +                   ,calledby ="atom_input")
                 ENDIF
                 IF (ncst0 > -9998  ) THEN
                   IF (ncorest /= ncst0) THEN
                     WRITE (6,'("  ==> core-states (ncst):",i3,
312 313
     &                              " =/= (econfig):",i3)') 
     &                        atoms%ncst,ncorest
314 315 316 317 318
                     CALL juDFT_error
     +                    ("econfig does not fit to the specified ncst"
     +                    ,calledby ="atom_input")
                   ENDIF
                 ELSE
319
                   atoms%ncst(n) = ncorest
320 321 322 323
                 ENDIF
              ENDIF
! ===> local orbitals
              IF (lo0 /= ' ') THEN
324
                WRITE (6,'(a6,i3,a7,i3,a3,a80)')
325 326
     &                     "nlod =",nlod0," llod =",llod," : ",lo0
                lo(n)      = lo0
327 328
                IF (nlod0 > atoms%nlod)  
     &             CALL juDFT_error("atom_input: too "
329 330
     &                              //"many lo",calledby="atom_input")

331 332
                atoms%nlo(n) = len_trim(lo(n))/2
                DO i = 1, atoms%nlo(n)
333 334 335
                  j = 2*i
                  DO l = 0, 3
                    IF (lo(n)(j:j) == lotype(l)) THEN
336
                      atoms%llo(i,n) = l
337 338 339 340 341 342
                    ENDIF
                  ENDDO
                  j = j - 1
                  READ (lo(n)(j:j),*) lonqn(i,n)
                ENDDO
                WRITE (6,'("   nlo(",i3,") = ",i2," llo = ",8i2)') n,
343 344 345
     &                    atoms%nlo(n),(atoms%llo(i,n),i=1,atoms%nlo(n))
                WRITE (6,'("   lonqn = ",8i2)') 
     &                    (lonqn(i,n),i=1,atoms%nlo(n))
346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366
              ENDIF
              idone(n)   = .true.
            ENDIF
          ENDDO
        ENDIF

!===> not an atom related namelist, we are done
      ELSE
        exit loop
      ENDIF

      nbuffer = 0
      ENDDO loop

 999  CONTINUE
      IF (fatalerror) 
     &   CALL juDFT_error("ERROR(S) reading input. Check output for "
     &                  //"details.",calledby="atom_input")

!----------- adjust the core-levels, lo's and the energy parameters ----

367 368
      coreqn(1:2,1:nstd,1:atoms%ntype) = 0
      coreocc(1:nstd,1:atoms%ntype) = -1.0
369

370
      nel = 0
371 372 373

      IF ( ANY(atoms%bmu(:) > 0.0) ) input%jspins=2 

374 375
      lmaxdTemp = atoms%lmaxd
      atoms%lmaxd = 3
376
      call enpara%init(atoms,input%jspins)
377
      DO n = 1, atoms%ntype
378 379

        CALL setcore_bystr(
380
     >                      n,nstd,atoms%ntype,l_buffer,
381 382 383 384 385 386 387 388 389 390 391 392
     X                      econfig,natomst,ncorest,
     <                      coreqn,coreocc)

        IF ( coreqn(1,1,n) /= 0 ) THEN
          DO i = 1, natomst
            IF (coreqn(2,i,n) < 0) THEN
               lval(i,n) = - coreqn(2,i,n) - 1
            ELSE
               lval(i,n) = coreqn(2,i,n)
            ENDIF
          ENDDO 

393 394
           d1  = mod(nint(atoms%zatom(n)),10)
           d10 = int( (nint(atoms%zatom(n)) + 0.5)/10 )
395
          aoff = iachar('1')-1
396 397 398 399 400
          IF(.NOT.input%l_inpXML) THEN
             fname = 'corelevels.'//achar(d10+aoff)//achar(d1+aoff)
             OPEN (27,file=fname,form='formatted')
             write(27,'(i3)') natomst
          END IF
401 402

          WRITE (6,*) '----------'
403
          electronsOnAtom = 0
404 405 406 407 408 409 410
          DO i = 1, ncorest
            WRITE(6,'("     core :",2i3,f6.1)') 
     &             coreqn(1,i,n),coreqn(2,i,n),coreocc(i,n)
            j = INT(coreocc(i,n) / 2)
            IF (coreocc(i,n) > 2*j) THEN
              j = - coreocc(i,n)
            ENDIF
411 412 413
            IF(.NOT.input%l_inpXML) THEN
               write(27,'(4i3)') coreqn(1,i,n),coreqn(2,i,n),j,j
            END IF
414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453
            xmlCoreStateNumber = 0
            SELECT CASE(coreqn(1,i,n))
               CASE (1)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 1   !(1s1/2)
               CASE (2)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 2   !(2s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 3    !(2p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 4   !(2p3/2)
               CASE (3)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 5   !(3s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 6    !(3p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 7   !(3p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 9    !(3d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 10  !(3d5/2)
               CASE (4)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 8   !(4s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 11   !(4p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 12  !(4p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 14   !(4d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 15  !(4d5/2)
                  IF(coreqn(2,i,n).EQ.3) xmlCoreStateNumber = 19   !(4f5/2)
                  IF(coreqn(2,i,n).EQ.-4) xmlCoreStateNumber = 20  !(4f7/2)
               CASE (5)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 13  !(5s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 16   !(5p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 17  !(5p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 21   !(5d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 22  !(5d5/2)
                  IF(coreqn(2,i,n).EQ.3) xmlCoreStateNumber = 26   !(5f5/2)
                  IF(coreqn(2,i,n).EQ.-4) xmlCoreStateNumber = 27  !(5f7/2)
               CASE (6)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 18  !(6s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 23   !(6p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 24  !(6p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 28   !(6d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 29  !(6d5/2)
               CASE (7)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 25  !(7s1/2)
            END SELECT
            IF(xmlCoreStateNumber.EQ.0) STOP 'Invalid core state!'
454
            xmlElectronStates(xmlCoreStateNumber,n) = coreState_const
455 456 457 458 459 460 461 462 463 464 465 466 467
            xmlPrintCoreStates(xmlCoreStateNumber,n) = 
     +         coreocc(i,n).NE.xmlCoreRefOccs(xmlCoreStateNumber)
            SELECT CASE(xmlCoreStateNumber)
               CASE (9:10,14:15,19:22,26:29)
                  up = MIN((xmlCoreRefOccs(xmlCoreStateNumber)/2),
     +                          coreocc(i,n))
                  dn = MAX(0.0,coreocc(i,n)-up)
               CASE DEFAULT
                  up = CEILING(coreocc(i,n)/2)
                  dn = FLOOR(coreocc(i,n)/2)
            END SELECT
            xmlCoreOccs(1,xmlCoreStateNumber,n) = up
            xmlCoreOccs(2,xmlCoreStateNumber,n) = dn
468
            electronsOnAtom = electronsOnAtom + up + dn
469 470 471 472 473
          ENDDO
          DO i = ncorest+1, natomst
            WRITE(6,'("  valence :",2i3,f6.1,i4,a1)') 
     &             coreqn(1,i,n),coreqn(2,i,n),coreocc(i,n),
     &                      coreqn(1,i,n),lotype(lval(i,n))
474
            nel = nel + coreocc(i,n) * atoms%neq(n)
475
            electronsOnAtom = electronsOnAtom + coreocc(i,n)
476

477 478 479 480
c           In d and f shells a magnetic alignment of the spins
c           is preferred in the valence bands
c           Hence the up and down occupation is chosen such that
c           the total spin is maximized
481

482 483 484 485 486 487 488 489
            IF ( abs(coreqn(2,i,n)+0.5) > 2.499 )
     +      THEN
              IF ( coreocc(i,n) > abs(coreqn(2,i,n)) ) THEN
                up = abs(coreqn(2,i,n))
                dn = coreocc(i,n) - abs(coreqn(2,i,n))
              ELSE
                up = coreocc(i,n)
                dn = 0
490 491 492
             END IF
             upreal=up
             dnreal=dn
493

494
c           in s and p states equal occupation of up and down states
495

496 497 498 499 500 501 502
            ELSE
              j = INT(coreocc(i,n) / 2)
              IF (coreocc(i,n) > 2*j) THEN
                j = - coreocc(i,n)
              ENDIF
              up = j
              dn = j
503 504
              upReal = coreocc(i,n) / 2.0
              dnReal = coreocc(i,n) / 2.0
505
            END IF
506 507 508 509
            IF(.NOT.input%l_inpXML) THEN
               WRITE(27,'(4i3,i4,a1)') coreqn(1,i,n),coreqn(2,i,n),
     &                             up,dn,coreqn(1,i,n),lotype(lval(i,n))
            END IF
510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561
            xmlCoreStateNumber = 0
            SELECT CASE(coreqn(1,i,n))
               CASE (1)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 1   !(1s1/2)
               CASE (2)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 2   !(2s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 3    !(2p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 4   !(2p3/2)
               CASE (3)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 5   !(3s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 6    !(3p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 7   !(3p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 9    !(3d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 10  !(3d5/2)
               CASE (4)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 8   !(4s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 11   !(4p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 12  !(4p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 14   !(4d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 15  !(4d5/2)
                  IF(coreqn(2,i,n).EQ.3) xmlCoreStateNumber = 19   !(4f5/2)
                  IF(coreqn(2,i,n).EQ.-4) xmlCoreStateNumber = 20  !(4f7/2)
               CASE (5)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 13  !(5s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 16   !(5p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 17  !(5p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 21   !(5d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 22  !(5d5/2)
                  IF(coreqn(2,i,n).EQ.3) xmlCoreStateNumber = 26   !(5f5/2)
                  IF(coreqn(2,i,n).EQ.-4) xmlCoreStateNumber = 27  !(5f7/2)
               CASE (6)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 18  !(6s1/2)
                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 23   !(6p1/2)
                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 24  !(6p3/2)
                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 28   !(6d3/2)
                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 29  !(6d5/2)
               CASE (7)
                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 25  !(7s1/2)
            END SELECT
            IF(xmlCoreStateNumber.EQ.0) STOP 'Invalid valence state!'
            xmlElectronStates(xmlCoreStateNumber,n) = valenceState_const
            xmlPrintCoreStates(xmlCoreStateNumber,n) = 
     +         coreocc(i,n).NE.xmlCoreRefOccs(xmlCoreStateNumber)
!            SELECT CASE(xmlCoreStateNumber)
!               CASE (9:10,14:15,19:22,26:29)
!                  up = MIN((xmlCoreRefOccs(xmlCoreStateNumber)/2),
!     +                          coreocc(i,n))
!                  dn = MAX(0.0,coreocc(i,n)-up)
!               CASE DEFAULT
!                  up = CEILING(coreocc(i,n)/2)
!                  dn = FLOOR(coreocc(i,n)/2)
!            END SELECT
562 563
            xmlCoreOccs(1,xmlCoreStateNumber,n) = upReal
            xmlCoreOccs(2,xmlCoreStateNumber,n) = dnReal
564 565
          ENDDO
          WRITE (6,*) '----------'
566 567

5392  FORMAT (' atom type: ',i5,' protons: ',f0.8,' electrons: ',f0.8)
568
          IF (ABS(electronsOnAtom-atoms%zatom(n)).GT.1e-13) THEN
569
             WRITE(*,*) 'Note: atom is charged. Is this Intended?'
570
             WRITE(*,5392) n, atoms%zatom(n), electronsOnAtom
571
             WRITE(6,*) 'Note: atom is charged. Is this Intended?'
572
             WRITE(6,5392) n, atoms%zatom(n), electronsOnAtom
573 574
          END IF

575
          IF(.NOT.input%l_inpXML) CLOSE(27)
576 577

          DO i = natomst,1,-1                    ! determine valence states
578 579
            IF (enpara%el0(lval(i,n),n,1) < -9999.8) THEN ! not processed already
              enpara%el0(lval(i,n),n,:) = REAL(coreqn(1,i,n))
580
              IF (i <= ncorest) THEN
581
                enpara%el0(lval(i,n),n,:) = coreqn(1,i,n) + 1.0 ! was already in the core
582 583 584 585
              ENDIF
            ENDIF
          ENDDO
          DO j = 0,3
586 587
            IF (enpara%el0(j,n,1) < -9999.8) THEN
              enpara%el0(j,n,:) = REAL(j+1)
588 589 590 591 592
            ENDIF
          ENDDO

        ELSE  ! determine defauts  as usual

593 594 595
          z = NINT(atoms%zatom(n))
          nlo0 = atoms%nlo(n)
          llo0 = atoms%llo(:,n)
596
          CALL atom_defaults(
597 598
     >                       n,atoms%ntype,atoms%nlod,z,atoms%neq,
     X                       ncst0,nel,atoms%nlo,atoms%llo)
599

600
          IF (atoms%ncst(n) == 0) atoms%ncst(n) = ncst0
601
          IF (lonqn(1,n) /= 0) THEN ! already set before
602 603 604 605
            DO i = 1,atoms%nlo(n)                       ! subtract lo-charge
              nel = nel - 2*(2*atoms%llo(i,n)+1)*atoms%neq(n)
              IF (atoms%llo(i,n) == 0) atoms%ncst(n) = atoms%ncst(n) + 1
              IF (atoms%llo(i,n) >  0) atoms%ncst(n) = atoms%ncst(n) + 2
606
            ENDDO
607 608 609 610 611 612
            atoms%nlo(n) = nlo0                         ! set old values
            atoms%llo(:,n) = llo0 
            DO i = 1,atoms%nlo(n)                       ! add old lo-charge
              nel = nel + 2*(2*atoms%llo(i,n)+1)*atoms%neq(n)   
              IF (atoms%llo(i,n) == 0) atoms%ncst(n) = atoms%ncst(n) - 1
              IF (atoms%llo(i,n) >  0) atoms%ncst(n) = atoms%ncst(n) - 2
613 614
            ENDDO
          ELSE
615 616
             lonqn(1:atoms%nlo(n),n) = 0 !LO check below should not be needed
                                         !for default setting of enparas
617 618
          ENDIF

619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644
          IF(juDFT_was_argument("-electronConfig")) THEN
             electronsLeft = NINT(atoms%zatom(n))
             DO i = 1, 29
                electronsLeft = electronsLeft - xmlCoreRefOccs(i)
                IF(electronsLeft.GT.-xmlCoreRefOccs(i)+eps) THEN
                   natomst = i
                   IF(electronsLeft.LT.-eps) THEN
                      xmlPrintCoreStates(i,n) = .TRUE.
                      SELECT CASE(i)
                         CASE (9:10,14:15,19:22,26:29)
                            up = MIN((xmlCoreRefOccs(i)/2),
     +                               -electronsLeft)
                            dn = MAX(0.0,(-electronsLeft)-up)
                         CASE DEFAULT
                            up = CEILING((-electronsLeft)/2)
                            dn = FLOOR((-electronsLeft)/2)
                      END SELECT
                      xmlCoreOccs(1,i,n) = up
                      xmlCoreOccs(2,i,n) = dn
                   END IF
                END IF
             END DO
             xmlElectronStates(1:atoms%ncst(n),n) = coreState_const
             xmlElectronStates(atoms%ncst(n)+1:natomst,n) = 
     +          valenceState_const
          END IF
645 646 647

        ENDIF

648 649 650
        IF (atoms%nlo(n) /= 0) THEN                    ! check for local orbitals
          DO i = 1, atoms%nlo(n)
            enpara%ello0(i,n,:) = REAL(lonqn(i,n))
651 652 653
            IF (lonqn(i,n) == enpara%qn_el(atoms%llo(i,n),n,1)) THEN  ! increase qn
              enpara%qn_el(atoms%llo(i,n),n,:) = 
     &           enpara%qn_el(atoms%llo(i,n),n,1) + 1          ! in LAPW's by 1
654 655 656
            ENDIF
          ENDDO
        ENDIF
657 658
        enpara%skiplo(n,:) = 0
        DO i = 1, atoms%nlo(n)
659
          enpara%qn_ello(i,n,:) = enpara%qn_el(atoms%llo(i,n),n,:) - 1
660
          enpara%skiplo(n,:) = enpara%skiplo(n,1) + (2*atoms%llo(i,n)+1)
661 662
        ENDDO

663 664 665 666 667 668 669 670 671
      ENDDO

      DO n = 1, atoms%ntype
! correct valence charge
         DO i = 1,atoms%nlo(n)
            IF (atoms%llo(i,n).GT.3) THEN
               nel = nel - 2*(2*atoms%llo(i,n)+1)*atoms%neq(n)   
               IF (atoms%llo(i,n) == 0) atoms%ncst(n) = atoms%ncst(n)+1
               IF (atoms%llo(i,n) >  0) atoms%ncst(n) = atoms%ncst(n)+2
672 673
            ELSE IF (enpara%qn_ello(i,n,1).GE.
     &               enpara%qn_el(atoms%llo(i,n),n,1)) THEN
674 675 676
               nel = nel - 2*(2*atoms%llo(i,n)+1)*atoms%neq(n)   
               IF (atoms%llo(i,n) == 0) atoms%ncst(n) = atoms%ncst(n)+1
               IF (atoms%llo(i,n) >  0) atoms%ncst(n) = atoms%ncst(n)+2
677 678 679 680
            END IF
         ENDDO
      ENDDO

681 682 683
      WRITE (6,'("Valence Electrons =",i5)') nel


684 685 686 687
      enpara%llochg = .FALSE.
      enpara%lchange = .FALSE.
      enpara%enmix = 1.0
      enpara%lchg_v = .TRUE.
688
      IF(juDFT_was_argument("-genEnpara")) THEN
689
         CALL enpara%write(atoms,input%jspins,input%film)
690
      END IF
691
      atoms%lmaxd = lmaxdTemp
692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751
      RETURN

!===> error handling

 911  CONTINUE
      WRITE (errfh,*) 'atom_input: ERROR reading input. ios  =',ios,
     &               ', line =',nline
      CALL juDFT_error("atom_input: ERROR reading input",
     &                calledby="atom_input")

 912  CONTINUE
      WRITE (errfh,*) 'atom_input: ERROR reading namelist.',
     &               ' ios =',ios,
     &               ' line =',nline
      WRITE (errfh,*) buffer(1:nbuffer)
      WRITE (errfh,*) 'The cause of this error may be ...'
      WRITE (errfh,*) '        a variable not defined in this namelist,'
      WRITE (errfh,*) '        wrong type of data for a variable.'
      CALL juDFT_error("atom_input: ERROR reading input",
     &                calledby="atom_input")

 913  CONTINUE
      WRITE (errfh,*) 'atom_input: ERROR reading record.',
     &               ' ios =',ios,
     &               ' line =',nline
      WRITE (errfh,*) buffer(1:nbuffer)
      CALL juDFT_error("atom_input: ERROR reading input",
     &                calledby="atom_input")

!----------------------------------------------------------------
      CONTAINS   ! INTERNAL subroutines
!----------------------------------------------------------------
      SUBROUTINE err( n )

      INTEGER, INTENT (IN) :: n

      WRITE(errfh,*)
      IF (n==1) THEN
        WRITE (errfh,*) 'atom_input: ERROR multiple namelists.',
     &               ' line =',nline
      ELSEIF (n==2) THEN
        WRITE (errfh,*) 'atom_input: ERROR unknown namelist.',
     &               ' line =',nline
      ELSEIF (n==3) THEN
        WRITE (errfh,*) 'atom_input: ERROR reading namelist.',
     &               ' line =',nline
      ELSE
        WRITE (errfh,*) 'atom_input: ERROR reading input.',
     &               ' line =',nline
      ENDIF
      WRITE (errfh,*) buffer(1:nbuffer)
      WRITE (errfh,*)
      fatalerror = .true.
      RETURN
      END SUBROUTINE err
!----------------------------------------------------------------
      END SUBROUTINE atom_input
!----------------------------------------------------------------
!================================================================
      SUBROUTINE read_allatoms(
752
     >                         bfh,
753 754 755 756 757 758 759 760
     X                         rmt,dx,jri,lmax,lnonsph,ncst,econfig,
     <                         bmu,ios)
!****************************************************************
!     reads in defaults for muffin-tin radius, mesh, etc.
!****************************************************************

      IMPLICIT NONE

761
      INTEGER, INTENT (IN)    :: bfh
762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779
      INTEGER, INTENT (INOUT) :: jri     ! mt radial mesh points
      INTEGER, INTENT (INOUT) :: lmax    ! max. l to include for density, overlap etc.
      INTEGER, INTENT (INOUT) :: lnonsph ! max. l for nonspherical MT-contributions
      INTEGER, INTENT (INOUT) :: ncst    ! # of core levels
      INTEGER, INTENT (OUT)   :: ios
      REAL, INTENT (INOUT)    :: rmt, dx ! muffin-tin radius and log. spacing
      REAL, INTENT (INOUT)    :: bmu     ! magnetic moment
      CHARACTER(len=l_buffer) :: econfig ! verbose electronic config

      NAMELIST /allatoms/ rmt,dx,jri,lmax,lnonsph,ncst,econfig,
     &                    bmu

      READ (bfh,allatoms,err=911,end=911,iostat=ios)

 911  CONTINUE
      END SUBROUTINE read_allatoms
!================================================================
      SUBROUTINE read_atom(
780
     >                     bfh,lotype,
781
     X                     id,z,rmt,jri,dx,lmax,lnonsph,ncst,econfig,
782
     <                     speciesName,bmu,lo,nlod,llod,ios )
783 784 785 786 787 788 789 790
!***********************************************************************
!     reads in muffin-tin radius, mesh, etc.
!***********************************************************************

      USE m_element, ONLY : z_namat
      IMPLICIT NONE

! ... arguments ...
791
      INTEGER, INTENT (IN)           :: bfh
792 793 794 795 796 797 798
      REAL, INTENT (OUT)             :: id,z,rmt,dx,bmu
      INTEGER                        :: lmax,lnonsph,ncst,jri,nlod,llod
      CHARACTER(len=l_buffer)        :: econfig
      CHARACTER(len=80)              :: lo
      CHARACTER(len=20), INTENT(OUT) :: speciesName
      INTEGER, INTENT (OUT)          :: ios
      CHARACTER(len=1), INTENT (IN)  :: lotype(0:3)
799 800 801 802 803

! ... internal variables ...
      INTEGER                  :: i,j,k,l,n
      REAL                     :: zz
      CHARACTER(len=2)         :: element
804
      CHARACTER(len=20)        :: name
805 806 807 808 809 810
      CHARACTER(len=80)        :: lo1

      CHARACTER(len=2) :: lotype2(0:3)
      DATA lotype2 /'sS','pP','dD','fF'/

      NAMELIST /atom/ id,z,rmt,dx,jri,lmax,lnonsph,ncst,
811
     &                econfig,bmu,lo,element,name
812 813 814 815

      id = -9999.9
      z  = -9999.9
      element = ' '
816 817
      speciesName = ''
      name = ''
818 819 820

      READ (bfh,atom,err=911,end=911,iostat=ios)

821 822
      speciesName = TRIM(ADJUSTL(name))

823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894
! -> determine which atom we are concerned with ...

      IF ((z < -9999.8).AND.(element.EQ.' ')) THEN
        WRITE (errfh,*)
     &       'ERROR! No element specified  in namelist atom...'
        WRITE (errfh,*) 'use z=.. or element=.. to define it!'
        ios = 3001
        RETURN
      ENDIF
      IF (id < -9999.8) THEN     ! if no id specified
        zz = REAL(z_namat(element))
        IF (z < 0.00) THEN
          IF (zz > eps) THEN
            id = zz              ! use element name
          ELSE
            id = z               ! or use "z" for id
          ENDIF
        ELSE
          IF (zz > 0 .AND. abs(zz-z)>0.5) THEN
            WRITE (warnfh,*)
            WRITE (warnfh,*) 'atom_input: WARNING! ',
     &       'z and z of specified element differ by more than 0.5. '
            WRITE (warnfh,*) '  z = ', z, 'element =',element
            WRITE (warnfh,*)
          ENDIF
        ENDIF
      ENDIF

!---> order local orbitals determine nlod, llod
      lo1 = adjustl(lo)
      lo = ' '
      i = 0
      n = 0
      DO l = 0, 3
        DO
          j = SCAN(lo1,lotype2(l))  ! search for 's' or 'S', 'p' or 'P' etc.
          IF (j > 0) THEN
            lo1(j:j) = ' '
            n = n + 1
            i = i + 1
            IF (j > 1) THEN
              k = SCAN(lo1(j-1:j-1),'123456') ! determine principal quantum number
              IF (k > 0) THEN
                lo(i:i) = lo1(j-1:j-1)
                lo1(j-1:j-1) = ' '
              ELSE
                lo(i:i) = '0'
              ENDIf
            ELSE
              lo(i:i) = '0'
            ENDIF
            i = i + 1
            lo(i:i) = lotype(l)
            nlod = max( nlod, n )
            llod = max( llod, l )
          ELSE
            EXIT
          ENDIF
        ENDDO
      ENDDO
      IF (len_trim(lo1) > 0) then
        WRITE (errfh,*) 'ERROR reading local orbital input...',lo1
        ios = 3002
      ENDIF

 911  CONTINUE
      END SUBROUTINE read_atom
!================================================================

      SUBROUTINE atom_defaults(
     >                         n,ntype,nlod,z,neq,
     X                         ncst2,nel,nlo,llo)
895
      USE m_juDFT
896 897 898 899 900 901 902 903
      IMPLICIT NONE

      INTEGER, INTENT (IN)    :: n,ntype,nlod,z
      INTEGER, INTENT (IN)    :: neq(ntype)
      INTEGER, INTENT (INOUT) :: nel,ncst2
      INTEGER, INTENT (INOUT) :: nlo(ntype),llo(nlod,ntype)
      

904 905 906 907 908 909 910
      INTEGER locore,lo
      INTEGER ncst1(0:103),nce(0:24)

!
! electrons associated with a given number of core-levels
!
      nce=-1;
911 912
      nce(0) = 0  ; nce(1) = 2  ; nce(2)= 4   ; nce(4) = 10  
      nce(5) = 12 ; nce(7) =18 ;nce(8)=20; 
913 914 915 916
      nce(9) = 28 ; nce(12) = 36; nce(14) = 46; nce(17) = 54
      nce(19) = 68; nce(21) = 78; nce(24) = 86


917 918
!
! number of core levels for each element; the INT(ncst1/100) number
919 920 921 922 923 924 925
! provides information about possible local orbitals: 
!    0  no LO
!    1  s-LO
!    2  p-LO
!    4  d-LO
!    8  f-LO
! Sums are allowed, i.e. 3 (s,p)-LO
926
!
927 928 929

!Defaults
      ncst1 =(/0,0,                                            0,        ! Va,H,He
930 931 932
     + 01, 01,                                         1, 1, 1, 1, 1, 1,        ! Li - Ne
     + 304,304,                                        4, 4, 4, 4, 4, 4,       ! Na - Ar
     + 307,307,307,307,307,307,307,307,207,207, 7,  7,409,409,409,409,
933
     +                                                          409, 9,  ! K - Kr
934
     + 312,312,312,312,312,312,312,212,212,212,312,212,414,414,414,414,
935 936
     +                                                         414,414,  ! Rb - Xe
     + 317,317,217,217,217,217,217,217,217,217,217,17, 17,17,17,17,17,    ! Cs - Lu
937
     +    1119,1119,319,319,219,219,219,219,219,421,421,421,421,421,421,  ! Hf - Rn
938 939
     + 324,324,224,224,224,24,24,24,24,24,24,24, 24,24,24,24,24/)   ! Fr - Lw

940
      if (judft_was_argument("-fast_defaults")) 
941
     + ncst1 =(/0,0,                                                0,  ! Va,H,He
942 943
     +     01, 01,                                  1, 1, 1, 1, 1, 1,  ! Li - Ne
     +     04, 04,                                  4, 4, 4, 4, 4, 4,  ! Na - Ar
944 945 946 947 948 949
     +    307,307,207,207, 7, 7, 7, 7, 7, 7, 7, 7,409, 9, 9, 9, 9, 9,  ! K - Kr
     +    312,312,212,212,12,12,12,12,12,12,12,12,414,14,14,14,14,14,  ! Rb - Xe
     +    317,317,217,217,17,17,17,17,17,17,17,17, 17,17,17,17,17,     ! Cs - Lu
     +                219,19,19,19,19,19,19,19,19,421,21,21,21,21,21,  ! Hf - Rn
     +    324,324,224,224,224,24,24,24,24,24,24,24, 24,24,24,24,24/)   ! Fr - Lw

950 951 952 953

!
!--> determine core levels
!
954 955 956 957 958 959 960 961 962 963
      ncst2 = mod(ncst1( z ),100)
      lo=int(ncst1(z)/100)
      
      nel=nel+(z - nce(ncst2))*neq(n)
      nlo(n) = 0 ; locore = 0

      IF (btest(lo,0)) THEN !s-lo
         locore=locore+2
         ncst2=ncst2-1
         nlo(n)=nlo(n)+1
964
         llo(nlo(n),n)=0
965
      ENDIF
966 967 968 969
      IF (btest(lo,1)) THEN !p-lo
         locore=locore+6
         ncst2=ncst2-2
         nlo(n)=nlo(n)+1
970
         llo(nlo(n),n)=1
971 972 973 974 975
      ENDIF
      IF (btest(lo,2)) THEN !d-lo
         locore=locore+10
         ncst2=ncst2-2
         nlo(n)=nlo(n)+1
976
         llo(nlo(n),n)=2
977 978 979 980 981
      ENDIF
      IF (btest(lo,3)) THEN !f-lo
         locore=locore+14
         ncst2=ncst2-2
         nlo(n)=nlo(n)+1
982
         llo(nlo(n),n)=3
983 984 985 986
      ENDIF
      
      nel = nel +  locore  * neq(n)
      
987 988 989
      END SUBROUTINE atom_defaults

      END MODULE m_atominput