r_inpXML.F90 91.6 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
MODULE m_rinpXML
8 9 10 11 12 13 14
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!
   !!! The routine r_inpXML reads in the inp.xml file
   !!!
   !!!                               GM'16
   !!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
15
CONTAINS
16
   SUBROUTINE r_inpXML(&
Daniel Wortmann's avatar
Daniel Wortmann committed
17
      atoms,vacuum,input,stars,sliceplot,banddos,DIMENSION,forcetheo,field,&
18
      cell,sym,xcpot,noco,oneD,hybrid,kpts,enpara,coreSpecInput,wann,&
19
      noel,namex,relcor,a1,a2,a3,dtild)
20 21 22 23 24 25

      USE iso_c_binding
      USE m_juDFT
      USE m_types
      USE m_types_forcetheo_extended
      USE m_symdata , ONLY : nammap, ord2, l_c2
26
      !USE m_rwsymfile
27
      !USE m_xmlIntWrapFort
28
      USE m_inv3
29
      !USE m_spg2set
30
      !USE m_closure, ONLY : check_close
31
      !USE m_symproperties
32 33
      USE m_calculator
      USE m_constants
34
      !USE m_inpeig
35 36
      USE m_sort
      USE m_types_xcpot_inbuild
37
#ifdef CPP_LIBXC
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
      USE xc_f03_lib_m
#endif
      IMPLICIT NONE

      TYPE(t_input),INTENT(INOUT)   :: input
      TYPE(t_sym),INTENT(INOUT)     :: sym
      TYPE(t_stars),INTENT(INOUT)   :: stars
      TYPE(t_atoms),INTENT(INOUT)   :: atoms
      TYPE(t_vacuum),INTENT(INOUT)   :: vacuum
      TYPE(t_kpts),INTENT(INOUT)     :: kpts
      TYPE(t_oneD),INTENT(INOUT)     :: oneD
      TYPE(t_hybrid),INTENT(INOUT)   :: hybrid
      TYPE(t_cell),INTENT(INOUT)     :: cell
      TYPE(t_banddos),INTENT(INOUT)  :: banddos
      TYPE(t_sliceplot),INTENT(INOUT):: sliceplot
      CLASS(t_xcpot),INTENT(INOUT),ALLOCATABLE :: xcpot
      TYPE(t_noco),INTENT(INOUT)     :: noco
55
      TYPE(t_dimension),INTENT(OUT)  :: dimension
56
      TYPE(t_enpara)   ,INTENT(OUT)  :: enpara
57
      TYPE(t_field), INTENT(INOUT)   :: field
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90
      CLASS(t_forcetheo),ALLOCATABLE,INTENT(OUT):: forcetheo
      TYPE(t_coreSpecInput),INTENT(OUT) :: coreSpecInput
      TYPE(t_wann)   ,INTENT(INOUT)  :: wann
      CHARACTER(len=3), ALLOCATABLE, INTENT(INOUT) :: noel(:)
      CHARACTER(len=4), INTENT(OUT)  :: namex
      CHARACTER(len=12), INTENT(OUT) :: relcor
      REAL, INTENT(OUT)              :: a1(3),a2(3),a3(3)
      REAL, INTENT(OUT)              :: dtild

      CHARACTER(len=8) :: name(10)

      !+odim
      INTEGER MM,vM,m_cyl
      LOGICAL invs1,zrfs1
      INTEGER chi,rot
      LOGICAL d1,band
      NAMELIST /odim/ d1,MM,vM,m_cyl,chi,rot,invs1,zrfs1
      !-odim
      ! ..
      ! ..  Local Variables
      REAL     :: scpos  ,zc
      INTEGER ieq,i,k,na,n,ii,vxc_id_c,vxc_id_x,exc_id_c,exc_id_x
      REAL s3,ah,a,hs2,rest,thetaj
      LOGICAL l_hyb,l_sym,ldum
      INTEGER :: ierr
      ! ..
      !...  Local Arrays
      !   CHARACTER :: helpchar(atoms%ntype)
      CHARACTER(len=  4) :: chntype
      CHARACTER(len= 41) :: chform
      CHARACTER(len=100) :: line

      CHARACTER(len=20) :: tempNumberString
91
      CHARACTER(len=150) :: format
92 93 94
      CHARACTER(len=20) :: mixingScheme
      CHARACTER(len=10) :: loType
      LOGICAL :: kptGamma, l_relcor,ldummy
95
      INTEGER :: iAtomType
96
      CHARACTER(len=100) :: xPosString, yPosString, zPosString
Daniel Wortmann's avatar
Daniel Wortmann committed
97 98
      CHARACTER(len=20),ALLOCATABLE :: speciesName(:)
      
99
      !   REAL :: tempTaual(3,atoms%nat)
100
      TYPE(t_econfig)  :: econf
101 102 103

      INTEGER            :: iType, iLO, iSpecies, lNumCount, nNumCount, iLLO, jsp, j, l, absSum, numTokens
      INTEGER            :: numberNodes, nodeSum, numSpecies, n2spg, n1, n2, ikpt, iqpt
104
      INTEGER            :: atomicNumber,  gridPoints, lmax, lnonsphr, lmaxAPW
105
      INTEGER            :: latticeDef, symmetryDef, nop48, firstAtomOfType, errorStatus
106
      INTEGER            :: loEDeriv, ntp1, ios, ntst, jrc, minNeigd
107 108 109 110 111 112 113
      INTEGER            :: nv, nv2, kq1, kq2, kq3, nprncTemp, kappaTemp, tempInt
      INTEGER            :: ldau_l(4), numVac, numU
      INTEGER            :: speciesEParams(0:3)
      INTEGER            :: mrotTemp(3,3,48)
      REAL               :: tauTemp(3,48)
      REAL               :: bk(3)
      LOGICAL            :: flipSpin, l_eV, invSym, l_qfix, relaxX, relaxY, relaxZ
Daniel Wortmann's avatar
Daniel Wortmann committed
114
      LOGICAL            :: coreConfigPresent, l_enpara, l_orbcomp, tempBool, l_nocoinp
115 116
      REAL               :: magMom, radius, logIncrement, qsc(3), latticeScale, dr
      REAL               :: aTemp, zp, rmtmax, sumWeight, ldau_u(4), ldau_j(4), tempReal
117
      REAL               :: ldau_phi(4),ldau_theta(4)
118 119 120 121
      REAL               :: weightScale, eParamUp, eParamDown
      LOGICAL            :: l_amf(4)
      REAL, PARAMETER    :: boltzmannConst = 3.1668114e-6 ! value is given in Hartree/Kelvin
      INTEGER            :: lcutm,lcutwf,hybSelect(4)
122
      REAL               :: evac0Temp(2,2)
123 124 125 126 127 128

      CHARACTER(LEN=200,KIND=c_char) :: schemaFilename, docFilename
      CHARACTER(LEN=255) :: valueString, lString, nString, token
      CHARACTER(LEN=255) :: xPathA, xPathB, xPathC, xPathD, xPathE
      CHARACTER(LEN=11)  :: latticeType
      CHARACTER(LEN=50)  :: versionString
129
      CHARACTER(LEN=150) :: kPointsPrefix
130

131
      INTEGER            :: altKPointSetIndex,  altKPointSetIndices(2)
132
      LOGICAL            :: ldaSpecies
Daniel Wortmann's avatar
Daniel Wortmann committed
133
      REAL               :: socscaleSpecies,b_field_mtspecies,vcaspecies
134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152

      INTEGER, ALLOCATABLE :: lNumbers(:), nNumbers(:), speciesLLO(:)
      INTEGER, ALLOCATABLE :: loOrderList(:)
      INTEGER, ALLOCATABLE :: speciesNLO(:)
      INTEGER, ALLOCATABLE :: multtab(:,:), invOps(:), optype(:)
      INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:)
      INTEGER, ALLOCATABLE :: speciesLOEDeriv(:)
      REAL,    ALLOCATABLE :: speciesLOeParams(:), speciesLLOReal(:)
      LOGICAL, ALLOCATABLE :: wannAtomList(:)

      ! Variables for MT radius testing:

      REAL                 :: dtild1,kmax1,dvac1
      LOGICAL              :: l_test
      INTEGER, ALLOCATABLE :: jri1(:), lmax1(:)
      REAL, ALLOCATABLE    :: rmt1(:), dx1(:)

      EXTERNAL prp_xcfft_box

153 154 155
      INTERFACE
         FUNCTION dropInputSchema() BIND(C, name="dropInputSchema")
            USE iso_c_binding
156
            INTEGER(c_int) dropInputSchema
157 158
         END FUNCTION dropInputSchema
      END INTERFACE
159 160 161 162

      errorStatus = 0
      errorStatus = dropInputSchema()
      IF(errorStatus.NE.0) THEN
163
         CALL juDFT_error('Error: Cannot print out FleurInputSchema.xsd')
164 165 166 167 168 169
      END IF

      schemaFilename = "FleurInputSchema.xsd"//C_NULL_CHAR
      docFilename = "inp.xml"//C_NULL_CHAR

      !TODO! these switches should be in the inp-file
170 171
      input%l_core_confpot=.TRUE. !former CPP_CORE
      input%l_useapw=.FALSE.   !former CPP_APW
172 173 174 175 176 177 178 179
      !WRITE(*,*) 'Start reading of inp.xml file'
      CALL xmlInitInterface()
      CALL xmlParseSchema(schemaFilename)
      CALL xmlParseDoc(docFilename)
      CALL xmlValidateDoc()
      CALL xmlInitXPath()

      ! Check version of inp.xml
180 181
      versionString = xml%GetAttributeValue('/fleurInput/@fleurInputVersion')
      IF((TRIM(ADJUSTL(versionString)).NE.'0.30')) THEN
182
         CALL juDFT_error('version number of inp.xml file is not compatible with this fleur version')
183 184
      END IF

185 186 187 188 189 190 191

      !Some types have their own readers
      CALL kpts%read_xml(xml)
      CALL sym%read_sym(xml)

      

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
      ! Get number of atoms, atom types, and atom species

      numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/relPos')
      numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/absPos')
      numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/filmPos')

      atoms%nat = numberNodes

      numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup')

      atoms%ntype = numberNodes

      numSpecies = xmlGetNumberOfNodes('/fleurInput/atomSpecies/species')

      ALLOCATE(atoms%nz(atoms%ntype))     !nz and zatom have the same content!
      ALLOCATE(atoms%zatom(atoms%ntype))  !nz and zatom have the same content!
      ALLOCATE(atoms%jri(atoms%ntype))
      ALLOCATE(atoms%dx(atoms%ntype))
      ALLOCATE(atoms%lmax(atoms%ntype))
      ALLOCATE(atoms%nlo(atoms%ntype))
212
      ALLOCATE(atoms%econf(atoms%ntype))
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
      ALLOCATE(atoms%lnonsph(atoms%ntype))
      ALLOCATE(atoms%nflip(atoms%ntype))
      ALLOCATE(atoms%l_geo(atoms%ntype))
      ALLOCATE(atoms%lda_u(4*atoms%ntype))
      ALLOCATE(atoms%bmu(atoms%ntype))
      ALLOCATE(atoms%relax(3,atoms%ntype))
      ALLOCATE(atoms%neq(atoms%ntype))
      ALLOCATE(atoms%taual(3,atoms%nat))
      ALLOCATE(atoms%label(atoms%nat))
      ALLOCATE(atoms%pos(3,atoms%nat))
      ALLOCATE(atoms%rmt(atoms%ntype))
      ALLOCATE(atoms%namex(atoms%ntype))
      ALLOCATE(atoms%icorr(atoms%ntype))
      ALLOCATE(atoms%igrd(atoms%ntype))
      ALLOCATE(atoms%krla(atoms%ntype))
      ALLOCATE(atoms%relcor(atoms%ntype))

      atoms%namex = ''
      atoms%icorr = -99

      ALLOCATE(atoms%ncv(atoms%ntype)) ! For what is this?
      ALLOCATE(atoms%ngopr(atoms%nat)) ! For what is this?
      ALLOCATE(atoms%lapw_l(atoms%ntype)) ! Where do I put this?
      ALLOCATE(atoms%invsat(atoms%nat)) ! Where do I put this?

      ALLOCATE(noco%l_relax(atoms%ntype),noco%b_con(2,atoms%ntype))
      ALLOCATE(noco%alphInit(atoms%ntype),noco%alph(atoms%ntype),noco%beta(atoms%ntype))
      ALLOCATE(noco%socscale(atoms%ntype))

Daniel Wortmann's avatar
Daniel Wortmann committed
242
   
243
    
244 245 246 247 248 249 250 251 252 253 254 255 256 257

      ALLOCATE (wannAtomList(atoms%nat))

      ! Read in constants

      xPathA = '/fleurInput/constants/constant'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      DO i = 1, numberNodes
         WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)), '[',i,']'
         tempReal = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@value'))
         valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@name')
         CALL ASSIGN_var(valueString,tempReal)
      END DO

258 259 260
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Comment section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
261 262 263 264 265

      input%comment = '        '
      xPathA = '/fleurInput/comment'
      valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
      DO i = 1, LEN(TRIM(ADJUSTL(valueString)))
Matthias Redies's avatar
Matthias Redies committed
266
         IF (valueString(i:i) == ACHAR(10)) valueString(i:i) = ' ' !remove line breaks
267
      END DO
Daniel Wortmann's avatar
Daniel Wortmann committed
268 269
      input%comment = TRIM(ADJUSTL(valueString))
    
270 271 272
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Start of calculationSetup section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
273 274 275 276 277 278 279 280 281 282

      ! Read general cutoff parameters

      input%rkmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Kmax'))
      stars%gmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Gmax'))

      stars%gmaxInit = stars%gmax

      xPathA = '/fleurInput/calculationSetup/cutoffs/@numbands'
      numberNodes = xmlGetNumberOfNodes(xPathA)
283
      DIMENSION%neigd = 0
284 285 286
      IF(numberNodes.EQ.1) THEN
         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
         IF(TRIM(ADJUSTL(valueString)).EQ.'all') THEN
287
            dimension%neigd = -1
288
         ELSE
289
            READ(valueString,*) DIMENSION%neigd
290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308
         END IF
      END IF

      ! Read SCF loop parametrization

      input%itmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@itmax'))
      input%minDistance = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@minDistance'))
      input%maxiter = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@maxIterBroyd'))

      valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@imix')))
      SELECT CASE (valueString)
      CASE ('straight')
         input%imix = 0
      CASE ('Broyden1')
         input%imix = 3
      CASE ('Broyden2')
         input%imix = 5
      CASE ('Anderson')
         input%imix = 7
Daniel Wortmann's avatar
Daniel Wortmann committed
309 310
      CASE ("Pulay")
         input%imix = 9
311 312 313 314 315 316
      CASE ("pPulay")
         input%imix = 11
      CASE ("rPulay")
         input%imix = 13
      CASE ("aPulay")
         input%imix = 15
317
      CASE DEFAULT
318
         CALL juDFT_error('Error: unknown mixing scheme selected!')
319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343
      END SELECT

      input%alpha = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@alpha'))
input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@preconditioning_param'))
      input%spinf = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@spinf'))

      ! Get parameters for core electrons

      input%ctail = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@ctail'))
      IF((TRIM(ADJUSTL(versionString)).EQ.'0.27')) THEN
         input%coretail_lmax = 99
      ELSE
       input%coretail_lmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@coretail_lmax'))
      END IF
      input%frcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@frcor'))
      input%kcrel = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@kcrel'))

      ! Read in magnetism parameters

      input%jspins = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@jspins'))
      noco%l_noco = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@l_noco'))
      input%swsp = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@swsp'))
      input%lflip = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@lflip'))
      input%fixed_moment=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@fixed_moment'))

344
  IF (ABS(input%fixed_moment)>1E-8.AND.(input%jspins==1.OR.noco%l_noco)) CALL judft_error("Fixed moment only in collinear calculations with two spins")
345

346 347 348 349 350 351 352 353 354 355 356 357 358
      ! Read in optional expert modes switches

      xPathA = '/fleurInput/calculationSetup/expertModes'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      input%gw = 0
      input%secvar = .FALSE.

      IF (numberNodes.EQ.1) THEN
         input%gw = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gw'))
         input%secvar = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@secvar'))
      END IF

359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391
      ! Check for alternative k point sets for the chosen FLEUR mode

      xPathA = '/fleurInput/output'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      banddos%band = .FALSE.
      IF (numberNodes.EQ.1) THEN
         banddos%band = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@band'))
      END IF

      altKPointSetIndices(:) = -1
      xPathA = '/fleurInput/calculationSetup/bzIntegration/altKPointSet'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      IF(numberNodes.NE.0) THEN
         DO i = 1, numberNodes
            WRITE(xPathA,*) '/fleurInput/calculationSetup/bzIntegration/altKPointSet[',i,']/@purpose'
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
            IF((altKPointSetIndices(2).EQ.-1).AND.(TRIM(ADJUSTL(valueString)).EQ.'GW')) THEN
               altKPointSetIndices(2) = i
            ELSE IF((altKPointSetIndices(1).EQ.-1).AND.(TRIM(ADJUSTL(valueString)).EQ.'bands')) THEN
               altKPointSetIndices(1) = i
            END IF
         END DO
      END IF

      altKPointSetIndex = -1
      IF(banddos%band) THEN
         altKPointSetIndex = altKPointSetIndices(1)
      ELSE IF (input%gw.EQ.2) THEN
         altKPointSetIndex = altKPointSetIndices(2)
      END IF

      IF (altKPointSetIndex.NE.-1) THEN
         WRITE(kPointsPrefix,*) '/fleurInput/calculationSetup/bzIntegration/altKPointSet[',altKPointSetIndex,']'
392
      END IF
393 394 395

      ! Read in Brillouin zone integration parameters

396
    
397 398 399 400 401 402 403 404 405 406 407 408
      valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/@mode')))
      SELECT CASE (valueString)
      CASE ('hist')
         input%gauss = .FALSE.
         input%tria = .FALSE.
      CASE ('gauss')
         input%gauss = .TRUE.
         input%tria = .FALSE.
      CASE ('tria')
         input%gauss = .FALSE.
         input%tria = .TRUE.
      CASE DEFAULT
409
         CALL juDFT_error('Invalid bzIntegration mode selected!')
410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426
      END SELECT

      nodeSum = 0
      xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingEnergy'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      nodeSum = nodeSum + numberNodes
      IF (numberNodes.EQ.1) THEN
         input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
      END IF
      xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingTemp'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      nodeSum = nodeSum + numberNodes
      IF (numberNodes.EQ.1) THEN
         input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
         input%tkb = boltzmannConst * input%tkb
      END IF
      IF(nodeSum.GE.2) THEN
427
         CALL juDFT_error('Error: Multiple fermi Smearing parameters provided in input file!')
428 429 430 431 432 433 434
      END IF

      xPathA = '/fleurInput/calculationSetup/bzIntegration/@valenceElectrons'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      IF (numberNodes.EQ.1) THEN
         input%zelec = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
      ELSE
435
         CALL juDFT_error('Error: Optionality of valence electrons in input file not yet implemented!')
436 437
      END IF

438
      IF (altKPointSetIndex.EQ.-1) THEN
439 440 441
         WRITE(kPointsPrefix,*) '/fleurInput/calculationSetup/bzIntegration'
      END IF

442
      call judft_error("BUG reading of kpts")
443
      ! Option kPointDensity
444
 
445
         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(kPointsPrefix))//'/kPointCount/specialPoint')
446
         IF(numberNodes.EQ.1) THEN
447
            CALL juDFT_error('Error: Single special k point provided. This does not make sense!')
448 449 450 451 452 453 454
         END IF
         kpts%numSpecialPoints = numberNodes
         IF(kpts%numSpecialPoints.GE.2) THEN
            DEALLOCATE(kpts%specialPoints)
            ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints))
            ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints))
            DO i = 1, kpts%numSpecialPoints
455
               WRITE(xPathA,*) TRIM(ADJUSTL(kPointsPrefix))//'/kPointCount/specialPoint[',i,']'
456 457 458 459 460 461 462 463 464
               valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
               kpts%specialPoints(1,i) = evaluatefirst(valueString)
               kpts%specialPoints(2,i) = evaluatefirst(valueString)
               kpts%specialPoints(3,i) = evaluatefirst(valueString)
               kpts%specialPointNames(i) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name')
            END DO
         END IF

      ! Option kPointList
465
      numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList')
466 467
      IF (numberNodes.EQ.1) THEN
         l_kpts = .TRUE.
468
         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/kPoint')
469 470 471 472 473 474
         kpts%nkpt = numberNodes
         kpts%l_gamma = .FALSE.
         ALLOCATE(kpts%bk(3,kpts%nkpt))
         ALLOCATE(kpts%wtkpt(kpts%nkpt))
         kpts%bk = 0.0
         kpts%wtkpt = 0.0
475 476
 
          weightScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/@weightScale'))
477 478

         DO i = 1, kpts%nkpt
479
            WRITE(xPathA,*) TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/kPoint[',i,']'
480 481
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
            READ(valueString,*) kpts%bk(1,i), kpts%bk(2,i), kpts%bk(3,i)
482
            kpts%bk(:,i)=kpts%bk(:,i)
483 484 485 486 487
            kpts%wtkpt(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@weight'))
            kpts%wtkpt(i) = kpts%wtkpt(i) / weightScale
         END DO
      END IF

488 489 490 491 492 493 494 495 496 497 498 499 500
      ! Option kPointListFile
      xPathA = TRIM(ADJUSTL(kPointsPrefix))//'/kPointListFile'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      IF (numberNodes.EQ.1) THEN
         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@filename')))
         OPEN (41,file=TRIM(ADJUSTL(valueString)),form='formatted',status='old')
            READ (41,*) kpts%nkpt
         CLOSE (41)
         ALLOCATE(kpts%bk(3,kpts%nkpt))
         ALLOCATE(kpts%wtkpt(kpts%nkpt))
         kpts%bk = 0.0
         kpts%wtkpt = 0.0
         kpts%l_gamma = .FALSE.
501
         l_kpts = .TRUE.
502
     
503 504 505 506 507 508
         ! We need to set input%film for inpeig. Unfortunately this is actually initialized at a later stage.
         ! So we do it here additionally.
         input%film = .FALSE.
         numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/filmLattice')
         IF (numberNodes.EQ.1) input%film = .TRUE.

509
         !CALL inpeig(atoms,cell,input,.FALSE.,kpts,kptsFilename=TRIM(ADJUSTL(valueString)))
510 511
      END IF

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
      ! Read in optional SOC parameters if present

      xPathA = '/fleurInput/calculationSetup/soc'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      noco%l_soc = .FALSE.
      noco%theta = 0.0
      noco%phi = 0.0

      IF (numberNodes.EQ.1) THEN
         noco%theta=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta'))
         noco%phi=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi'))
         noco%l_soc = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_soc'))
         noco%l_spav = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spav'))
      END IF

      ! Read in optional noco parameters if present

      xPathA = '/fleurInput/calculationSetup/nocoParams'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      noco%l_ss = .FALSE.
      noco%l_mperp = .FALSE.
      noco%l_constr = .FALSE.
      noco%mix_b = 0.0
      noco%qss = 0.0

      noco%l_relax(:) = .FALSE.
      noco%alphInit(:) = 0.0
      noco%alph(:) = 0.0
      noco%beta(:) = 0.0
      noco%b_con(:,:) = 0.0

      IF ((noco%l_noco).AND.(numberNodes.EQ.0)) THEN
546
         CALL juDFT_error('Error: l_noco is true but no noco parameters set in xml input file!')
547 548 549 550 551 552
      END IF

      IF (numberNodes.EQ.1) THEN
         noco%l_ss = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_ss'))
         noco%l_mperp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mperp'))
         noco%l_constr = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_constr'))
553
         noco%l_mtNocoPot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mtNocoPot'))
554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601

         noco%mix_b = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mix_b'))

         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qss')))
         READ(valueString,*) noco%qss(1), noco%qss(2), noco%qss(3)

         !WRITE(*,*) 'Note: TODO: Calculation of q points!'

         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/qsc')
         IF (numberNodes.EQ.1) THEN
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qsc')))
            READ(valueString,*) qsc(1), qsc(2), qsc(3)
            DO i = 1, 3
               noco%qss(i) = noco%qss(i) / qsc(i)
            END DO
            !WRITE(*,*) 'Note: TODO: Integrate qsc directly into qss in input file!'
            !WRITE(*,*) '(no problem for users)'
         END IF
      END IF

      ! Read in optional 1D parameters if present

      xPathA = '/fleurInput/calculationSetup/oneDParams'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      oneD%odd%d1 = .FALSE.

      IF (numberNodes.EQ.1) THEN
         oneD%odd%d1 = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@d1'))
         oneD%odd%M = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@MM'))
         oneD%odd%mb = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vM'))
         oneD%odd%m_cyl = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@m_cyl'))
         oneD%odd%chi = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@chi'))
         oneD%odd%rot = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@rot'))
         oneD%odd%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs1'))
         oneD%odd%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs1'))
      END IF

      ! Read in optional geometry optimization parameters

      xPathA = '/fleurInput/calculationSetup/geometryOptimization'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      input%l_f = .FALSE.
      input%qfix = 0

      IF (numberNodes.EQ.1) THEN
         input%l_f = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_f'))
602
         input%forcealpha = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@forcealpha'))
603 604
         input%epsdisp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsdisp'))
         input%epsforce = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsforce'))
605 606
         input%forcemix = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@forcemix'))
         input%force_converged = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@force_converged'))
607
         input%qfix = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@qfix'))
608 609 610 611 612 613 614 615 616 617 618 619
      END IF

      ! Read in optional general LDA+U parameters

      xPathA = '/fleurInput/calculationSetup/ldaU'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      IF (numberNodes.EQ.1) THEN
         input%ldauLinMix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_linMix'))
         input%ldauMixParam = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mixParam'))
         input%ldauSpinf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spinf'))
      END IF

620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643
      ! Read in RDMFT parameters

      input%l_rdmft = .FALSE.
      input%rdmftOccEps = 0.00001
      input%rdmftStatesBelow = 5
      input%rdmftStatesAbove = 5
      input%rdmftFunctional = -1

      xPathA = '/fleurInput/calculationSetup/rdmft'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      IF (numberNodes.EQ.1) THEN
         input%l_rdmft = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_rdmft'))
         input%rdmftOccEps = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@occEps'))
         input%rdmftStatesBelow = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@statesBelow'))
         input%rdmftStatesAbove = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@statesAbove'))
         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@functional')))
         SELECT CASE (valueString)
            CASE ('Muller')
               input%rdmftFunctional = 1
            CASE DEFAULT
               STOP 'Error: unknown RDMFT functional selected!'
         END SELECT
      END IF

644 645 646 647 648 649
      ! Read in optional q point mesh for spin spirals

      xPathA = '/fleurInput/calculationSetup/spinSpiralQPointMesh'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      !   IF ((noco%l_ss).AND.(numberNodes.EQ.0)) THEN
Matthias Redies's avatar
Matthias Redies committed
650
      !      call juDFT_error('Error: l_ss is true but no q point mesh set in xml input file!')
651 652 653
      !   END IF

      ! Read in optional E-Field parameters
654 655
     
      xPathA = '/fleurInput/calculationSetup/fields'
656
      numberNodes = xmlGetNumberOfNodes(xPathA)
Daniel Wortmann's avatar
Daniel Wortmann committed
657 658 659 660
      field%b_field=0.0
      field%l_b_field=.FALSE.
      field%efield%sigma=0.0
      
661 662

      IF (numberNodes.EQ.1) THEN
663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681
         IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@b_field')>0) THEN
            field%b_field=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'//@b_field'))
            field%l_b_field=.true.
         ENDIF
         field%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma'))
         field%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1'))
         field%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2'))
         field%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge'))
         field%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho'))
         field%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp'))
         field%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet'))
         field%efield%l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV'))

         numberNodes=xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/shape')
         ALLOCATE(field%efield%shapes(numberNodes))
         DO i=1,numberNodes
            WRITE(xPathB,"(a,a,i0,a)") TRIM(ADJUSTL(xpathA)),'/shape[',i,']'
            field%efield%shapes(i)=TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
         ENDDO
Daniel Wortmann's avatar
Daniel Wortmann committed
682 683
      ELSE
         ALLOCATE(field%efield%shapes(0))
684 685 686 687 688 689 690 691 692 693 694 695
      END IF

      ! Read in optional energy parameter limits

      xPathA = '/fleurInput/calculationSetup/energyParameterLimits'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      IF (numberNodes.EQ.1) THEN
         input%ellow = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ellow'))
         input%elup = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@elup'))
      END IF

696 697 698
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! End of calculationSetup section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
699

700 701 702
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Start of cell section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
703 704 705 706 707 708 709 710 711 712 713 714

      ! Read in lattice parameters

      a1 = 0.0
      a2 = 0.0
      a3 = 0.0
      cell%z1 = 0.0
      dtild = 0.0
      input%film = .FALSE.
      latticeType = 'bulkLattice'
      latticeDef = 0
      symmetryDef = 0
715
    
716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734
      numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/filmLattice')

      IF (numberNodes.EQ.1) THEN
         input%film = .TRUE.
         latticeType = 'filmLattice'
      END IF

      xPathA = '/fleurInput/cell/'//latticeType
      numberNodes = xmlGetNumberOfNodes(xPathA)

      IF (numberNodes.EQ.1) THEN
         latticeScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@scale'))
         input%scaleCell = latticeScale

         IF(input%film) THEN
            cell%z1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dVac'))
            dtild = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dTilda'))
            vacuum%dvac = cell%z1
            a3(3) = dtild
735
            evac0Temp = eVac0Default_const
736 737 738 739 740 741 742 743 744 745 746 747
            xPathB = TRIM(ADJUSTL(xPathA))//'/vacuumEnergyParameters'
            numberNodes = xmlGetNumberOfNodes(xPathB)
            IF(numberNodes.GE.1) THEN
               DO i = 1, numberNodes
                  xPathC = ''
                  WRITE(xPathC,'(a,i0,a)') TRIM(ADJUSTL(xPathB))//'[',i,']'
                  numVac = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@vacuum'))
                  eParamUp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinUp'))
                  eParamDown = eParamUp
                  IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathC))//'/@spinDown').GE.1) THEN
                     eParamDown = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinDown'))
                  END IF
748 749
                  evac0Temp(numVac,1) = eParamUp
                  IF(input%jspins.GT.1) evac0Temp(numVac,2) = eParamDown
Matthias Redies's avatar
Matthias Redies committed
750
                  IF(i == 1) THEN
751 752
                     evac0Temp(3-numVac,1) = eParamUp
                     IF(input%jspins.GT.1) evac0Temp(3-numVac,2) = eParamDown
753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841
                  END IF
               END DO
            END IF
         END IF

         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a1')
         IF (numberNodes.EQ.1) THEN
            latticeDef = 1
            input%scaleA1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1/@scale'))
            a1(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1'))
            a1(1) = a1(1) * input%scaleA1
            numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a2')
            IF (numberNodes.EQ.1) THEN
               latticeDef = 2
               input%scaleA2 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2/@scale'))
               a2(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2'))
               a2(2) = a2(2) * input%scaleA2
            END IF
            IF(.NOT.input%film) THEN
               input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale'))
               a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c'))
               a3(3) = a3(3) * input%scaleC
            END IF
         END IF

         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/row-1')
         IF (numberNodes.EQ.1) THEN
            latticeDef = 3
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-1')))
            a1(1) = evaluateFirst(valueString)
            a1(2) = evaluateFirst(valueString)
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-2')))
            a2(1) = evaluateFirst(valueString)
            a2(2) = evaluateFirst(valueString)
            IF(.NOT.input%film) THEN
               input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale'))
               a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c'))
               a3(3) = a3(3) * input%scaleC
            END IF
         END IF

         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix')
         IF (numberNodes.EQ.1) THEN
            latticeDef = 4
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-1')))
            a1(1) = evaluateFirst(valueString)
            a1(2) = evaluateFirst(valueString)
            IF(.NOT.input%film) THEN
               a1(3) = evaluateFirst(valueString)
            END IF
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-2')))
            a2(1) = evaluateFirst(valueString)
            a2(2) = evaluateFirst(valueString)
            IF(.NOT.input%film) THEN
               a2(3) = evaluateFirst(valueString)
            END IF
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-3')))
            IF(.NOT.input%film) THEN
               a3(1) = evaluateFirst(valueString)
               a3(2) = evaluateFirst(valueString)
               a3(3) = evaluateFirst(valueString)
            ELSE
               !WRITE(*,*) 'Note: For film calculations only the upper left 2x2 part of the Bravais matrix is considered.'
            END IF
         END IF
      END IF ! Note: Further ways to define lattices might be added later. (1D lattice,...)

      ! Construction of amat requires additional information about the lattice
      ! and is done later (scroll down)!

      ! Read in symmetry parameters

      sym%namgrp = 'any'

      xPathA = '/fleurInput/cell/symmetry'
      numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/symmetry')

      IF (numberNodes.EQ.1) THEN
         sym%symSpecType = 2
         symmetryDef = 1
         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spgrp')))
         READ(valueString,*) sym%namgrp
         sym%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs'))
         sym%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs'))
         sym%invs2 = sym%invs.AND.sym%zrfs

         IF (sym%namgrp.EQ.'any ') THEN
            sym%nop = 48
            ! Read in sym.out file if sym%namgrp='any' set.
842 843
            !CALL rw_symfile('r',94,'sym.out',48,cell%bmat,&
            !   &                        mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor)
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
            IF (ALLOCATED(sym%mrot)) THEN
               DEALLOCATE(sym%mrot)
            END IF
            ALLOCATE(sym%mrot(3,3,sym%nop))
            IF (ALLOCATED(sym%tau)) THEN
               DEALLOCATE(sym%tau)
            END IF
            ALLOCATE(sym%tau(3,sym%nop))

            DO k = 1, sym%nop
               DO i = 1, 3
                  DO j = 1, 3
                     sym%mrot(j,i,k) = mrotTemp(j,i,k)
                  END DO
                  sym%tau(i,k) = tauTemp(i,k)
               END DO
            END DO
         END IF
      END IF

      xPathA = '/fleurInput/cell/symmetryFile'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      IF (numberNodes.EQ.1) THEN
         symmetryDef = 2
         sym%symSpecType = 1
         sym%nop = 48
         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@filename')))

873 874
         !CALL rw_symfile('r',94,TRIM(ADJUSTL(valueString)),48,cell%bmat,&
         !   &                     mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor)
875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969

         IF (ALLOCATED(sym%mrot)) THEN
            DEALLOCATE(sym%mrot)
         END IF
         ALLOCATE(sym%mrot(3,3,sym%nop))
         IF (ALLOCATED(sym%tau)) THEN
            DEALLOCATE(sym%tau)
         END IF
         ALLOCATE(sym%tau(3,sym%nop))

         sym%invs = .FALSE.
         sym%zrfs = .FALSE.

         DO k = 1, sym%nop
            absSum = 0
            DO i = 1, 3
               DO j = 1, 3
                  sym%mrot(j,i,k) = mrotTemp(j,i,k)
                  absSum = absSum + ABS(sym%mrot(j,i,k))
               END DO
               sym%tau(i,k) = tauTemp(i,k)
            END DO
            IF (absSum.EQ.3) THEN
               IF (ALL(sym%tau(:,k).EQ.0.0)) THEN
                  IF ((sym%mrot(1,1,k).EQ.-1).AND.(sym%mrot(2,2,k).EQ.-1).AND.(sym%mrot(3,3,k).EQ.-1)) sym%invs = .TRUE.
                  IF ((sym%mrot(1,1,k).EQ.1).AND.(sym%mrot(2,2,k).EQ.1).AND.(sym%mrot(3,3,k).EQ.-1)) sym%zrfs = .TRUE.
               END IF
            END IF
         END DO

         sym%invs2 = sym%invs.AND.sym%zrfs
      END IF

      xPathA = '/fleurInput/cell/symmetryOperations'
      numberNodes = xmlGetNumberOfNodes(xPathA)

      IF (numberNodes.EQ.1) THEN
         sym%symSpecType = 3
         symmetryDef = 3

         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/symOp')
         sym%nop = numberNodes

         IF (ALLOCATED(sym%mrot)) DEALLOCATE(sym%mrot)
         ALLOCATE(sym%mrot(3,3,sym%nop))
         IF (ALLOCATED(sym%tau)) DEALLOCATE(sym%tau)
         ALLOCATE(sym%tau(3,sym%nop))
         sym%symor = .TRUE.
         DO i = 1, sym%nop
            WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-1'
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
            READ(valueString,*) sym%mrot(1,1,i), sym%mrot(1,2,i), sym%mrot(1,3,i), sym%tau(1,i)

            WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-2'
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
            READ(valueString,*) sym%mrot(2,1,i), sym%mrot(2,2,i), sym%mrot(2,3,i), sym%tau(2,i)

            WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-3'
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
            READ(valueString,*) sym%mrot(3,1,i), sym%mrot(3,2,i), sym%mrot(3,3,i), sym%tau(3,i)

            IF ((sym%tau(1,i)**2 + sym%tau(2,i)**2 + sym%tau(3,i)**2).GT.1.e-8) THEN
               sym%symor = .FALSE.
            END IF
            DO j = 1,3
               IF (ABS(sym%tau(j,i)-0.33333) < 0.00001) THEN
                  sym%tau(j,i) = 1./3.
               ENDIF
               IF (ABS(sym%tau(j,i)+0.33333) < 0.00001) THEN
                  sym%tau(j,i) = -1./3.
               ENDIF
               IF (ABS(sym%tau(j,i)-0.66667) < 0.00001) THEN
                  sym%tau(j,i) = 2./3.
               ENDIF
               IF (ABS(sym%tau(j,i)+0.66667) < 0.00001) THEN
                  sym%tau(j,i) = -2./3.
               ENDIF
            ENDDO
         END DO
      END IF

      ! Calculate missing symmetry and cell properties and check consistency of parameters.


      IF (latticeScale.EQ.0.0) latticeScale = 1.0
      IF (.NOT.input%film) vacuum%dvac = a3(3)
      vacuum%dvac = latticeScale*vacuum%dvac
      dtild = latticeScale*dtild

      cell%amat(:,1) = a1(:) * latticeScale
      cell%amat(:,2) = a2(:) * latticeScale
      cell%amat(:,3) = a3(:) * latticeScale

      CALL inv3(cell%amat,cell%bmat,cell%omtil)
      cell%bmat(:,:) = tpi_const*cell%bmat(:,:)
970
      cell%omtil = ABS(cell%omtil)
971 972 973 974 975 976 977 978 979 980 981 982

      IF (input%film.AND..NOT.oneD%odd%d1) THEN
         cell%vol = (cell%omtil/cell%amat(3,3))*vacuum%dvac
         cell%area = cell%omtil/cell%amat(3,3)
         !-odim
      ELSE IF (oneD%odd%d1) THEN
         cell%area = tpi_const*cell%amat(3,3)
         cell%vol = pi_const*(vacuum%dvac**2)*cell%amat(3,3)/4.0
         !+odim
      ELSE
         cell%vol = cell%omtil
         cell%area = cell%amat(1,1)*cell%amat(2,2)-cell%amat(1,2)*cell%amat(2,1)
Matthias Redies's avatar
Matthias Redies committed
983
         IF (cell%area < 1.0e-7) THEN
984 985 986 987 988 989
               cell%area = 1.
         END IF
      END IF

      ! Construction of missing symmetry information
      IF ((symmetryDef.EQ.2).OR.(symmetryDef.EQ.3)) THEN
990
         CALL sym%init(cell,input%film)
991
      END IF
992
 
993
      ALLOCATE (sym%invarop(atoms%nat,sym%nop),sym%invarind(atoms%nat))
994
      ALLOCATE (sym%invsatnr(atoms%nat))
995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009

      !some settings for film calculations
      vacuum%nmzd = 250
      vacuum%nmzxyd = 100
      vacuum%nvac = 2
      IF (sym%zrfs.OR.sym%invs) vacuum%nvac = 1
      IF (oneD%odd%d1) vacuum%nvac = 1
      cell%z1 = vacuum%dvac/2
      vacuum%nmz = vacuum%nmzd
      vacuum%delz = 25.0/vacuum%nmz
      IF (oneD%odd%d1) vacuum%delz = 20.0 / vacuum%nmz
      IF (vacuum%nmz.GT.vacuum%nmzd) CALL juDFT_error("nmzd",calledby ="inped")
      vacuum%nmzxy = vacuum%nmzxyd
      IF (vacuum%nmzxy.GT.vacuum%nmzxyd) CALL juDFT_error("nmzxyd",calledby ="inped")

1010 1011 1012
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! End of cell section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1013

1014 1015 1016
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Start of XC functional section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1017 1018 1019 1020 1021

      !Read in libxc parameters if present
      xPathA = '/fleurInput/xcFunctional/LibXCID'
      xPathB = '/fleurInput/xcFunctional/LibXCName'

1022 1023 1024
      IF(xmlGetNumberOfNodes(xPathA) == 1 .AND. xmlGetNumberOfNodes(xPathB) == 1) THEN
         CALL judft_error("LibXC is given both by Name and ID and is therefore overdetermined", calledby="r_inpXML")
      ENDIF
1025 1026 1027 1028 1029 1030 1031 1032


      ! LibXCID 
      IF (xmlGetNumberOfNodes(xPathA) == 1) THEN
#ifdef CPP_LIBXC
         vxc_id_x=evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@exchange'))
         vxc_id_c=evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@correlation'))

1033
         IF(xmlGetNumberOfNodes(TRIM(xPathA) // '/@etot_exchange') == 1) THEN
1034
            exc_id_x = evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@etot_exchange'))
1035
         ELSE
1036
            exc_id_x = vxc_id_x
1037
         ENDIF
1038
         
1039
         IF(xmlGetNumberOfNodes(TRIM(xPathA) // '/@exc_correlation') == 1) THEN
1040
            exc_id_c = evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@exc_correlation'))
1041
         ELSE
1042
            exc_id_c = vxc_id_c
1043
         ENDIF
1044
#else
1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055
         CALL judft_error("To use libxc functionals you have to compile with libXC support")
#endif
      ! LibXCName 
      ELSEIF (xmlGetNumberOfNodes(TRIM(xPathB)) == 1) THEN
#ifdef CPP_LIBXC
         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@exchange')))
         vxc_id_x =  xc_f03_functional_get_number(TRIM(valueString))

         valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@correlation')))
         vxc_id_c =  xc_f03_functional_get_number(TRIM(valueString))
         
Matthias Redies's avatar
Matthias Redies committed
1056 1057
         IF(xmlGetNumberOfNodes(TRIM(xPathB) // '/@etot_exchange') == 1) THEN
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@etot_exchange')))
1058
            exc_id_x =  xc_f03_functional_get_number(TRIM(valueString))
1059
         ELSE
1060
            exc_id_x = vxc_id_x
1061
         ENDIF
1062
         
Matthias Redies's avatar
Matthias Redies committed
1063 1064
         IF(xmlGetNumberOfNodes(TRIM(xPathB) // '/@etot_correlation') == 1) THEN
            valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@etot_correlation')))
1065
            exc_id_c =  xc_f03_functional_get_number(TRIM(valueString))
1066
         ELSE
1067
            exc_id_c = vxc_id_c
1068
         ENDIF
1069 1070 1071 1072
#else
         CALL judft_error("To use libxc functionals you have to compile with libXC support")
#endif
      ELSE
Matthias Redies's avatar
Matthias Redies committed
1073 1074
         vxc_id_x=0; vxc_id_c=0;
         exc_id_x=0; exc_id_c=0;
1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099
      ENDIF

      ! Read in xc functional parameters
      valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/@name')))))
      namex(1:4) = valueString(1:4)
      l_relcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/@relativisticCorrections'))

      relcor = 'non-relativi'
      IF (l_relcor) THEN
         relcor = 'relativistic'
      END IF

      !now initialize the xcpot variable
      CALL setXCParameters(atoms,valueString,l_relcor,input%jspins,vxc_id_x,vxc_id_c,exc_id_x, exc_id_c, xcpot)

      xPathA = '/fleurInput/calculationSetup/cutoffs/@GmaxXC'
      numberNodes = xmlGetNumberOfNodes(xPathA)
      xcpot%gmaxxc = stars%gmax
      IF(numberNodes.EQ.1) THEN
         xcpot%gmaxxc = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
      END IF
      hybrid%l_hybrid=xcpot%is_hybrid()

      ALLOCATE(hybrid%lcutm1(atoms%ntype),hybrid%lcutwf(atoms%ntype),hybrid%select1(4,atoms%ntype))

Daniel Wortmann's avatar
Daniel Wortmann committed
1100
  
1101 1102 1103 1104
      hybrid%gcutm1 = input%rkmax - 0.5
      hybrid%tolerance1 = 1.0e-4
      hybrid%ewaldlambda = 3
      hybrid%lexp = 16
1105
      hybrid%bands1 = DIMENSION%neigd
1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117

      numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/prodBasis')
      IF (numberNodes==0) THEN
         IF (hybrid%l_hybrid) CALL judft_error("Mixed product basis input missing in inp.xml")
      ELSE
         hybrid%gcutm1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@gcutm'))
         hybrid%tolerance1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@tolerance'))
         hybrid%ewaldlambda=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@ewaldlambda'))
         hybrid%lexp=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@lexp'))
         hybrid%bands1=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@bands'))
      ENDIF

1118 1119 1120
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! End of XC functional section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1121

1122 1123 1124
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!! Start of species section
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1125 1126

      ALLOCATE (speciesNLO(numSpecies))
Daniel Wortmann's avatar
Daniel Wortmann committed
1127 1128 1129
      
      ALLOCATE(speciesName(numSpecies))
      ALLOCATE(atoms%speciesName(atoms%ntype))
1130 1131 1132 1133 1134 1135 1136 1137 1138 1139

      atoms%lapw_l(:) = -1
      atoms%n_u = 0

      DEALLOCATE(noel)
      ALLOCATE(noel(atoms%ntype))

      DO iSpecies = 1, numSpecies
         ! Attributes of species
         WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']'
Daniel Wortmann's avatar
Daniel Wortmann committed
1140
         speciesName(iSpecies) = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name')))
1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169
         atomicNumber = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomicNumber'))
         magMom = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@magMom'))
         flipSpin = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@flipSpin'))

         ! Attributes of mtSphere element of species
         radius = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@radius'))
         gridPoints = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@gridPoints'))
         logIncrement = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@logIncrement'))

         ! Attributes of atomicCutoffs element of species
         lmax = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmax'))
         lnonsphr = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lnonsphr'))
         lmaxAPW = -1
         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW')
         IF (numberNodes.EQ.1) THEN
            lmaxAPW = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW'))
         END IF

         numU = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/ldaU')
        IF (numU.GT.4) CALL juDFT_error("Too many U parameters provided for a certain species (maximum is 4).",calledby ="r_inpXML")
         ldau_l = -1
         ldau_u = 0.0
         ldau_j = 0.0
         l_amf = .FALSE.
         DO i = 1, numU
            WRITE(xPathB,*) i
            ldau_l(i) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l'))
            ldau_u(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@U'))
            ldau_j(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@J'))
1170 1171 1172
            ldau_phi(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@phi'))
            ldau_theta(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@theta'))
            l_amf(i) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l_amf'))
1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185
         END DO

         speciesNLO(iSpecies) = 0
         WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo'
         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
         DO iLO = 1, numberNodes
            WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l'
            WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n'
            lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))
            nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC)))
            CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount)
            CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount)
            IF(lNumCount.NE.nNumCount) THEN
1186
               CALL judft_error('Error in LO input: l quantum number count does not equal n quantum number count')
1187 1188 1189 1190
            END IF
            speciesNLO(iSpecies) = speciesNLO(iSpecies) + lNumCount
            DEALLOCATE (lNumbers, nNumbers)
         END DO
1191 1192 1193 1194 1195 1196 1197 1198
         ! Special switches for species
         vcaspecies=0.0
         WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/special'
         numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
         IF (numberNodes==1) THEN
            vcaSpecies   = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vca_charge'))))
         ENDIF
       
1199 1200
         DO iType = 1, atoms%ntype
            WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species'
Daniel Wortmann's avatar
Daniel Wortmann committed
1201 1202
            atoms%speciesName(iType) = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
            IF(TRIM(ADJUSTL(speciesName(iSpecies))).EQ.TRIM(ADJUSTL(atoms%speciesName(iType)))) THEN
1203
               atoms%nz(iType) = atomicNumber
Daniel Wortmann's avatar
Daniel Wortmann committed
1204
               atoms%zatom(iType) = atoms%nz(iType)
1205 1206 1207 1208
               IF (atoms%nz(iType).EQ.0) THEN
                  WRITE(*,*) 'Note: Replacing atomic number 0 by 1.0e-10 on atom type ', iType
                  atoms%zatom(iType) = 1.0e-10
               END IF
Daniel Wortmann's avatar
Daniel Wortmann committed
1209
               atoms%zatom(iType)=atoms%zatom(iType)+vcaspecies
1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228
               noel(iType) = namat_const(atoms%nz(iType))
               atoms%rmt(iType) = radius
               atoms%jri(iType) = gridPoints
               atoms%dx(iType) = logIncrement
               atoms%lmax(iType) = lmax
               atoms%nlo(iType) = speciesNLO(iSpecies)
               atoms%lnonsph(iType) = lnonsphr
               atoms%lapw_l(iType) = lmaxAPW
               IF (flipSpin) THEN
                  atoms%nflip(iType) = 1
               ELSE
                  atoms%nflip(iType) = 0
               ENDIF
               atoms%bmu(iType) = magMom
               DO i = 1, numU
                  atoms%n_u = atoms%n_u + 1
                  atoms%lda_u(atoms%n_u)%l = ldau_l(i)
                  atoms%lda_u(atoms%n_u)%u = ldau_u(i)
                  atoms%lda_u(atoms%n_u)%j = ldau_j(i)
1229 1230
                  atoms%lda_u(atoms%n_u)%phi = ldau_phi(i)
                  atoms%lda_u(atoms%n_u)%theta = ldau_theta(i)
1231 1232 1233 1234 1235 1236 1237
                  atoms%lda_u(atoms%n_u)%l_amf = l_amf(i)
                  atoms%lda_u(atoms%n_u)%atomType = iType
               END DO
            END IF
         END DO
      END DO

1238
      atoms%lmaxd = MAXVAL(atoms%lmax(:))
1239 1240 1241
      atoms%llod  = 0
      atoms%nlod = 0
      DO iType = 1, atoms%ntype
1242
         atoms%nlod = MAX(atoms%nlod,atoms%nlo(iType))
1243
      END DO