LCOV - code coverage report
Current view: top level - io - r_inpXML.F90 (source / functions) Hit Total Coverage
Test: combined.info Lines: 961 1448 66.4 %
Date: 2019-09-08 04:53:50 Functions: 4 7 57.1 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------
       2             : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
       3             : ! This file is part of FLEUR and available as free software under the conditions
       4             : ! of the MIT license as expressed in the LICENSE file in more detail.
       5             : !--------------------------------------------------------------------------------
       6             : 
       7             : MODULE m_rinpXML
       8             :    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       9             :    !!!
      10             :    !!! The routine r_inpXML reads in the inp.xml file
      11             :    !!!
      12             :    !!!                               GM'16
      13             :    !!!
      14             :    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      15             : CONTAINS
      16          48 :    SUBROUTINE r_inpXML(&
      17             :       atoms,obsolete,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,xmlElectronStates,&
      20             :       xmlPrintCoreStates,xmlCoreOccs,atomTypeSpecies,speciesRepAtomType,&
      21             :       l_kpts)
      22             : 
      23             :       USE iso_c_binding
      24             :       USE m_juDFT
      25             :       USE m_types
      26             :       USE m_types_forcetheo_extended
      27             :       USE m_symdata , ONLY : nammap, ord2, l_c2
      28             :       USE m_rwsymfile
      29             :       USE m_xmlIntWrapFort
      30             :       USE m_inv3
      31             :       USE m_spg2set
      32             :       USE m_closure, ONLY : check_close
      33             :       USE m_symproperties
      34             :       USE m_calculator
      35             :       USE m_constants
      36             :       USE m_inpeig
      37             :       USE m_inpnoco
      38             :       USE m_sort
      39             :       USE m_types_xcpot_inbuild
      40             : #ifdef CPP_LIBXC
      41             :       USE xc_f03_lib_m
      42             : #endif
      43             :       IMPLICIT NONE
      44             : 
      45             :       TYPE(t_input),INTENT(INOUT)   :: input
      46             :       TYPE(t_sym),INTENT(INOUT)     :: sym
      47             :       TYPE(t_stars),INTENT(INOUT)   :: stars
      48             :       TYPE(t_atoms),INTENT(INOUT)   :: atoms
      49             :       TYPE(t_vacuum),INTENT(INOUT)   :: vacuum
      50             :       TYPE(t_obsolete),INTENT(INOUT) :: obsolete
      51             :       TYPE(t_kpts),INTENT(INOUT)     :: kpts
      52             :       TYPE(t_oneD),INTENT(INOUT)     :: oneD
      53             :       TYPE(t_hybrid),INTENT(INOUT)   :: hybrid
      54             :       TYPE(t_cell),INTENT(INOUT)     :: cell
      55             :       TYPE(t_banddos),INTENT(INOUT)  :: banddos
      56             :       TYPE(t_sliceplot),INTENT(INOUT):: sliceplot
      57             :       CLASS(t_xcpot),INTENT(INOUT),ALLOCATABLE :: xcpot
      58             :       TYPE(t_noco),INTENT(INOUT)     :: noco
      59             :       TYPE(t_dimension),INTENT(OUT)  :: dimension
      60             :       TYPE(t_enpara)   ,INTENT(OUT)  :: enpara
      61             :       TYPE(t_field), INTENT(INOUT)   :: field
      62             :       CLASS(t_forcetheo),ALLOCATABLE,INTENT(OUT):: forcetheo
      63             :       TYPE(t_coreSpecInput),INTENT(OUT) :: coreSpecInput
      64             :       TYPE(t_wann)   ,INTENT(INOUT)  :: wann
      65             :       LOGICAL, INTENT(OUT)           :: l_kpts
      66             :       INTEGER,          ALLOCATABLE, INTENT(INOUT) :: xmlElectronStates(:,:)
      67             :       INTEGER,          ALLOCATABLE, INTENT(INOUT) :: atomTypeSpecies(:)
      68             :       INTEGER,          ALLOCATABLE, INTENT(INOUT) :: speciesRepAtomType(:)
      69             :       REAL,             ALLOCATABLE, INTENT(INOUT) :: xmlCoreOccs(:,:,:)
      70             :       LOGICAL,          ALLOCATABLE, INTENT(INOUT) :: xmlPrintCoreStates(:,:)
      71             :       CHARACTER(len=3), ALLOCATABLE, INTENT(INOUT) :: noel(:)
      72             :       CHARACTER(len=4), INTENT(OUT)  :: namex
      73             :       CHARACTER(len=12), INTENT(OUT) :: relcor
      74             :       REAL, INTENT(OUT)              :: a1(3),a2(3),a3(3)
      75             :       REAL, INTENT(OUT)              :: dtild
      76             : 
      77             :       CHARACTER(len=8) :: name(10)
      78             : 
      79             :       !+odim
      80             :       INTEGER MM,vM,m_cyl
      81             :       LOGICAL invs1,zrfs1
      82             :       INTEGER chi,rot
      83             :       LOGICAL d1,band
      84             :       NAMELIST /odim/ d1,MM,vM,m_cyl,chi,rot,invs1,zrfs1
      85             :       !-odim
      86             :       ! ..
      87             :       ! ..  Local Variables
      88             :       REAL     :: scpos  ,zc
      89             :       INTEGER ieq,i,k,na,n,ii,vxc_id_c,vxc_id_x,exc_id_c,exc_id_x
      90             :       REAL s3,ah,a,hs2,rest,thetaj
      91             :       LOGICAL l_hyb,l_sym,ldum
      92             :       INTEGER :: ierr
      93             :       ! ..
      94             :       !...  Local Arrays
      95             :       !   CHARACTER :: helpchar(atoms%ntype)
      96             :       CHARACTER(len=  4) :: chntype
      97             :       CHARACTER(len= 41) :: chform
      98             :       CHARACTER(len=100) :: line
      99             : 
     100             :       CHARACTER(len=20) :: tempNumberString
     101             :       CHARACTER(len=150) :: format
     102             :       CHARACTER(len=20) :: mixingScheme
     103             :       CHARACTER(len=10) :: loType
     104             :       LOGICAL :: kptGamma, l_relcor,ldummy
     105             :       INTEGER :: iAtomType, startCoreStates, endCoreStates
     106             :       CHARACTER(len=100) :: xPosString, yPosString, zPosString
     107             :       CHARACTER(len=200) :: coreStatesString
     108             :       !   REAL :: tempTaual(3,atoms%nat)
     109             :       REAL             :: coreStateOccs(29,2)
     110             :       INTEGER          :: coreStateNprnc(29), coreStateKappa(29)
     111             :       INTEGER          :: speciesXMLElectronStates(29)
     112             :       REAL             :: speciesXMLCoreOccs(2,29)
     113             :       LOGICAL          :: speciesXMLPrintCoreStates(29)
     114             : 
     115             :       INTEGER            :: iType, iLO, iSpecies, lNumCount, nNumCount, iLLO, jsp, j, l, absSum, numTokens
     116             :       INTEGER            :: numberNodes, nodeSum, numSpecies, n2spg, n1, n2, ikpt, iqpt
     117             :       INTEGER            :: atomicNumber, coreStates, gridPoints, lmax, lnonsphr, lmaxAPW
     118             :       INTEGER            :: latticeDef, symmetryDef, nop48, firstAtomOfType, errorStatus
     119             :       INTEGER            :: loEDeriv, ntp1, ios, ntst, jrc, minNeigd, providedCoreStates, providedStates
     120             :       INTEGER            :: nv, nv2, kq1, kq2, kq3, nprncTemp, kappaTemp, tempInt
     121             :       INTEGER            :: ldau_l(4), numVac, numU
     122             :       INTEGER            :: speciesEParams(0:3)
     123             :       INTEGER            :: mrotTemp(3,3,48)
     124             :       REAL               :: tauTemp(3,48)
     125             :       REAL               :: bk(3)
     126             :       LOGICAL            :: flipSpin, l_eV, invSym, l_qfix, relaxX, relaxY, relaxZ
     127             :       LOGICAL            :: coreConfigPresent, l_enpara, l_orbcomp, tempBool, l_nocoinp
     128             :       REAL               :: magMom, radius, logIncrement, qsc(3), latticeScale, dr
     129             :       REAL               :: aTemp, zp, rmtmax, sumWeight, ldau_u(4), ldau_j(4), tempReal
     130             :       REAL               :: ldau_phi(4),ldau_theta(4)
     131             :       REAL               :: weightScale, eParamUp, eParamDown
     132             :       LOGICAL            :: l_amf(4)
     133             :       REAL, PARAMETER    :: boltzmannConst = 3.1668114e-6 ! value is given in Hartree/Kelvin
     134             :       INTEGER            :: lcutm,lcutwf,hybSelect(4)
     135             :       REAL               :: evac0Temp(2,2)
     136             : 
     137             :       CHARACTER(LEN=200,KIND=c_char) :: schemaFilename, docFilename
     138             :       CHARACTER(LEN=255) :: valueString, lString, nString, token
     139             :       CHARACTER(LEN=255) :: xPathA, xPathB, xPathC, xPathD, xPathE
     140             :       CHARACTER(LEN=11)  :: latticeType
     141             :       CHARACTER(LEN=50)  :: versionString
     142             :       CHARACTER(LEN=150) :: kPointsPrefix
     143             : 
     144             :       INTEGER            :: altKPointSetIndex,  altKPointSetIndices(2)
     145             :       LOGICAL            :: ldaSpecies
     146             :       REAL               :: socscaleSpecies,b_field_mtspecies,vcaspecies
     147             : 
     148          72 :       INTEGER, ALLOCATABLE :: lNumbers(:), nNumbers(:), speciesLLO(:)
     149          24 :       INTEGER, ALLOCATABLE :: loOrderList(:)
     150          24 :       INTEGER, ALLOCATABLE :: speciesNLO(:)
     151          72 :       INTEGER, ALLOCATABLE :: multtab(:,:), invOps(:), optype(:)
     152             :       INTEGER, ALLOCATABLE :: lmx1(:), nq1(:), nlhtp1(:)
     153          24 :       INTEGER, ALLOCATABLE :: speciesLOEDeriv(:)
     154          48 :       REAL,    ALLOCATABLE :: speciesLOeParams(:), speciesLLOReal(:)
     155          24 :       LOGICAL, ALLOCATABLE :: wannAtomList(:)
     156             : 
     157             :       ! Variables for MT radius testing:
     158             : 
     159             :       REAL                 :: dtild1,kmax1,dvac1
     160             :       LOGICAL              :: l_test
     161             :       INTEGER, ALLOCATABLE :: jri1(:), lmax1(:)
     162             :       REAL, ALLOCATABLE    :: rmt1(:), dx1(:)
     163             : 
     164             :       EXTERNAL prp_xcfft_box
     165             : 
     166             :       INTERFACE
     167             :          FUNCTION dropInputSchema() BIND(C, name="dropInputSchema")
     168             :             USE iso_c_binding
     169             :             INTEGER(c_int) dropInputSchema
     170             :          END FUNCTION dropInputSchema
     171             :       END INTERFACE
     172             : 
     173          24 :       errorStatus = 0
     174          24 :       errorStatus = dropInputSchema()
     175          24 :       IF(errorStatus.NE.0) THEN
     176           0 :          CALL juDFT_error('Error: Cannot print out FleurInputSchema.xsd')
     177             :       END IF
     178             : 
     179          24 :       schemaFilename = "FleurInputSchema.xsd"//C_NULL_CHAR
     180          24 :       docFilename = "inp.xml"//C_NULL_CHAR
     181             : 
     182             :       !TODO! these switches should be in the inp-file
     183          24 :       input%l_core_confpot=.TRUE. !former CPP_CORE
     184          24 :       input%l_useapw=.FALSE.   !former CPP_APW
     185             :       !WRITE(*,*) 'Start reading of inp.xml file'
     186          24 :       CALL xmlInitInterface()
     187          24 :       CALL xmlParseSchema(schemaFilename)
     188          24 :       CALL xmlParseDoc(docFilename)
     189          24 :       CALL xmlValidateDoc()
     190          24 :       CALL xmlInitXPath()
     191             : 
     192             :       ! Check version of inp.xml
     193          24 :       versionString = xmlGetAttributeValue('/fleurInput/@fleurInputVersion')
     194             :       IF((TRIM(ADJUSTL(versionString)).NE.'0.27').AND.(TRIM(ADJUSTL(versionString)).NE.'0.28').AND.&
     195          24 :          (TRIM(ADJUSTL(versionString)).NE.'0.29').AND.(TRIM(ADJUSTL(versionString)).NE.'0.30')) THEN
     196           0 :          CALL juDFT_error('version number of inp.xml file is not compatible with this fleur version')
     197             :       END IF
     198             : 
     199             :       ! Get number of atoms, atom types, and atom species
     200             : 
     201          24 :       numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/relPos')
     202          24 :       numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/absPos')
     203          24 :       numberNodes = numberNodes + xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup/filmPos')
     204             : 
     205          24 :       atoms%nat = numberNodes
     206             : 
     207          24 :       numberNodes = xmlGetNumberOfNodes('/fleurInput/atomGroups/atomGroup')
     208             : 
     209          24 :       atoms%ntype = numberNodes
     210             : 
     211          24 :       numSpecies = xmlGetNumberOfNodes('/fleurInput/atomSpecies/species')
     212             : 
     213          24 :       ALLOCATE(atoms%nz(atoms%ntype))     !nz and zatom have the same content!
     214          24 :       ALLOCATE(atoms%zatom(atoms%ntype))  !nz and zatom have the same content!
     215          24 :       ALLOCATE(atoms%jri(atoms%ntype))
     216          24 :       ALLOCATE(atoms%dx(atoms%ntype))
     217          24 :       ALLOCATE(atoms%lmax(atoms%ntype))
     218          24 :       ALLOCATE(atoms%nlo(atoms%ntype))
     219          24 :       ALLOCATE(atoms%ncst(atoms%ntype))
     220          24 :       ALLOCATE(atoms%lnonsph(atoms%ntype))
     221          24 :       ALLOCATE(atoms%nflip(atoms%ntype))
     222          24 :       ALLOCATE(atoms%l_geo(atoms%ntype))
     223          24 :       ALLOCATE(atoms%lda_u(4*atoms%ntype))
     224          24 :       ALLOCATE(atoms%bmu(atoms%ntype))
     225          24 :       ALLOCATE(atoms%relax(3,atoms%ntype))
     226          24 :       ALLOCATE(atoms%neq(atoms%ntype))
     227          24 :       ALLOCATE(atoms%taual(3,atoms%nat))
     228          24 :       ALLOCATE(atoms%label(atoms%nat))
     229          24 :       ALLOCATE(atoms%pos(3,atoms%nat))
     230          24 :       ALLOCATE(atoms%rmt(atoms%ntype))
     231          24 :       ALLOCATE(atoms%numStatesProvided(atoms%ntype))
     232          24 :       ALLOCATE(atoms%namex(atoms%ntype))
     233          24 :       ALLOCATE(atoms%icorr(atoms%ntype))
     234          24 :       ALLOCATE(atoms%igrd(atoms%ntype))
     235          24 :       ALLOCATE(atoms%krla(atoms%ntype))
     236          24 :       ALLOCATE(atoms%relcor(atoms%ntype))
     237             : 
     238          79 :       atoms%namex = ''
     239          24 :       atoms%icorr = -99
     240             : 
     241          24 :       ALLOCATE(atoms%ncv(atoms%ntype)) ! For what is this?
     242          24 :       ALLOCATE(atoms%ngopr(atoms%nat)) ! For what is this?
     243          24 :       ALLOCATE(atoms%lapw_l(atoms%ntype)) ! Where do I put this?
     244          24 :       ALLOCATE(atoms%invsat(atoms%nat)) ! Where do I put this?
     245             : 
     246          24 :       ALLOCATE(noco%l_relax(atoms%ntype),noco%b_con(2,atoms%ntype))
     247          24 :       ALLOCATE(noco%alphInit(atoms%ntype),noco%alph(atoms%ntype),noco%beta(atoms%ntype))
     248          24 :       ALLOCATE(noco%socscale(atoms%ntype))
     249             : 
     250          24 :       DEALLOCATE(atomTypeSpecies,speciesRepAtomType)
     251          24 :       ALLOCATE(atomTypeSpecies(atoms%ntype))
     252          24 :       ALLOCATE(speciesRepAtomType(numSpecies))
     253          79 :       atomTypeSpecies = -1
     254          57 :       speciesRepAtomType = -1
     255             : 
     256          24 :       DEALLOCATE(xmlElectronStates,xmlPrintCoreStates,xmlCoreOccs)
     257          24 :       ALLOCATE(xmlElectronStates(29,atoms%ntype))
     258          24 :       ALLOCATE(xmlPrintCoreStates(29,atoms%ntype))
     259          24 :       ALLOCATE(xmlCoreOccs(2,29,atoms%ntype))
     260          79 :       xmlElectronStates = noState_const
     261          79 :       xmlPrintCoreStates = .FALSE.
     262          79 :       xmlCoreOccs = 0.0
     263             : 
     264          24 :       ALLOCATE (kpts%ntetra(4,kpts%ntet),kpts%voltet(kpts%ntet))
     265             : 
     266          24 :       ALLOCATE (wannAtomList(atoms%nat))
     267             : 
     268             :       ! Read in constants
     269             : 
     270          24 :       xPathA = '/fleurInput/constants/constant'
     271          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     272          24 :       DO i = 1, numberNodes
     273           0 :          WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)), '[',i,']'
     274           0 :          tempReal = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@value'))
     275           0 :          valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@name')
     276          24 :          CALL ASSIGN_var(valueString,tempReal)
     277             :       END DO
     278             : 
     279             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     280             : !!! Comment section
     281             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     282             : 
     283         264 :       input%comment = '        '
     284          24 :       xPathA = '/fleurInput/comment'
     285          24 :       valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
     286        1630 :       DO i = 1, LEN(TRIM(ADJUSTL(valueString)))
     287        1630 :          IF (valueString(i:i) == ACHAR(10)) valueString(i:i) = ' ' !remove line breaks
     288             :       END DO
     289          24 :       valueString = TRIM(ADJUSTL(valueString))
     290         264 :       DO i = 1, 10
     291         240 :          j = (i-1) * 8 + 1
     292         264 :          input%comment(i) = valueString(j:j+7)
     293             :       END DO
     294             : 
     295             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     296             : !!! Start of calculationSetup section
     297             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     298             : 
     299             :       ! Read general cutoff parameters
     300             : 
     301          24 :       input%rkmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Kmax'))
     302          24 :       stars%gmax = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/cutoffs/@Gmax'))
     303             : 
     304          24 :       stars%gmaxInit = stars%gmax
     305             : 
     306          24 :       xPathA = '/fleurInput/calculationSetup/cutoffs/@numbands'
     307          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     308          24 :       DIMENSION%neigd = 0
     309          24 :       IF(numberNodes.EQ.1) THEN
     310          24 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
     311          24 :          IF(TRIM(ADJUSTL(valueString)).EQ.'all') THEN
     312           0 :             dimension%neigd = -1
     313             :          ELSE
     314          24 :             READ(valueString,*) DIMENSION%neigd
     315             :          END IF
     316             :       END IF
     317             : 
     318             :       ! Read SCF loop parametrization
     319             : 
     320          24 :       input%itmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@itmax'))
     321          24 :       input%minDistance = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@minDistance'))
     322          24 :       input%maxiter = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@maxIterBroyd'))
     323             : 
     324          24 :       valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@imix')))
     325           0 :       SELECT CASE (valueString)
     326             :       CASE ('straight')
     327           0 :          input%imix = 0
     328             :       CASE ('Broyden1')
     329           0 :          input%imix = 3
     330             :       CASE ('Broyden2')
     331           0 :          input%imix = 5
     332             :       CASE ('Anderson')
     333          24 :          input%imix = 7
     334             :       CASE ("Pulay")
     335           0 :          input%imix = 9
     336             :       CASE ("pPulay")
     337           0 :          input%imix = 11
     338             :       CASE ("rPulay")
     339           0 :          input%imix = 13
     340             :       CASE ("aPulay")
     341           0 :          input%imix = 15
     342             :       CASE DEFAULT
     343          24 :          CALL juDFT_error('Error: unknown mixing scheme selected!')
     344             :       END SELECT
     345             : 
     346          24 :       input%alpha = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@alpha'))
     347          24 :       input%preconditioning_param = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@precondParam'))
     348          24 :       input%spinf = evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/scfLoop/@spinf'))
     349             : 
     350             :       ! Get parameters for core electrons
     351             : 
     352          24 :       input%ctail = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@ctail'))
     353          24 :       IF((TRIM(ADJUSTL(versionString)).EQ.'0.27')) THEN
     354          13 :          input%coretail_lmax = 99
     355             :       ELSE
     356          11 :        input%coretail_lmax = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@coretail_lmax'))
     357             :       END IF
     358          24 :       input%frcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@frcor'))
     359          24 :       input%kcrel = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/coreElectrons/@kcrel'))
     360             : 
     361             :       ! Read in magnetism parameters
     362             : 
     363          24 :       input%jspins = evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@jspins'))
     364          24 :       noco%l_noco = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@l_noco'))
     365          24 :       input%swsp = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@swsp'))
     366          24 :       input%lflip = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@lflip'))
     367          24 :       input%fixed_moment=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/magnetism/@fixed_moment'))
     368             : 
     369          24 :   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")
     370             : 
     371             :       ! Read in optional expert modes switches
     372             : 
     373          24 :       xPathA = '/fleurInput/calculationSetup/expertModes'
     374          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     375             : 
     376          24 :       input%gw = 0
     377          24 :       input%secvar = .FALSE.
     378             : 
     379          24 :       IF (numberNodes.EQ.1) THEN
     380          24 :          input%gw = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gw'))
     381          24 :          input%secvar = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@secvar'))
     382             :       END IF
     383             : 
     384             :       ! Check for alternative k point sets for the chosen FLEUR mode
     385             : 
     386          24 :       xPathA = '/fleurInput/output'
     387          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     388          24 :       banddos%band = .FALSE.
     389          24 :       IF (numberNodes.EQ.1) THEN
     390          24 :          banddos%band = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@band'))
     391             :       END IF
     392             : 
     393          72 :       altKPointSetIndices(:) = -1
     394          24 :       xPathA = '/fleurInput/calculationSetup/bzIntegration/altKPointSet'
     395          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     396          24 :       IF(numberNodes.NE.0) THEN
     397           0 :          DO i = 1, numberNodes
     398           0 :             WRITE(xPathA,*) '/fleurInput/calculationSetup/bzIntegration/altKPointSet[',i,']/@purpose'
     399           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
     400           0 :             IF((altKPointSetIndices(2).EQ.-1).AND.(TRIM(ADJUSTL(valueString)).EQ.'GW')) THEN
     401           0 :                altKPointSetIndices(2) = i
     402           0 :             ELSE IF((altKPointSetIndices(1).EQ.-1).AND.(TRIM(ADJUSTL(valueString)).EQ.'bands')) THEN
     403           0 :                altKPointSetIndices(1) = i
     404             :             END IF
     405             :          END DO
     406             :       END IF
     407             : 
     408          24 :       altKPointSetIndex = -1
     409          24 :       IF(banddos%band) THEN
     410           1 :          altKPointSetIndex = altKPointSetIndices(1)
     411          23 :       ELSE IF (input%gw.EQ.2) THEN
     412           0 :          altKPointSetIndex = altKPointSetIndices(2)
     413             :       END IF
     414             : 
     415          24 :       IF (altKPointSetIndex.NE.-1) THEN
     416           0 :          WRITE(kPointsPrefix,*) '/fleurInput/calculationSetup/bzIntegration/altKPointSet[',altKPointSetIndex,']'
     417             :       END IF
     418             : 
     419             :       ! Read in Brillouin zone integration parameters
     420             : 
     421          96 :       kpts%nkpt3 = 0
     422          24 :       l_kpts = .FALSE.
     423             : 
     424          24 :       valueString = TRIM(ADJUSTL(xmlGetAttributeValue('/fleurInput/calculationSetup/bzIntegration/@mode')))
     425          23 :       SELECT CASE (valueString)
     426             :       CASE ('hist')
     427          23 :          input%gauss = .FALSE.
     428          23 :          input%tria = .FALSE.
     429             :       CASE ('gauss')
     430           0 :          input%gauss = .TRUE.
     431           0 :          input%tria = .FALSE.
     432             :       CASE ('tria')
     433           1 :          input%gauss = .FALSE.
     434           1 :          input%tria = .TRUE.
     435             :       CASE DEFAULT
     436          24 :          CALL juDFT_error('Invalid bzIntegration mode selected!')
     437             :       END SELECT
     438             : 
     439          24 :       nodeSum = 0
     440          24 :       xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingEnergy'
     441          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     442          24 :       nodeSum = nodeSum + numberNodes
     443          24 :       IF (numberNodes.EQ.1) THEN
     444          24 :          input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
     445             :       END IF
     446          24 :       xPathA = '/fleurInput/calculationSetup/bzIntegration/@fermiSmearingTemp'
     447          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     448          24 :       nodeSum = nodeSum + numberNodes
     449          24 :       IF (numberNodes.EQ.1) THEN
     450           0 :          input%tkb = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
     451           0 :          input%tkb = boltzmannConst * input%tkb
     452             :       END IF
     453          24 :       IF(nodeSum.GE.2) THEN
     454           0 :          CALL juDFT_error('Error: Multiple fermi Smearing parameters provided in input file!')
     455             :       END IF
     456             : 
     457          24 :       xPathA = '/fleurInput/calculationSetup/bzIntegration/@valenceElectrons'
     458          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     459          24 :       IF (numberNodes.EQ.1) THEN
     460          24 :          input%zelec = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
     461             :       ELSE
     462           0 :          CALL juDFT_error('Error: Optionality of valence electrons in input file not yet implemented!')
     463             :       END IF
     464             : 
     465          24 :       IF (altKPointSetIndex.EQ.-1) THEN
     466          24 :          WRITE(kPointsPrefix,*) '/fleurInput/calculationSetup/bzIntegration'
     467             :       END IF
     468             : 
     469             :       ! Option kPointDensity
     470          96 :       kpts%kPointDensity(:) = 0.0
     471          24 :       xPathA = TRIM(ADJUSTL(kPointsPrefix))//'/kPointDensity'
     472          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     473          24 :       IF (numberNodes.EQ.1) THEN
     474           0 :          l_kpts = .FALSE.
     475           0 :          kpts%kPointDensity(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denX'))
     476           0 :          kpts%kPointDensity(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denY'))
     477           0 :          kpts%kPointDensity(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denZ'))
     478           0 :          kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma'))
     479           0 :          kpts%specificationType = 4
     480             :       END IF
     481             : 
     482             :       ! Option kPointMesh
     483          24 :       xPathA = TRIM(ADJUSTL(kPointsPrefix))//'/kPointMesh'
     484          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     485          24 :       IF (numberNodes.EQ.1) THEN
     486           6 :          l_kpts = .FALSE.
     487           6 :          kpts%nkpt3(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nx'))
     488           6 :          kpts%nkpt3(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ny'))
     489           6 :          kpts%nkpt3(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nz'))
     490           6 :          kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma'))
     491           6 :          kpts%nkpt = kpts%nkpt3(1) * kpts%nkpt3(2) * kpts%nkpt3(3)
     492           6 :          kpts%specificationType = 2
     493             :       END IF
     494             : 
     495             :       ! Option kPointCount
     496          24 :       xPathA = TRIM(ADJUSTL(kPointsPrefix))//'/kPointCount'
     497          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     498          24 :       IF (numberNodes.EQ.1) THEN
     499          17 :          l_kpts = .FALSE.
     500          17 :          kpts%nkpt = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@count'))
     501          17 :          kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma'))
     502          17 :          kpts%nkpt = kpts%nkpt
     503          17 :          kpts%specificationType = 1
     504             : 
     505          17 :          ALLOCATE(kpts%bk(3,kpts%nkpt))
     506          17 :          ALLOCATE(kpts%wtkpt(kpts%nkpt))
     507         151 :          kpts%bk = 0.0
     508         151 :          kpts%wtkpt = 0.0
     509          17 :          kpts%posScale = 1.0
     510             : 
     511          17 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(kPointsPrefix))//'/kPointCount/specialPoint')
     512          17 :          IF(numberNodes.EQ.1) THEN
     513           0 :             CALL juDFT_error('Error: Single special k point provided. This does not make sense!')
     514             :          END IF
     515          17 :          kpts%numSpecialPoints = numberNodes
     516          17 :          IF(kpts%numSpecialPoints.GE.2) THEN
     517           1 :             DEALLOCATE(kpts%specialPoints)
     518           1 :             ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints))
     519           1 :             ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints))
     520           3 :             DO i = 1, kpts%numSpecialPoints
     521           2 :                WRITE(xPathA,*) TRIM(ADJUSTL(kPointsPrefix))//'/kPointCount/specialPoint[',i,']'
     522           2 :                valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
     523           2 :                kpts%specialPoints(1,i) = evaluatefirst(valueString)
     524           2 :                kpts%specialPoints(2,i) = evaluatefirst(valueString)
     525           2 :                kpts%specialPoints(3,i) = evaluatefirst(valueString)
     526           3 :                kpts%specialPointNames(i) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name')
     527             :             END DO
     528             :          END IF
     529             :       ELSE
     530           7 :          DEALLOCATE(kpts%specialPoints)
     531           7 :          ALLOCATE(kpts%specialPoints(3,kpts%numSpecialPoints))
     532           7 :          ALLOCATE(kpts%specialPointNames(kpts%numSpecialPoints))
     533             :       END IF
     534             : 
     535             :       ! Option kPointList
     536          24 :       numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList')
     537          24 :       IF (numberNodes.EQ.1) THEN
     538           1 :          l_kpts = .TRUE.
     539           1 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/kPoint')
     540           1 :          kpts%nkpt = numberNodes
     541           1 :          kpts%l_gamma = .FALSE.
     542           1 :          ALLOCATE(kpts%bk(3,kpts%nkpt))
     543           1 :          ALLOCATE(kpts%wtkpt(kpts%nkpt))
     544           2 :          kpts%bk = 0.0
     545           2 :          kpts%wtkpt = 0.0
     546           1 :          kpts%specificationType = 3
     547             : 
     548           1 :          kpts%posScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/@posScale'))
     549           1 :          weightScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/@weightScale'))
     550             : 
     551           2 :          DO i = 1, kpts%nkpt
     552           1 :             WRITE(xPathA,*) TRIM(ADJUSTL(kPointsPrefix))//'/kPointList/kPoint[',i,']'
     553           1 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
     554           1 :             READ(valueString,*) kpts%bk(1,i), kpts%bk(2,i), kpts%bk(3,i)
     555           1 :             kpts%bk(:,i)=kpts%bk(:,i)/kpts%posScale
     556           1 :             kpts%wtkpt(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@weight'))
     557           2 :             kpts%wtkpt(i) = kpts%wtkpt(i) / weightScale
     558             :          END DO
     559           1 :          kpts%posScale=1.0
     560             :       END IF
     561             : 
     562             :       ! Option kPointListFile
     563          24 :       xPathA = TRIM(ADJUSTL(kPointsPrefix))//'/kPointListFile'
     564          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     565          24 :       IF (numberNodes.EQ.1) THEN
     566           0 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@filename')))
     567           0 :          OPEN (41,file=TRIM(ADJUSTL(valueString)),form='formatted',status='old')
     568           0 :             READ (41,*) kpts%nkpt
     569           0 :          CLOSE (41)
     570           0 :          ALLOCATE(kpts%bk(3,kpts%nkpt))
     571           0 :          ALLOCATE(kpts%wtkpt(kpts%nkpt))
     572           0 :          kpts%bk = 0.0
     573           0 :          kpts%wtkpt = 0.0
     574           0 :          kpts%l_gamma = .FALSE.
     575           0 :          l_kpts = .TRUE.
     576           0 :          kpts%specificationType = 3
     577           0 :          kpts%posScale=1.0
     578             : 
     579             :          ! We need to set input%film for inpeig. Unfortunately this is actually initialized at a later stage.
     580             :          ! So we do it here additionally.
     581           0 :          input%film = .FALSE.
     582           0 :          numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/filmLattice')
     583           0 :          IF (numberNodes.EQ.1) input%film = .TRUE.
     584             : 
     585           0 :          CALL inpeig(atoms,cell,input,.FALSE.,kpts,kptsFilename=TRIM(ADJUSTL(valueString)))
     586             :       END IF
     587             : 
     588             :       ! Read in optional SOC parameters if present
     589             : 
     590          24 :       xPathA = '/fleurInput/calculationSetup/soc'
     591          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     592             : 
     593          24 :       noco%l_soc = .FALSE.
     594          24 :       noco%theta = 0.0
     595          24 :       noco%phi = 0.0
     596             : 
     597          24 :       IF (numberNodes.EQ.1) THEN
     598          24 :          noco%theta=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta'))
     599          24 :          noco%phi=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi'))
     600          24 :          noco%l_soc = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_soc'))
     601          24 :          noco%l_spav = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spav'))
     602             :       END IF
     603             : 
     604             :       ! Read in optional noco parameters if present
     605             : 
     606          24 :       xPathA = '/fleurInput/calculationSetup/nocoParams'
     607          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     608             : 
     609          24 :       noco%l_ss = .FALSE.
     610          24 :       noco%l_mperp = .FALSE.
     611          24 :       noco%l_constr = .FALSE.
     612          24 :       noco%mix_b = 0.0
     613          96 :       noco%qss = 0.0
     614             : 
     615          24 :       noco%l_relax(:) = .FALSE.
     616          24 :       noco%alphInit(:) = 0.0
     617          24 :       noco%alph(:) = 0.0
     618          24 :       noco%beta(:) = 0.0
     619          24 :       noco%b_con(:,:) = 0.0
     620             : 
     621          24 :       IF ((noco%l_noco).AND.(numberNodes.EQ.0)) THEN
     622           0 :          CALL juDFT_error('Error: l_noco is true but no noco parameters set in xml input file!')
     623             :       END IF
     624             : 
     625          24 :       IF (numberNodes.EQ.1) THEN
     626           7 :          noco%l_ss = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_ss'))
     627           7 :          noco%l_mperp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mperp'))
     628           7 :          noco%l_constr = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_constr'))
     629           7 :          noco%l_mtNocoPot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_mtNocoPot'))
     630             : 
     631           7 :          noco%mix_b = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mix_b'))
     632             : 
     633           7 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qss')))
     634           7 :          READ(valueString,*) noco%qss(1), noco%qss(2), noco%qss(3)
     635             : 
     636             :          !WRITE(*,*) 'Note: TODO: Calculation of q points!'
     637             : 
     638           7 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/qsc')
     639           7 :          IF (numberNodes.EQ.1) THEN
     640           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/qsc')))
     641           0 :             READ(valueString,*) qsc(1), qsc(2), qsc(3)
     642           0 :             DO i = 1, 3
     643           0 :                noco%qss(i) = noco%qss(i) / qsc(i)
     644             :             END DO
     645             :             !WRITE(*,*) 'Note: TODO: Integrate qsc directly into qss in input file!'
     646             :             !WRITE(*,*) '(no problem for users)'
     647             :          END IF
     648             :       END IF
     649             : 
     650             :       ! Read in optional 1D parameters if present
     651             : 
     652          24 :       xPathA = '/fleurInput/calculationSetup/oneDParams'
     653          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     654             : 
     655          24 :       oneD%odd%d1 = .FALSE.
     656             : 
     657          24 :       IF (numberNodes.EQ.1) THEN
     658           0 :          oneD%odd%d1 = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@d1'))
     659           0 :          oneD%odd%M = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@MM'))
     660           0 :          oneD%odd%mb = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vM'))
     661           0 :          oneD%odd%m_cyl = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@m_cyl'))
     662           0 :          oneD%odd%chi = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@chi'))
     663           0 :          oneD%odd%rot = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@rot'))
     664           0 :          oneD%odd%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs1'))
     665           0 :          oneD%odd%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs1'))
     666             :       END IF
     667             : 
     668             :       ! Read in optional geometry optimization parameters
     669             : 
     670          24 :       xPathA = '/fleurInput/calculationSetup/geometryOptimization'
     671          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     672             : 
     673          24 :       input%l_f = .FALSE.
     674          24 :       input%qfix = 0
     675          24 :       input%forcemix = 2 ! BFGS is default.
     676             : 
     677          24 :       IF (numberNodes.EQ.1) THEN
     678           1 :          input%l_f = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_f'))
     679           1 :          input%forcealpha = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@forcealpha'))
     680           1 :          input%epsdisp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsdisp'))
     681           1 :          input%epsforce = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@epsforce'))
     682             : 
     683           1 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@forcemix')))
     684           0 :          SELECT CASE (valueString)
     685             :          CASE ('straight')
     686           0 :             input%forcemix = 0
     687             :          CASE ('CG')
     688           0 :             input%forcemix = 1
     689             :          CASE ('BFGS')
     690           1 :             input%forcemix = 2
     691             :          CASE DEFAULT
     692           1 :             CALL juDFT_error('Error: unknown force mixing scheme selected!', calledby='r_inpXML')
     693             :          END SELECT
     694             : 
     695           1 :          input%force_converged = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@force_converged'))
     696           1 :          input%qfix = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@qfix'))
     697             :       END IF
     698             : 
     699             :       ! Read in optional general LDA+U parameters
     700             : 
     701          24 :       xPathA = '/fleurInput/calculationSetup/ldaU'
     702          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     703          24 :       IF (numberNodes.EQ.1) THEN
     704           4 :          input%ldauLinMix = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_linMix'))
     705           4 :          input%ldauMixParam = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mixParam'))
     706           4 :          input%ldauSpinf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spinf'))
     707             :       END IF
     708             : 
     709             :       ! Read in RDMFT parameters
     710             : 
     711          24 :       input%l_rdmft = .FALSE.
     712          24 :       input%rdmftOccEps = 0.00001
     713          24 :       input%rdmftStatesBelow = 5
     714          24 :       input%rdmftStatesAbove = 5
     715          24 :       input%rdmftFunctional = -1
     716             : 
     717          24 :       xPathA = '/fleurInput/calculationSetup/rdmft'
     718          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     719          24 :       IF (numberNodes.EQ.1) THEN
     720           0 :          input%l_rdmft = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@l_rdmft'))
     721           0 :          input%rdmftOccEps = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@occEps'))
     722           0 :          input%rdmftStatesBelow = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@statesBelow'))
     723           0 :          input%rdmftStatesAbove = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@statesAbove'))
     724           0 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@functional')))
     725           0 :          SELECT CASE (valueString)
     726             :             CASE ('Muller')
     727           0 :                input%rdmftFunctional = 1
     728             :             CASE DEFAULT
     729           0 :                STOP 'Error: unknown RDMFT functional selected!'
     730             :          END SELECT
     731             :       END IF
     732             : 
     733             :       ! Read in optional q point mesh for spin spirals
     734             : 
     735          24 :       xPathA = '/fleurInput/calculationSetup/spinSpiralQPointMesh'
     736          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     737             : 
     738             :       !   IF ((noco%l_ss).AND.(numberNodes.EQ.0)) THEN
     739             :       !      call juDFT_error('Error: l_ss is true but no q point mesh set in xml input file!')
     740             :       !   END IF
     741             : 
     742             :       ! Read in optional E-Field parameters
     743             :      
     744          24 :       xPathA = '/fleurInput/calculationSetup/fields'
     745          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     746          24 :       field%b_field=0.0
     747          24 :       field%l_b_field=.FALSE.
     748          24 :       field%efield%sigma=0.0
     749             :       
     750             : 
     751          24 :       IF (numberNodes.EQ.1) THEN
     752           0 :          IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@b_field')>0) THEN
     753           0 :             field%b_field=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'//@b_field'))
     754           0 :             IF (ABS(field%b_field)>1.E-15) field%l_b_field=.true.
     755             :          ENDIF
     756           0 :          field%efield%zsigma = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zsigma'))
     757           0 :          field%efield%sig_b(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_1'))
     758           0 :          field%efield%sig_b(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sig_b_2'))
     759           0 :          field%efield%plot_charge = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_charge'))
     760           0 :          field%efield%plot_rho = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plot_rho'))
     761           0 :          field%efield%autocomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@autocomp'))
     762           0 :          field%efield%dirichlet = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dirichlet'))
     763           0 :          field%efield%l_eV = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eV'))
     764             : 
     765           0 :          numberNodes=xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/shape')
     766           0 :          ALLOCATE(field%efield%shapes(numberNodes))
     767           0 :          DO i=1,numberNodes
     768           0 :             WRITE(xPathB,"(a,a,i0,a)") TRIM(ADJUSTL(xpathA)),'/shape[',i,']'
     769           0 :             field%efield%shapes(i)=TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
     770             :          ENDDO
     771             :       ELSE
     772          24 :          ALLOCATE(field%efield%shapes(0))
     773             :       END IF
     774             : 
     775             :       ! Read in optional energy parameter limits
     776             : 
     777          24 :       xPathA = '/fleurInput/calculationSetup/energyParameterLimits'
     778          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     779             : 
     780          24 :       IF (numberNodes.EQ.1) THEN
     781          24 :          input%ellow = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ellow'))
     782          24 :          input%elup = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@elup'))
     783             :       END IF
     784             : 
     785             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     786             : !!! End of calculationSetup section
     787             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     788             : 
     789             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     790             : !!! Start of cell section
     791             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     792             : 
     793             :       ! Read in lattice parameters
     794             : 
     795          24 :       a1 = 0.0
     796          24 :       a2 = 0.0
     797          24 :       a3 = 0.0
     798          24 :       cell%z1 = 0.0
     799          24 :       dtild = 0.0
     800          24 :       input%film = .FALSE.
     801          24 :       latticeType = 'bulkLattice'
     802          24 :       latticeDef = 0
     803          24 :       symmetryDef = 0
     804          24 :       cell%latnam = 'any'
     805             : 
     806          24 :       numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/filmLattice')
     807             : 
     808          24 :       IF (numberNodes.EQ.1) THEN
     809           6 :          input%film = .TRUE.
     810           6 :          latticeType = 'filmLattice'
     811             :       END IF
     812             : 
     813          24 :       xPathA = '/fleurInput/cell/'//latticeType
     814          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     815             : 
     816          24 :       IF (numberNodes.EQ.1) THEN
     817          24 :          latticeScale = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@scale'))
     818          24 :          input%scaleCell = latticeScale
     819          24 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@latnam')))
     820          24 :          READ(valueString,*) cell%latnam
     821             : 
     822          24 :          IF(input%film) THEN
     823           6 :             cell%z1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dVac'))
     824           6 :             dtild = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dTilda'))
     825           6 :             vacuum%dvac = cell%z1
     826           6 :             a3(3) = dtild
     827          18 :             evac0Temp = eVac0Default_const
     828           6 :             xPathB = TRIM(ADJUSTL(xPathA))//'/vacuumEnergyParameters'
     829           6 :             numberNodes = xmlGetNumberOfNodes(xPathB)
     830           6 :             IF(numberNodes.GE.1) THEN
     831          12 :                DO i = 1, numberNodes
     832           6 :                   xPathC = ''
     833           6 :                   WRITE(xPathC,'(a,i0,a)') TRIM(ADJUSTL(xPathB))//'[',i,']'
     834           6 :                   numVac = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@vacuum'))
     835           6 :                   eParamUp = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinUp'))
     836           6 :                   eParamDown = eParamUp
     837           6 :                   IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathC))//'/@spinDown').GE.1) THEN
     838           6 :                      eParamDown = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathC))//'/@spinDown'))
     839             :                   END IF
     840           6 :                   evac0Temp(numVac,1) = eParamUp
     841           6 :                   IF(input%jspins.GT.1) evac0Temp(numVac,2) = eParamDown
     842          12 :                   IF(i == 1) THEN
     843           6 :                      evac0Temp(3-numVac,1) = eParamUp
     844           6 :                      IF(input%jspins.GT.1) evac0Temp(3-numVac,2) = eParamDown
     845             :                   END IF
     846             :                END DO
     847             :             END IF
     848             :          END IF
     849             : 
     850          24 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a1')
     851          24 :          IF (numberNodes.EQ.1) THEN
     852          19 :             latticeDef = 1
     853          19 :             input%scaleA1 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1/@scale'))
     854          19 :             a1(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a1'))
     855          19 :             a1(1) = a1(1) * input%scaleA1
     856          19 :             numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/a2')
     857          19 :             IF (numberNodes.EQ.1) THEN
     858           0 :                latticeDef = 2
     859           0 :                input%scaleA2 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2/@scale'))
     860           0 :                a2(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/a2'))
     861           0 :                a2(2) = a2(2) * input%scaleA2
     862             :             END IF
     863          19 :             IF(.NOT.input%film) THEN
     864          13 :                input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale'))
     865          13 :                a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c'))
     866          13 :                a3(3) = a3(3) * input%scaleC
     867             :             END IF
     868             :          END IF
     869             : 
     870          24 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/row-1')
     871          24 :          IF (numberNodes.EQ.1) THEN
     872           0 :             latticeDef = 3
     873           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-1')))
     874           0 :             a1(1) = evaluateFirst(valueString)
     875           0 :             a1(2) = evaluateFirst(valueString)
     876           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/row-2')))
     877           0 :             a2(1) = evaluateFirst(valueString)
     878           0 :             a2(2) = evaluateFirst(valueString)
     879           0 :             IF(.NOT.input%film) THEN
     880           0 :                input%scaleC = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c/@scale'))
     881           0 :                a3(3) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/c'))
     882           0 :                a3(3) = a3(3) * input%scaleC
     883             :             END IF
     884             :          END IF
     885             : 
     886          24 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix')
     887          24 :          IF (numberNodes.EQ.1) THEN
     888           5 :             latticeDef = 4
     889           5 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-1')))
     890           5 :             a1(1) = evaluateFirst(valueString)
     891           5 :             a1(2) = evaluateFirst(valueString)
     892           5 :             IF(.NOT.input%film) THEN
     893           5 :                a1(3) = evaluateFirst(valueString)
     894             :             END IF
     895           5 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-2')))
     896           5 :             a2(1) = evaluateFirst(valueString)
     897           5 :             a2(2) = evaluateFirst(valueString)
     898           5 :             IF(.NOT.input%film) THEN
     899           5 :                a2(3) = evaluateFirst(valueString)
     900             :             END IF
     901           5 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/bravaisMatrix/row-3')))
     902           5 :             IF(.NOT.input%film) THEN
     903           5 :                a3(1) = evaluateFirst(valueString)
     904           5 :                a3(2) = evaluateFirst(valueString)
     905           5 :                a3(3) = evaluateFirst(valueString)
     906             :             ELSE
     907             :                !WRITE(*,*) 'Note: For film calculations only the upper left 2x2 part of the Bravais matrix is considered.'
     908             :             END IF
     909             :          END IF
     910             :       END IF ! Note: Further ways to define lattices might be added later. (1D lattice,...)
     911             : 
     912             :       ! Construction of amat requires additional information about the lattice
     913             :       ! and is done later (scroll down)!
     914             : 
     915             :       ! Read in symmetry parameters
     916             : 
     917          24 :       sym%namgrp = 'any'
     918             : 
     919          24 :       xPathA = '/fleurInput/cell/symmetry'
     920          24 :       numberNodes = xmlGetNumberOfNodes('/fleurInput/cell/symmetry')
     921             : 
     922          24 :       IF (numberNodes.EQ.1) THEN
     923           4 :          sym%symSpecType = 2
     924           4 :          symmetryDef = 1
     925           4 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@spgrp')))
     926           4 :          READ(valueString,*) sym%namgrp
     927           4 :          sym%invs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@invs'))
     928           4 :          sym%zrfs = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@zrfs'))
     929           4 :          sym%invs2 = sym%invs.AND.sym%zrfs
     930             : 
     931           4 :          IF (sym%namgrp.EQ.'any ') THEN
     932           0 :             sym%nop = 48
     933             :             ! Read in sym.out file if sym%namgrp='any' set.
     934             :             CALL rw_symfile('r',94,'sym.out',48,cell%bmat,&
     935           0 :                &                        mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor)
     936           0 :             IF (ALLOCATED(sym%mrot)) THEN
     937           0 :                DEALLOCATE(sym%mrot)
     938             :             END IF
     939           0 :             ALLOCATE(sym%mrot(3,3,sym%nop))
     940           0 :             IF (ALLOCATED(sym%tau)) THEN
     941           0 :                DEALLOCATE(sym%tau)
     942             :             END IF
     943           0 :             ALLOCATE(sym%tau(3,sym%nop))
     944             : 
     945           0 :             DO k = 1, sym%nop
     946           0 :                DO i = 1, 3
     947           0 :                   DO j = 1, 3
     948           0 :                      sym%mrot(j,i,k) = mrotTemp(j,i,k)
     949             :                   END DO
     950           0 :                   sym%tau(i,k) = tauTemp(i,k)
     951             :                END DO
     952             :             END DO
     953             :          ELSE
     954           4 :             n2spg = 0
     955          84 :             DO i = 1, 20
     956          84 :                IF (sym%namgrp.EQ.nammap(i)) n2spg = i
     957             :             END DO
     958           4 :             IF (n2spg == 0 ) THEN
     959           0 :                WRITE (*,*) 'Spacegroup ',sym%namgrp,' not known! Choose one of:'
     960           0 :                WRITE (*,'(20(a4,1x))') (nammap(i),i=1,20)
     961           0 :                CALL juDFT_error("Could not determine spacegroup!", calledby = "r_inpXML")
     962             :             END IF
     963           4 :             IF ((n2spg.GE.13).AND.(n2spg.LE.17)) THEN
     964           0 :                IF (.NOT.((cell%latnam.EQ.'hx3').OR.(cell%latnam.EQ.'hex'))) THEN
     965           0 :                   CALL juDFT_error("Use only hex or hx3 with p3, p3m1, p31m, p6 or p6m!", calledby ="r_inpXML")
     966             :                END IF
     967             :             END IF
     968           4 :             sym%nop = ord2(n2spg)
     969           4 :             IF (sym%invs) THEN
     970           4 :                sym%nop = 2*sym%nop
     971           4 :                IF (sym%zrfs.AND.(.NOT.l_c2(n2spg))) sym%nop = 2*sym%nop
     972             :             ELSE
     973           0 :                IF (sym%zrfs) sym%nop = 2*sym%nop
     974             :             END IF
     975           4 :             IF (ALLOCATED(sym%mrot)) THEN
     976           0 :                DEALLOCATE(sym%mrot)
     977             :             END IF
     978           4 :             ALLOCATE(sym%mrot(3,3,sym%nop))
     979           4 :             IF (ALLOCATED(sym%tau)) THEN
     980           0 :                DEALLOCATE(sym%tau)
     981             :             END IF
     982           4 :             ALLOCATE(sym%tau(3,sym%nop))
     983             :             CALL spg2set(sym%nop,sym%zrfs,sym%invs,sym%namgrp,cell%latnam,&
     984           4 :                &                     sym%mrot,sym%tau,sym%nop2,sym%symor)
     985             :          END IF
     986             :       END IF
     987             : 
     988          24 :       xPathA = '/fleurInput/cell/symmetryFile'
     989          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
     990             : 
     991          24 :       IF (numberNodes.EQ.1) THEN
     992          20 :          symmetryDef = 2
     993          20 :          sym%symSpecType = 1
     994          20 :          sym%nop = 48
     995          20 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@filename')))
     996             : 
     997             :          CALL rw_symfile('r',94,TRIM(ADJUSTL(valueString)),48,cell%bmat,&
     998          20 :             &                     mrotTemp,tauTemp,sym%nop,sym%nop2,sym%symor)
     999             : 
    1000          20 :          IF (ALLOCATED(sym%mrot)) THEN
    1001           0 :             DEALLOCATE(sym%mrot)
    1002             :          END IF
    1003          20 :          ALLOCATE(sym%mrot(3,3,sym%nop))
    1004          20 :          IF (ALLOCATED(sym%tau)) THEN
    1005           0 :             DEALLOCATE(sym%tau)
    1006             :          END IF
    1007          20 :          ALLOCATE(sym%tau(3,sym%nop))
    1008             : 
    1009          20 :          sym%invs = .FALSE.
    1010          20 :          sym%zrfs = .FALSE.
    1011             : 
    1012         247 :          DO k = 1, sym%nop
    1013         227 :             absSum = 0
    1014         908 :             DO i = 1, 3
    1015        4767 :                DO j = 1, 3
    1016        2043 :                   sym%mrot(j,i,k) = mrotTemp(j,i,k)
    1017        2724 :                   absSum = absSum + ABS(sym%mrot(j,i,k))
    1018             :                END DO
    1019         908 :                sym%tau(i,k) = tauTemp(i,k)
    1020             :             END DO
    1021         247 :             IF (absSum.EQ.3) THEN
    1022         151 :                IF (ALL(sym%tau(:,k).EQ.0.0)) THEN
    1023         149 :                   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.
    1024         149 :                   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.
    1025             :                END IF
    1026             :             END IF
    1027             :          END DO
    1028             : 
    1029          20 :          sym%invs2 = sym%invs.AND.sym%zrfs
    1030             :       END IF
    1031             : 
    1032          24 :       xPathA = '/fleurInput/cell/symmetryOperations'
    1033          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
    1034             : 
    1035          24 :       IF (numberNodes.EQ.1) THEN
    1036           0 :          sym%symSpecType = 3
    1037           0 :          symmetryDef = 3
    1038             : 
    1039           0 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/symOp')
    1040           0 :          sym%nop = numberNodes
    1041             : 
    1042           0 :          IF (ALLOCATED(sym%mrot)) DEALLOCATE(sym%mrot)
    1043           0 :          ALLOCATE(sym%mrot(3,3,sym%nop))
    1044           0 :          IF (ALLOCATED(sym%tau)) DEALLOCATE(sym%tau)
    1045           0 :          ALLOCATE(sym%tau(3,sym%nop))
    1046           0 :          sym%symor = .TRUE.
    1047           0 :          DO i = 1, sym%nop
    1048           0 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-1'
    1049           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
    1050           0 :             READ(valueString,*) sym%mrot(1,1,i), sym%mrot(1,2,i), sym%mrot(1,3,i), sym%tau(1,i)
    1051             : 
    1052           0 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-2'
    1053           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
    1054           0 :             READ(valueString,*) sym%mrot(2,1,i), sym%mrot(2,2,i), sym%mrot(2,3,i), sym%tau(2,i)
    1055             : 
    1056           0 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA))//'/symOp[',i,']/row-3'
    1057           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))))
    1058           0 :             READ(valueString,*) sym%mrot(3,1,i), sym%mrot(3,2,i), sym%mrot(3,3,i), sym%tau(3,i)
    1059             : 
    1060           0 :             IF ((sym%tau(1,i)**2 + sym%tau(2,i)**2 + sym%tau(3,i)**2).GT.1.e-8) THEN
    1061           0 :                sym%symor = .FALSE.
    1062             :             END IF
    1063           0 :             DO j = 1,3
    1064           0 :                IF (ABS(sym%tau(j,i)-0.33333) < 0.00001) THEN
    1065           0 :                   sym%tau(j,i) = 1./3.
    1066             :                ENDIF
    1067           0 :                IF (ABS(sym%tau(j,i)+0.33333) < 0.00001) THEN
    1068           0 :                   sym%tau(j,i) = -1./3.
    1069             :                ENDIF
    1070           0 :                IF (ABS(sym%tau(j,i)-0.66667) < 0.00001) THEN
    1071           0 :                   sym%tau(j,i) = 2./3.
    1072             :                ENDIF
    1073           0 :                IF (ABS(sym%tau(j,i)+0.66667) < 0.00001) THEN
    1074           0 :                   sym%tau(j,i) = -2./3.
    1075             :                ENDIF
    1076             :             ENDDO
    1077             :          END DO
    1078             :       END IF
    1079             : 
    1080             :       ! Calculate missing symmetry and cell properties and check consistency of parameters.
    1081             : 
    1082             :       ! Construction of amat
    1083          19 :       SELECT CASE (latticeDef)
    1084             :       CASE (1)
    1085          19 :          IF (cell%latnam.EQ.'squ') THEN
    1086          13 :             a2(2) = a1(1)
    1087           6 :          ELSE IF (cell%latnam.EQ.'c-b') THEN
    1088           0 :             aTemp = a1(1)
    1089           0 :             a1(1) = aTemp*0.5*SQRT(2.0)
    1090           0 :             a1(2) = -aTemp*0.5
    1091           0 :             a2(1) = aTemp*0.5*SQRT(2.0)
    1092           0 :             a2(2) = aTemp*0.5
    1093           6 :          ELSE IF (cell%latnam.EQ.'hex') THEN
    1094           2 :             aTemp = 0.5*a1(1)
    1095           2 :             a1(1) = aTemp*SQRT(3.0)
    1096           2 :             a1(2) = -aTemp
    1097           2 :             a2(1) = a1(1)
    1098           2 :             a2(2) = aTemp
    1099           4 :          ELSE IF (cell%latnam.EQ.'hx3') THEN
    1100           4 :             aTemp = 0.5*a1(1)
    1101           4 :             a1(1) = aTemp
    1102           4 :             a1(2) = -aTemp*SQRT(3.0)
    1103           4 :             a2(1) = a1(1)
    1104           4 :             a2(2) = -a1(2)
    1105           0 :          ELSE IF (cell%latnam.EQ.'fcc') THEN
    1106           0 :             aTemp = a1(1)
    1107           0 :             a1(1) =       0.0 ; a1(2) = 0.5*aTemp ; a1(3) = 0.5*aTemp
    1108           0 :             a2(1) = 0.5*aTemp ; a2(2) =       0.0 ; a2(3) = 0.5*aTemp
    1109           0 :             a3(1) = 0.5*aTemp ; a3(2) = 0.5*aTemp ; a3(3) =       0.0
    1110           0 :          ELSE IF (cell%latnam.EQ.'bcc') THEN
    1111           0 :             aTemp = a1(1)
    1112           0 :             a1(1) =-0.5*aTemp ; a1(2) = 0.5*aTemp ; a1(3) = 0.5*aTemp
    1113           0 :             a2(1) = 0.5*aTemp ; a2(2) =-0.5*aTemp ; a2(3) = 0.5*aTemp
    1114           0 :             a3(1) = 0.5*aTemp ; a3(2) = 0.5*aTemp ; a3(3) =-0.5*aTemp
    1115             :          ELSE
    1116           0 :             CALL juDFT_error("latnam is incompatible to parametrization of lattice (1)", calledby ="r_inpXML")
    1117             :          END IF
    1118             :       CASE (2)
    1119           0 :          IF ((cell%latnam.EQ.'c-r').OR.(cell%latnam.EQ.'p-r')) THEN
    1120           0 :             IF (cell%latnam.EQ.'c-r') THEN
    1121           0 :                a1(2) = -a2(2)
    1122           0 :                a2(1) =  a1(1)
    1123             :             END IF
    1124             :          ELSE
    1125           0 :             CALL juDFT_error("latnam is incompatible to parametrization of lattice (2)", calledby ="r_inpXML")
    1126             :          END IF
    1127             :       CASE (3)
    1128           5 :          IF (.NOT.(cell%latnam.EQ.'obl')) THEN
    1129           0 :             CALL juDFT_error("latnam is incompatible to parametrization of lattice (3)", calledby ="r_inpXML")
    1130             :          END IF
    1131             :       CASE (4)
    1132           5 :          IF (.NOT.(cell%latnam.EQ.'any')) THEN
    1133           0 :             CALL juDFT_error("latnam is incompatible to parametrization of lattice (4)", calledby ="r_inpXML")
    1134             :          END IF
    1135             :       CASE DEFAULT
    1136          24 :          CALL juDFT_error("Illegal lattice definition", calledby ="r_inpXML")
    1137             :       END SELECT
    1138             : 
    1139          24 :       IF (latticeScale.EQ.0.0) latticeScale = 1.0
    1140          24 :       IF (.NOT.input%film) vacuum%dvac = a3(3)
    1141          24 :       vacuum%dvac = latticeScale*vacuum%dvac
    1142          24 :       dtild = latticeScale*dtild
    1143             : 
    1144          96 :       cell%amat(:,1) = a1(:) * latticeScale
    1145          96 :       cell%amat(:,2) = a2(:) * latticeScale
    1146          96 :       cell%amat(:,3) = a3(:) * latticeScale
    1147             : 
    1148          24 :       CALL inv3(cell%amat,cell%bmat,cell%omtil)
    1149          96 :       cell%bmat(:,:) = tpi_const*cell%bmat(:,:)
    1150          24 :       cell%omtil = ABS(cell%omtil)
    1151             : 
    1152          24 :       IF (input%film.AND..NOT.oneD%odd%d1) THEN
    1153           6 :          cell%vol = (cell%omtil/cell%amat(3,3))*vacuum%dvac
    1154           6 :          cell%area = cell%amat(1,1)*cell%amat(2,2)-cell%amat(1,2)*cell%amat(2,1)
    1155             :          !-odim
    1156          18 :       ELSE IF (oneD%odd%d1) THEN
    1157           0 :          cell%area = tpi_const*cell%amat(3,3)
    1158           0 :          cell%vol = pi_const*(vacuum%dvac**2)*cell%amat(3,3)/4.0
    1159             :          !+odim
    1160             :       ELSE
    1161          18 :          cell%vol = cell%omtil
    1162          18 :          cell%area = cell%amat(1,1)*cell%amat(2,2)-cell%amat(1,2)*cell%amat(2,1)
    1163          18 :          IF (cell%area < 1.0e-7) THEN
    1164           5 :             IF (cell%latnam.EQ.'any') THEN
    1165           5 :                cell%area = 1.
    1166             :             ELSE
    1167           0 :                CALL juDFT_error("area = 0",calledby ="r_inpXML")
    1168             :             END IF
    1169             :          END IF
    1170             :       END IF
    1171             : 
    1172             :       ! Construction of missing symmetry information
    1173          24 :       IF ((symmetryDef.EQ.2).OR.(symmetryDef.EQ.3)) THEN
    1174          20 :          nop48 = 48
    1175          20 :          ALLOCATE (invOps(sym%nop),multtab(sym%nop,sym%nop),optype(nop48))
    1176             :          CALL check_close(sym%nop,sym%mrot,sym%tau,&
    1177          20 :             &                      multtab,invOps,optype)
    1178             : 
    1179             :          CALL symproperties(nop48,optype,input%film,sym%nop,multtab,cell%amat,&
    1180             :             &                        sym%symor,sym%mrot,sym%tau,&
    1181          20 :             &                        invSym,sym%invs,sym%zrfs,sym%invs2,sym%nop,sym%nop2)
    1182          20 :          DEALLOCATE(invOps,multtab,optype)
    1183          20 :          IF (.NOT.input%film) sym%nop2=sym%nop
    1184          20 :          IF (input%film) THEN
    1185          58 :             DO n = 1, sym%nop
    1186         214 :                DO i = 1, 3
    1187         208 :                   IF (ABS(sym%tau(i,n)) > 0.00001) THEN
    1188           0 :                      CALL juDFT_error("nonsymmorphic symmetries not yet implemented for films!",calledby ="r_inpXML")
    1189             :                   ENDIF
    1190             :                END DO
    1191             :             END DO
    1192             :          END IF
    1193             :       END IF
    1194          24 :       sym%invs2 = sym%invs.AND.sym%zrfs
    1195             : 
    1196          24 :       ALLOCATE (sym%invarop(atoms%nat,sym%nop),sym%invarind(atoms%nat))
    1197          24 :       ALLOCATE (sym%multab(sym%nop,sym%nop),sym%invtab(sym%nop))
    1198          24 :       ALLOCATE (sym%invsatnr(atoms%nat),sym%d_wgn(-3:3,-3:3,3,sym%nop))
    1199             : 
    1200             :       !some settings for film calculations
    1201          24 :       vacuum%nmzd = 250
    1202          24 :       vacuum%nmzxyd = 100
    1203          24 :       vacuum%nvac = 2
    1204          24 :       IF (sym%zrfs.OR.sym%invs) vacuum%nvac = 1
    1205          24 :       IF (oneD%odd%d1) vacuum%nvac = 1
    1206          24 :       cell%z1 = vacuum%dvac/2
    1207          24 :       vacuum%nmz = vacuum%nmzd
    1208          24 :       vacuum%delz = 25.0/vacuum%nmz
    1209          24 :       IF (oneD%odd%d1) vacuum%delz = 20.0 / vacuum%nmz
    1210             :       IF (vacuum%nmz.GT.vacuum%nmzd) CALL juDFT_error("nmzd",calledby ="inped")
    1211          24 :       vacuum%nmzxy = vacuum%nmzxyd
    1212             :       IF (vacuum%nmzxy.GT.vacuum%nmzxyd) CALL juDFT_error("nmzxyd",calledby ="inped")
    1213             : 
    1214             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1215             : !!! End of cell section
    1216             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1217             : 
    1218             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1219             : !!! Start of XC functional section
    1220             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1221             : 
    1222             :       !Read in libxc parameters if present
    1223          24 :       xPathA = '/fleurInput/xcFunctional/LibXCID'
    1224          24 :       xPathB = '/fleurInput/xcFunctional/LibXCName'
    1225             : 
    1226          24 :       IF(xmlGetNumberOfNodes(xPathA) == 1 .AND. xmlGetNumberOfNodes(xPathB) == 1) THEN
    1227           0 :          CALL judft_error("LibXC is given both by Name and ID and is therefore overdetermined", calledby="r_inpXML")
    1228             :       ENDIF
    1229             : 
    1230             : 
    1231             :       ! LibXCID 
    1232          24 :       IF (xmlGetNumberOfNodes(xPathA) == 1) THEN
    1233             : #ifdef CPP_LIBXC
    1234           1 :          vxc_id_x=evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@exchange'))
    1235           1 :          vxc_id_c=evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@correlation'))
    1236             : 
    1237           1 :          IF(xmlGetNumberOfNodes(TRIM(xPathA) // '/@etot_exchange') == 1) THEN
    1238           0 :             exc_id_x = evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@etot_exchange'))
    1239             :          ELSE
    1240           1 :             exc_id_x = vxc_id_x
    1241             :          ENDIF
    1242             :          
    1243           1 :          IF(xmlGetNumberOfNodes(TRIM(xPathA) // '/@exc_correlation') == 1) THEN
    1244           0 :             exc_id_c = evaluateFirstOnly(xmlGetAttributeValue(xPathA // '/@exc_correlation'))
    1245             :          ELSE
    1246           1 :             exc_id_c = vxc_id_c
    1247             :          ENDIF
    1248             : #else
    1249             :          CALL judft_error("To use libxc functionals you have to compile with libXC support")
    1250             : #endif
    1251             :       ! LibXCName 
    1252          23 :       ELSEIF (xmlGetNumberOfNodes(TRIM(xPathB)) == 1) THEN
    1253             : #ifdef CPP_LIBXC
    1254           1 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@exchange')))
    1255           1 :          vxc_id_x =  xc_f03_functional_get_number(TRIM(valueString))
    1256             : 
    1257           1 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@correlation')))
    1258           1 :          vxc_id_c =  xc_f03_functional_get_number(TRIM(valueString))
    1259             :          
    1260           1 :          IF(xmlGetNumberOfNodes(TRIM(xPathB) // '/@etot_exchange') == 1) THEN
    1261           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@etot_exchange')))
    1262           0 :             exc_id_x =  xc_f03_functional_get_number(TRIM(valueString))
    1263             :          ELSE
    1264           1 :             exc_id_x = vxc_id_x
    1265             :          ENDIF
    1266             :          
    1267           1 :          IF(xmlGetNumberOfNodes(TRIM(xPathB) // '/@etot_correlation') == 1) THEN
    1268           0 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(xPathB) // '/@etot_correlation')))
    1269           0 :             exc_id_c =  xc_f03_functional_get_number(TRIM(valueString))
    1270             :          ELSE
    1271           1 :             exc_id_c = vxc_id_c
    1272             :          ENDIF
    1273             : #else
    1274             :          CALL judft_error("To use libxc functionals you have to compile with libXC support")
    1275             : #endif
    1276             :       ELSE
    1277          22 :          vxc_id_x=0; vxc_id_c=0;
    1278          22 :          exc_id_x=0; exc_id_c=0;
    1279             :       ENDIF
    1280             : 
    1281             :       ! Read in xc functional parameters
    1282          24 :       valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL('/fleurInput/xcFunctional/@name')))))
    1283          24 :       namex(1:4) = valueString(1:4)
    1284          24 :       l_relcor = evaluateFirstBoolOnly(xmlGetAttributeValue('/fleurInput/xcFunctional/@relativisticCorrections'))
    1285             : 
    1286          24 :       relcor = 'non-relativi'
    1287          24 :       IF (l_relcor) THEN
    1288           0 :          relcor = 'relativistic'
    1289             :       END IF
    1290             : 
    1291             :       !now initialize the xcpot variable
    1292          24 :       CALL setXCParameters(atoms,valueString,l_relcor,input%jspins,vxc_id_x,vxc_id_c,exc_id_x, exc_id_c, xcpot)
    1293             : 
    1294          24 :       xPathA = '/fleurInput/calculationSetup/cutoffs/@GmaxXC'
    1295          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
    1296          24 :       xcpot%gmaxxc = stars%gmax
    1297          24 :       IF(numberNodes.EQ.1) THEN
    1298          24 :          xcpot%gmaxxc = evaluateFirstOnly(xmlGetAttributeValue(xPathA))
    1299             :       END IF
    1300          24 :       hybrid%l_hybrid=xcpot%is_hybrid()
    1301             : 
    1302          24 :       ALLOCATE(hybrid%lcutm1(atoms%ntype),hybrid%lcutwf(atoms%ntype),hybrid%select1(4,atoms%ntype))
    1303             : 
    1304          24 :       obsolete%lwb=.FALSE.
    1305          24 :       IF (xcpot%needs_grad()) THEN
    1306          20 :          obsolete%ndvgrd=6
    1307          20 :          obsolete%chng=-0.1e-11
    1308             :       END IF
    1309             : 
    1310          24 :       IF (xcpot%needs_grad()) THEN
    1311          20 :          obsolete%ndvgrd = MAX(obsolete%ndvgrd,3)
    1312             :       END IF
    1313             : 
    1314          24 :       hybrid%gcutm1 = input%rkmax - 0.5
    1315          24 :       hybrid%tolerance1 = 1.0e-4
    1316          24 :       hybrid%ewaldlambda = 3
    1317          24 :       hybrid%lexp = 16
    1318          24 :       hybrid%bands1 = DIMENSION%neigd
    1319             : 
    1320          24 :       numberNodes = xmlGetNumberOfNodes('/fleurInput/calculationSetup/prodBasis')
    1321          24 :       IF (numberNodes==0) THEN
    1322          24 :          IF (hybrid%l_hybrid) CALL judft_error("Mixed product basis input missing in inp.xml")
    1323             :       ELSE
    1324           0 :          hybrid%gcutm1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@gcutm'))
    1325           0 :          hybrid%tolerance1=evaluateFirstOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@tolerance'))
    1326           0 :          hybrid%ewaldlambda=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@ewaldlambda'))
    1327           0 :          hybrid%lexp=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@lexp'))
    1328           0 :          hybrid%bands1=evaluateFirstIntOnly(xmlGetAttributeValue('/fleurInput/calculationSetup/prodBasis/@bands'))
    1329             :       ENDIF
    1330             : 
    1331             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1332             : !!! End of XC functional section
    1333             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1334             : 
    1335             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1336             : !!! Start of species section
    1337             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1338             : 
    1339          24 :       ALLOCATE (speciesNLO(numSpecies))
    1340          24 :       ALLOCATE(atoms%speciesName(numSpecies))
    1341             : 
    1342          24 :       atoms%numStatesProvided = 0
    1343          24 :       atoms%lapw_l(:) = -1
    1344          24 :       atoms%n_u = 0
    1345             : 
    1346          24 :       DEALLOCATE(noel)
    1347          24 :       ALLOCATE(noel(atoms%ntype))
    1348             : 
    1349          57 :       DO iSpecies = 1, numSpecies
    1350             :          ! Attributes of species
    1351          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']'
    1352          33 :          atoms%speciesName(iSpecies) = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@name')))
    1353          33 :          atomicNumber = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomicNumber'))
    1354          33 :          coreStates = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@coreStates'))
    1355          33 :          magMom = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@magMom'))
    1356          33 :          flipSpin = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@flipSpin'))
    1357             : 
    1358             :          ! Attributes of mtSphere element of species
    1359          33 :          radius = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@radius'))
    1360          33 :          gridPoints = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@gridPoints'))
    1361          33 :          logIncrement = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/mtSphere/@logIncrement'))
    1362             : 
    1363             :          ! Attributes of atomicCutoffs element of species
    1364          33 :          lmax = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmax'))
    1365          33 :          lnonsphr = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lnonsphr'))
    1366          33 :          lmaxAPW = -1
    1367          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW')
    1368          33 :          IF (numberNodes.EQ.1) THEN
    1369           0 :             lmaxAPW = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/atomicCutoffs/@lmaxAPW'))
    1370             :          END IF
    1371             : 
    1372          33 :          numU = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/ldaU')
    1373          33 :         IF (numU.GT.4) CALL juDFT_error("Too many U parameters provided for a certain species (maximum is 4).",calledby ="r_inpXML")
    1374         165 :          ldau_l = -1
    1375          33 :          ldau_u = 0.0
    1376          33 :          ldau_j = 0.0
    1377          33 :          l_amf = .FALSE.
    1378          45 :          DO i = 1, numU
    1379          12 :             WRITE(xPathB,*) i
    1380          12 :             ldau_l(i) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l'))
    1381          12 :             ldau_u(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@U'))
    1382          12 :             ldau_j(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@J'))
    1383          12 :             ldau_phi(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@phi'))
    1384          12 :             ldau_theta(i) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@theta'))
    1385          45 :             l_amf(i) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/ldaU['//TRIM(ADJUSTL(xPathB))//']/@l_amf'))
    1386             :          END DO
    1387             : 
    1388          33 :          speciesNLO(iSpecies) = 0
    1389          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo'
    1390          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1391          42 :          DO iLO = 1, numberNodes
    1392           9 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l'
    1393           9 :             WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n'
    1394           9 :             lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))
    1395           9 :             nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC)))
    1396           9 :             CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount)
    1397           9 :             CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount)
    1398           9 :             IF(lNumCount.NE.nNumCount) THEN
    1399           0 :                CALL judft_error('Error in LO input: l quantum number count does not equal n quantum number count')
    1400             :             END IF
    1401           9 :             speciesNLO(iSpecies) = speciesNLO(iSpecies) + lNumCount
    1402          42 :             DEALLOCATE (lNumbers, nNumbers)
    1403             :          END DO
    1404             :          ! Special switches for species
    1405          33 :          vcaspecies=0.0
    1406          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/special'
    1407          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1408          33 :          IF (numberNodes==1) THEN
    1409           0 :             vcaSpecies   = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vca_charge'))))
    1410             :          ENDIF
    1411             :        
    1412         138 :          DO iType = 1, atoms%ntype
    1413          81 :             WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species'
    1414          81 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
    1415         114 :             IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN
    1416          55 :                atoms%nz(iType) = atomicNumber
    1417          55 :                atoms%zatom(iType) = atoms%nz(iType)
    1418          55 :                IF (atoms%nz(iType).EQ.0) THEN
    1419           0 :                   WRITE(*,*) 'Note: Replacing atomic number 0 by 1.0e-10 on atom type ', iType
    1420           0 :                   atoms%zatom(iType) = 1.0e-10
    1421             :                END IF
    1422          55 :                atoms%zatom(iType)=atoms%zatom(iType)+vcaspecies
    1423          55 :                noel(iType) = namat_const(atoms%nz(iType))
    1424          55 :                atoms%rmt(iType) = radius
    1425          55 :                atoms%jri(iType) = gridPoints
    1426          55 :                atoms%dx(iType) = logIncrement
    1427          55 :                atoms%lmax(iType) = lmax
    1428          55 :                atoms%nlo(iType) = speciesNLO(iSpecies)
    1429          55 :                atoms%ncst(iType) = coreStates
    1430          55 :                atoms%lnonsph(iType) = lnonsphr
    1431          55 :                atoms%lapw_l(iType) = lmaxAPW
    1432          55 :                IF (flipSpin) THEN
    1433          55 :                   atoms%nflip(iType) = -1
    1434             :                ELSE
    1435           0 :                   atoms%nflip(iType) = 0
    1436             :                ENDIF
    1437          55 :                atoms%bmu(iType) = magMom
    1438          67 :                DO i = 1, numU
    1439          12 :                   atoms%n_u = atoms%n_u + 1
    1440          12 :                   atoms%lda_u(atoms%n_u)%l = ldau_l(i)
    1441          12 :                   atoms%lda_u(atoms%n_u)%u = ldau_u(i)
    1442          12 :                   atoms%lda_u(atoms%n_u)%j = ldau_j(i)
    1443          12 :                   atoms%lda_u(atoms%n_u)%phi = ldau_phi(i)
    1444          12 :                   atoms%lda_u(atoms%n_u)%theta = ldau_theta(i)
    1445          12 :                   atoms%lda_u(atoms%n_u)%l_amf = l_amf(i)
    1446          67 :                   atoms%lda_u(atoms%n_u)%atomType = iType
    1447             :                END DO
    1448          55 :                atomTypeSpecies(iType) = iSpecies
    1449          55 :                IF(speciesRepAtomType(iSpecies).EQ.-1) speciesRepAtomType(iSpecies) = iType
    1450             :             END IF
    1451             :          END DO
    1452             :       END DO
    1453             : 
    1454          24 :       atoms%lmaxd = MAXVAL(atoms%lmax(:))
    1455          24 :       atoms%llod  = 0
    1456          24 :       atoms%nlod = 0
    1457          79 :       DO iType = 1, atoms%ntype
    1458          79 :          atoms%nlod = MAX(atoms%nlod,atoms%nlo(iType))
    1459             :       END DO
    1460          24 :       atoms%nlod = MAX(atoms%nlod,2) ! for chkmt
    1461          24 :       ALLOCATE(atoms%llo(atoms%nlod,atoms%ntype)); atoms%llo=-1
    1462          24 :       ALLOCATE(atoms%ulo_der(atoms%nlod,atoms%ntype))
    1463          24 :       ALLOCATE(atoms%l_dulo(atoms%nlod,atoms%ntype)) ! For what is this?
    1464             : 
    1465          24 :       DIMENSION%nstd = 29
    1466             : 
    1467          24 :       ALLOCATE(atoms%coreStateOccs(DIMENSION%nstd,2,atoms%ntype)); atoms%coreStateOccs=0.0
    1468          24 :       ALLOCATE(atoms%coreStateNprnc(DIMENSION%nstd,atoms%ntype))
    1469          24 :       ALLOCATE(atoms%coreStateKappa(DIMENSION%nstd,atoms%ntype))
    1470             : 
    1471          24 :       CALL enpara%init(atoms,input%jspins)
    1472          72 :       enpara%evac0(:,:) = evac0Temp(:,:)
    1473             : 
    1474          57 :       DO iSpecies = 1, numSpecies
    1475          33 :          ALLOCATE(speciesLLO(speciesNLO(iSpecies)))
    1476          33 :          ALLOCATE(speciesLOeParams(speciesNLO(iSpecies)))
    1477          33 :          ALLOCATE(speciesLOEDeriv(speciesNLO(iSpecies)))
    1478             : 
    1479             :          ! Attributes of energyParameters element of species
    1480          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']'
    1481          33 :          speciesEParams(0) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@s'))
    1482          33 :          speciesEParams(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@p'))
    1483          33 :          speciesEParams(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@d'))
    1484          33 :          speciesEParams(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/energyParameters/@f'))
    1485             : 
    1486             :          ! Parameters for hybrid functionals
    1487          33 :          IF (hybrid%l_hybrid) THEN
    1488           0 :             WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/prodBasis'
    1489           0 :             numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1490           0 :             IF (numberNodes.NE.1) CALL judft_error("Parameters for mixed basis are missing for some specified")
    1491           0 :             lcutm =evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lcutm'))
    1492           0 :             lcutwf=evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lcutwf'))
    1493           0 :             xPathA=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@select')
    1494           0 :             hybSelect(1) = NINT(evaluateFirst(xPathA))
    1495           0 :             hybSelect(2) = NINT(evaluateFirst(xPathA))
    1496           0 :             hybSelect(3) = NINT(evaluateFirst(xPathA))
    1497           0 :             hybSelect(4) = NINT(evaluateFirst(xPathA))
    1498             :          ENDIF
    1499             : 
    1500             :          ! Special switches for species
    1501          33 :          ldaspecies=.FALSE.
    1502          33 :          socscalespecies=1.0
    1503          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/special'
    1504          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1505          33 :          IF (numberNodes==1) THEN
    1506           0 :             ldaSpecies = evaluateFirstBoolOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lda'))))
    1507           0 :             socscaleSpecies   = evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socscale'))))
    1508           0 :             IF (xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/@b_field_mt')>0) THEN
    1509           0 :                b_field_mtSpecies=evaluateFirstOnly(TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@b_field_mt'))))
    1510           0 :                IF (ABS(b_field_mtSpecies)>1E-15) field%l_b_field=.TRUE.
    1511             :             ENDIF
    1512             :          ENDIF
    1513             :          ! Explicitely provided core configurations
    1514             : 
    1515          33 :          coreConfigPresent = .FALSE.
    1516          33 :          providedCoreStates = 0
    1517          33 :          providedStates = 0
    1518          33 :          coreStateOccs = 0.0
    1519          33 :          speciesXMLElectronStates = noState_const
    1520         990 :          speciesXMLCoreOccs = -1.0
    1521          33 :          speciesXMLPrintCoreStates = .FALSE.
    1522          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/electronConfig'
    1523          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1524          33 :          IF (numberNodes.EQ.1) THEN
    1525           0 :             coreConfigPresent = .TRUE.
    1526           0 :             valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/coreConfig')
    1527           0 :             token = popFirstStringToken(valueString)
    1528           0 :             DO WHILE (token.NE.' ')
    1529           0 :                IF (token(1:1).EQ.'[') THEN
    1530           0 :                   DO i = 1, 6
    1531           0 :                      IF (TRIM(ADJUSTL(token)).EQ.nobleGasConfigList_const(i)) THEN
    1532           0 :                         IF (providedCoreStates+nobleGasNumStatesList_const(i).GT.29) THEN
    1533           0 :                            CALL judft_error('Error: Too many core states provided in xml input file!')
    1534             :                         END IF
    1535           0 :                         DO j = providedCoreStates+1, providedCoreStates+nobleGasNumStatesList_const(i)
    1536           0 :                            coreStateOccs(j-providedCoreStates,:) = coreStateNumElecsList_const(j)
    1537           0 :                            coreStateNprnc(j-providedCoreStates) = coreStateNprncList_const(j)
    1538           0 :                            coreStateKappa(j-providedCoreStates) = coreStateKappaList_const(j)
    1539           0 :                            speciesXMLElectronStates(j) = coreState_const
    1540             :                         END DO
    1541           0 :                         providedCoreStates = providedCoreStates + nobleGasNumStatesList_const(i)
    1542             :                      END IF
    1543             :                   END DO
    1544             :                ELSE
    1545           0 :                   DO i = 1, 29
    1546           0 :                      IF (TRIM(ADJUSTL(token)).EQ.coreStateList_const(i)) THEN
    1547           0 :                         providedCoreStates = providedCoreStates + 1
    1548           0 :                         IF (providedCoreStates.GT.29) THEN
    1549           0 :                            CALL judft_error('Error: Too many core states provided in xml input file!')
    1550             :                         END IF
    1551           0 :                         coreStateOccs(providedCoreStates,:) = coreStateNumElecsList_const(i)
    1552           0 :                         coreStateNprnc(providedCoreStates) = coreStateNprncList_const(i)
    1553           0 :                         coreStateKappa(providedCoreStates) = coreStateKappaList_const(i)
    1554           0 :                         speciesXMLElectronStates(i) = coreState_const
    1555             :                      END IF
    1556             :                   END DO
    1557             :                END IF
    1558           0 :                token = popFirstStringToken(valueString)
    1559             :             END DO
    1560           0 :             numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/valenceConfig')
    1561           0 :             providedStates = providedCoreStates
    1562           0 :             IF(numberNodes.EQ.1) THEN
    1563           0 :                valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/valenceConfig')
    1564           0 :                token = popFirstStringToken(valueString)
    1565           0 :                DO WHILE (token.NE.' ')
    1566           0 :                   DO i = 1, 29
    1567           0 :                      IF (TRIM(ADJUSTL(token)).EQ.coreStateList_const(i)) THEN
    1568           0 :                         providedStates = providedStates + 1
    1569           0 :                         IF (providedStates.GT.29) THEN
    1570           0 :                            CALL judft_error('Error: Too many valence states provided in xml input file!')
    1571             :                         END IF
    1572           0 :                         coreStateOccs(providedStates,:) = coreStateNumElecsList_const(i)
    1573           0 :                         coreStateNprnc(providedStates) = coreStateNprncList_const(i)
    1574           0 :                         coreStateKappa(providedStates) = coreStateKappaList_const(i)
    1575           0 :                         speciesXMLElectronStates(i) = valenceState_const
    1576             :                      END IF
    1577             :                   END DO
    1578           0 :                   token = popFirstStringToken(valueString)
    1579             :                END DO
    1580             :             END IF
    1581             :          END IF
    1582             : 
    1583             :          ! Explicitely provided core occupations
    1584             : 
    1585          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/electronConfig/stateOccupation'
    1586          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1587          33 :          IF (numberNodes.GE.1) THEN
    1588           0 :             IF (.NOT.coreConfigPresent) THEN
    1589           0 :                WRITE(*,*) 'Note: This just has to be implemented:'
    1590           0 :                CALL judft_error('Error: Core occupation given while core config not set!')
    1591             :             END IF
    1592           0 :             DO i = 1, numberNodes
    1593           0 :                WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',i,']'
    1594           0 :                valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@state')
    1595           0 :                nprncTemp = 0
    1596           0 :                kappaTemp = 0
    1597           0 :                DO j = 1, 29
    1598           0 :                   IF (TRIM(ADJUSTL(valueString)).EQ.coreStateList_const(j)) THEN
    1599           0 :                      nprncTemp = coreStateNprncList_const(j)
    1600           0 :                      kappaTemp = coreStateKappaList_const(j)
    1601           0 :                      speciesXMLPrintCoreStates(j) = .TRUE.
    1602           0 :                      speciesXMLCoreOccs(1,j) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinUp'))
    1603           0 :                      speciesXMLCoreOccs(2,j) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinDown'))
    1604             :                   END IF
    1605             :                END DO
    1606           0 :                DO j = 1, providedStates
    1607           0 :                   IF ((nprncTemp.EQ.coreStateNprnc(j)).AND.(kappaTemp.EQ.coreStateKappa(j))) THEN
    1608           0 :                      coreStateOccs(j,1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinUp'))
    1609           0 :                      coreStateOccs(j,2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@spinDown'))
    1610             :                   END IF
    1611             :                END DO
    1612             :             END DO
    1613             :          END IF
    1614             : 
    1615             :          ! local orbitals
    1616             : 
    1617          33 :          WRITE(xPathA,*) '/fleurInput/atomSpecies/species[',iSpecies,']/lo'
    1618          33 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA)))
    1619          33 :          iLLO = 1
    1620          42 :          DO iLO = 1, numberNodes
    1621           9 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@l'
    1622           9 :             WRITE(xPathC,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@n'
    1623           9 :             WRITE(xPathD,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@type'
    1624           9 :             WRITE(xPathE,*) TRIM(ADJUSTL(xPathA)),'[',iLO,']/@eDeriv'
    1625           9 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathD)))))
    1626           9 :             lString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))
    1627           9 :             nString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathC)))
    1628           9 :             CALL getIntegerSequenceFromString(TRIM(ADJUSTL(lString)), lNumbers, lNumCount)
    1629           9 :             CALL getIntegerSequenceFromString(TRIM(ADJUSTL(nString)), nNumbers, nNumCount)
    1630           9 :             IF(lNumCount.NE.nNumCount) THEN
    1631           0 :                CALL judft_error('Error in LO input: l quantum number count does not equal n quantum number count')
    1632             :             END IF
    1633           9 :             loEDeriv = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathE))))
    1634          18 :             DO i = 1, lNumCount
    1635           9 :                speciesLLO(iLLO) = lNumbers(i)
    1636           9 :                speciesLOeParams(iLLO) = nNumbers(i)
    1637           9 :                IF(TRIM(ADJUSTL(valueString)).EQ.'HELO') THEN
    1638           0 :                   speciesLOeParams(iLLO) = -speciesLOeParams(iLLO)
    1639             :                END IF
    1640           9 :                speciesLOEDeriv(iLLO) = loEDeriv
    1641          18 :                iLLO = iLLO + 1
    1642             :             END DO
    1643          42 :             DEALLOCATE (lNumbers, nNumbers)
    1644             :          END DO
    1645             : 
    1646             :          ! sort LOs according to l quantum number
    1647             : 
    1648          33 :          ALLOCATE (loOrderList(speciesNLO(iSpecies)),speciesLLOReal(speciesNLO(iSpecies)))
    1649          51 :          DO iLLO = 1, speciesNLO(iSpecies)
    1650          42 :             speciesLLOReal(iLLO) = speciesLLO(iLLO)
    1651             :          END DO
    1652          33 :          CALL sort(loOrderList(:speciesNLO(iSpecies)),speciesLLOReal(:speciesNLO(iSpecies)))
    1653          33 :          DEALLOCATE(speciesLLOReal)
    1654             : 
    1655             :          ! apply species parameters to atom groups
    1656             : 
    1657         114 :          DO iType = 1, atoms%ntype
    1658          81 :             WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']/@species'
    1659          81 :             valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
    1660         114 :             IF(TRIM(ADJUSTL(atoms%speciesName(iSpecies))).EQ.TRIM(ADJUSTL(valueString))) THEN
    1661          55 :                atoms%numStatesProvided(iType) = providedStates
    1662          55 :                IF (coreConfigPresent) THEN
    1663           0 :                   IF (providedCoreStates.NE.atoms%ncst(iType)) THEN
    1664           0 :                      WRITE(6,*) " providedCoreStates:",providedCoreStates
    1665           0 :                      WRITE(6,*) "atoms%ncst(iType)  :",atoms%ncst(iType)
    1666           0 :                      CALL judft_error('Wrong number of core states provided!')
    1667             :                   END IF
    1668           0 :                   DO k = 1, providedStates !atoms%ncst(iType)
    1669           0 :                      atoms%coreStateOccs(k,1,iType) = coreStateOccs(k,1)
    1670           0 :                      atoms%coreStateOccs(k,2,iType) = coreStateOccs(k,2)
    1671           0 :                      atoms%coreStateNprnc(k,iType) = coreStateNprnc(k)
    1672           0 :                      atoms%coreStateKappa(k,iType) = coreStateKappa(k)
    1673           0 :                      xmlElectronStates(k,iType) = speciesXMLElectronStates(k)
    1674           0 :                      xmlPrintCoreStates(k,iType) = speciesXMLPrintCoreStates(k)
    1675           0 :                      xmlCoreOccs (1,k,iType) = speciesXMLCoreOccs(1,k)
    1676           0 :                      xmlCoreOccs (2,k,iType) = speciesXMLCoreOccs(2,k)
    1677             :                   END DO
    1678             :                END IF
    1679          73 :                DO iLLO = 1, speciesNLO(iSpecies)
    1680          18 :                   atoms%llo(iLLO,iType) = speciesLLO(loOrderList(iLLO))
    1681          18 :                   atoms%ulo_der(iLLO,iType) = speciesLOEDeriv(loOrderList(iLLO))
    1682          18 :                   atoms%llod = MAX(ABS(atoms%llo(iLLO,iType)),atoms%llod)
    1683         103 :                   DO jsp = 1, input%jspins
    1684          30 :                      enpara%ello0(iLLO,iType,jsp) = speciesLOeParams(loOrderList(iLLO))
    1685          30 :                      IF (enpara%ello0(iLLO,iType,jsp)==NINT(enpara%ello0(iLLO,iType,jsp))) THEN
    1686          30 :                         enpara%qn_ello(iLLO,iType,jsp)=NINT(enpara%ello0(iLLO,iType,jsp))
    1687          30 :                         enpara%ello0(iLLO,iType,jsp)=0
    1688             :                      ELSE
    1689           0 :                         enpara%qn_ello(iLLO,iType,jsp)=0
    1690             :                      ENDIF
    1691          48 :                      enpara%skiplo(iType,jsp)=enpara%skiplo(iType,jsp)+(2*atoms%llo(iLLO,itype)+1)
    1692             :                   END DO
    1693             :                END DO
    1694             :                ! Energy parameters
    1695         130 :                DO jsp = 1, input%jspins
    1696         675 :                   DO l = 0, 3
    1697         300 :                      enpara%el0(l,iType,jsp) = speciesEParams(l)
    1698         375 :                      IF (enpara%el0(l,iType,jsp)==NINT(enpara%el0(l,iType,jsp))) THEN
    1699         300 :                         enpara%qn_el(l,iType,jsp)=NINT(enpara%el0(l,iType,jsp))
    1700         300 :                         enpara%el0(l,iType,jsp)=0
    1701             :                      ELSE
    1702           0 :                         enpara%qn_el(l,iType,jsp)=0
    1703             :                      ENDIF
    1704             :                   END DO
    1705         545 :                   DO l = 4,atoms%lmax(iType)
    1706         490 :                      enpara%el0(l,iType,jsp) = enpara%el0(3,iType,jsp)
    1707             :                   END DO
    1708             :                END DO
    1709             :                !Hybrid functional stuff
    1710          55 :                hybrid%lcutm1(iType) = 4
    1711          55 :                hybrid%lcutwf(iType) = atoms%lmax(iType) - atoms%lmax(iType) / 10
    1712          55 :                hybrid%select1(:,iType) = (/4, 0, 4, 2 /)
    1713          55 :                IF (hybrid%l_hybrid) THEN
    1714           0 :                   hybrid%lcutm1(iType)=lcutm
    1715           0 :                   hybrid%lcutwf(iType)=lcutwf
    1716           0 :                   hybrid%select1(:,iType)=hybSelect
    1717             :                ENDIF
    1718             :                ! Explicit xc functional
    1719             :                SELECT TYPE(xcpot)
    1720             :                TYPE IS(t_xcpot_inbuild)
    1721          51 :                   xcpot%lda_atom(iType)=ldaSpecies
    1722             :                END SELECT
    1723          55 :                noco%socscale(iType)=socscaleSpecies
    1724          55 :                IF (field%l_b_field) THEN
    1725           0 :                   IF (.NOT.ALLOCATED(field%b_field_mt)) THEN
    1726           0 :                      ALLOCATE(field%b_field_mt(atoms%ntype))
    1727           0 :                      field%b_field_mt=0.0
    1728             :                   ENDIF
    1729           0 :                   field%b_field_mt(itype)=b_field_mtSpecies
    1730             :                ENDIF
    1731             :             END IF
    1732             :          END DO
    1733          33 :          DEALLOCATE(loOrderList)
    1734          57 :          DEALLOCATE(speciesLLO,speciesLOeParams,speciesLOEDeriv)
    1735             :       END DO
    1736             : 
    1737             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1738             : !!! End of species section
    1739             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1740             : 
    1741             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1742             : !!! Start of atomGroup section
    1743             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1744             : 
    1745          24 :       banddos%l_orb = .FALSE.
    1746          24 :       banddos%orbCompAtom = 0
    1747          24 :       atoms%l_geo = .FALSE.
    1748          24 :       atoms%relax = 0
    1749          24 :       na = 0
    1750          24 :       firstAtomOfType = 1
    1751          79 :       DO iType = 1, atoms%ntype
    1752          55 :          WRITE(xPathA,*) '/fleurInput/atomGroups/atomGroup[',iType,']'
    1753             : 
    1754             :          ! Read in force parameters
    1755          55 :          xPathB = TRIM(ADJUSTL(xPathA))//'/force'
    1756          55 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB)))
    1757          55 :          IF (numberNodes.GE.1) THEN
    1758          55 :             atoms%l_geo(iType) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@calculate'))
    1759          55 :             valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@relaxXYZ')
    1760          55 :             READ(valueString,'(3l1)') relaxX, relaxY, relaxZ
    1761          55 :             IF (relaxX) atoms%relax(1,iType) = 1
    1762          55 :             IF (relaxY) atoms%relax(2,iType) = 1
    1763          55 :             IF (relaxZ) atoms%relax(3,iType) = 1
    1764             :          END IF
    1765             : 
    1766             :          ! Obtain number of equivalent atoms
    1767          55 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos')
    1768          55 :          numberNodes = numberNodes + xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/absPos')
    1769          55 :          numberNodes = numberNodes + xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/filmPos')
    1770          55 :          atoms%neq(iType) = numberNodes
    1771             : 
    1772          55 :          IF (iType.GE.2) THEN
    1773             :             firstAtomOfType = firstAtomOfType + atoms%neq(iType-1)
    1774             :          END IF
    1775             : 
    1776             :          ! Read in atom positions
    1777          55 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/relPos')
    1778         103 :          DO i = 1, numberNodes
    1779          48 :             na = na + 1
    1780          48 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'/relPos[',i,']'
    1781          48 :             IF(xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@label').NE.0) THEN
    1782          15 :                atoms%label(na) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@label')
    1783             :             ELSE
    1784          33 :                WRITE(atoms%label(na),'(i0)') na
    1785             :             END IF
    1786          48 :             valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))
    1787          48 :             atoms%taual(1,na) = evaluatefirst(valueString)
    1788          48 :             atoms%taual(2,na) = evaluatefirst(valueString)
    1789          48 :             atoms%taual(3,na) = evaluatefirst(valueString)
    1790          48 :             atoms%pos(:,na) = MATMUL(cell%amat,atoms%taual(:,na))
    1791          48 :             l_orbcomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp'))
    1792          48 :             IF(l_orbcomp) THEN
    1793           0 :                IF(banddos%l_orb) THEN
    1794           0 :                   CALL juDFT_error("Multiple orbcomp flags set.", calledby = "r_inpXML")
    1795             :                END IF
    1796           0 :                banddos%l_orb = .TRUE.
    1797           0 :                banddos%orbCompAtom = na
    1798             :             END IF
    1799         103 :             wannAtomList(na) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@wannier'))
    1800             :          END DO
    1801             : 
    1802          55 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/absPos')
    1803          55 :          DO i = 1, numberNodes
    1804           0 :             na = na + 1
    1805          55 :             CALL judft_error('absPos not yet implemented!')
    1806             :          END DO
    1807             : 
    1808          55 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathA))//'/filmPos')
    1809          81 :          DO i = 1, numberNodes
    1810          26 :             na = na + 1
    1811          26 :             WRITE(xPathB,*) TRIM(ADJUSTL(xPathA)),'/filmPos[',i,']'
    1812          26 :             IF(xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB))//'/@label').NE.0) THEN
    1813          24 :                atoms%label(na) = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@label')
    1814             :             ELSE
    1815           2 :                WRITE(atoms%label(na),'(i0)') na
    1816             :             END IF
    1817          26 :             valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))
    1818          26 :             atoms%taual(1,na) = evaluatefirst(valueString)
    1819          26 :             atoms%taual(2,na) = evaluatefirst(valueString)
    1820          26 :             atoms%taual(3,na) = evaluatefirst(valueString) / cell%amat(3,3)
    1821          26 :             atoms%pos(:,na) = MATMUL(cell%amat,atoms%taual(:,na))
    1822          26 :             l_orbcomp = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@orbcomp'))
    1823          26 :             IF(l_orbcomp) THEN
    1824           0 :                IF(banddos%l_orb) THEN
    1825           0 :                   CALL juDFT_error("Multiple orbcomp flags set.", calledby = "r_inpXML")
    1826             :                END IF
    1827           0 :                banddos%l_orb = .TRUE.
    1828           0 :                banddos%orbCompAtom = na
    1829           0 :                banddos%alpha=0.0;banddos%beta=0.0;banddos%gamma=0.0
    1830           0 :                WRITE(*,*) "Orbcomp-Rotation feature not fully implemented. Please create an issue on gitlab if you need it :-)"
    1831             :             END IF
    1832          81 :             wannAtomList(na) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@wannier'))
    1833             :          END DO
    1834             : 
    1835             :          !Read in atom group specific noco parameters
    1836          55 :          xPathB = TRIM(ADJUSTL(xPathA))//'/nocoParams'
    1837          55 :          numberNodes = xmlGetNumberOfNodes(TRIM(ADJUSTL(xPathB)))
    1838          79 :          IF (numberNodes.GE.1) THEN
    1839          13 :             noco%l_relax(iType) = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@l_relax'))
    1840          13 :             noco%alphInit(iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@alpha'))
    1841          13 :             noco%alph(iType) = noco%alphInit(iType)
    1842          13 :             noco%beta(iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@beta'))
    1843          13 :             noco%b_con(1,iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@b_cons_x'))
    1844          13 :             noco%b_con(2,iType) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathB))//'/@b_cons_y'))
    1845             :          END IF
    1846             :       END DO
    1847             : 
    1848             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1849             : !!! End of atomGroup section
    1850             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1851             : 
    1852             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1853             : !!! Start of force-theorem section
    1854             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1855             : 
    1856          24 :       xPathA = '/fleurInput/forceTheorem'
    1857          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
    1858          24 :       IF (numberNodes.EQ.1) THEN
    1859             :          !Magnetic anisotropy...
    1860           0 :          xPathA = '/fleurInput/forceTheorem/MAE'
    1861           0 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1862           0 :          IF (numberNodes.EQ.1) THEN
    1863           0 :             lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta')
    1864           0 :             nString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi')
    1865           0 :             ALLOCATE(t_forcetheo_mae::forcetheo)
    1866             :             SELECT TYPE(forcetheo)
    1867             :             TYPE IS(t_forcetheo_mae) !this is ok, we just allocated the type...
    1868           0 :                CALL forcetheo%init(cell,sym,lString,nString)
    1869             :             END SELECT
    1870             :          ENDIF
    1871             :          !spin-spiral dispersion
    1872           0 :          xPathA = '/fleurInput/forceTheorem/spinSpiralDispersion'
    1873           0 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1874           0 :          IF (numberNodes.EQ.1) THEN
    1875           0 :             ALLOCATE(t_forcetheo_ssdisp::forcetheo)
    1876             :             SELECT TYPE(forcetheo)
    1877             :             TYPE IS(t_forcetheo_ssdisp) !this is ok, we just allocated the type...
    1878           0 :                CALL forcetheo%init(priv_read_q_list(xPathA))
    1879             :             END SELECT
    1880             :          ENDIF
    1881             :          !dmi
    1882           0 :          xPathA = '/fleurInput/forceTheorem/DMI'
    1883           0 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1884           0 :          IF (numberNodes.EQ.1) THEN
    1885           0 :             lString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@theta')
    1886           0 :             nString=xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@phi')
    1887           0 :             ALLOCATE(t_forcetheo_dmi::forcetheo)
    1888             :             SELECT TYPE(forcetheo)
    1889             :             TYPE IS(t_forcetheo_dmi) !this is ok, we just allocated the type...
    1890           0 :                CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),lstring,nstring)
    1891             :             END SELECT
    1892             :          ENDIF
    1893             :          !jij
    1894           0 :          xPathA = '/fleurInput/forceTheorem/Jij'
    1895           0 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1896           0 :          IF (numberNodes.EQ.1) THEN
    1897           0 :             thetaj=evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@thetaj'))
    1898           0 :             ALLOCATE(t_forcetheo_jij::forcetheo)
    1899             :             SELECT TYPE(forcetheo)
    1900             :             TYPE IS(t_forcetheo_jij) !this is ok, we just allocated the type...
    1901           0 :                CALL forcetheo%init(priv_read_q_list(TRIM(ADJUSTL(xPathA))//'/qVectors'),thetaj,atoms)
    1902             :             END SELECT
    1903             :          ENDIF
    1904             : 
    1905             :       ELSE
    1906          24 :          ALLOCATE(t_forcetheo::forcetheo) !default no forcetheorem type
    1907             :       ENDIF
    1908             : 
    1909             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1910             : !!! End of force-theorem section
    1911             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1912             : 
    1913             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1914             : !!! Start of output section
    1915             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    1916             : 
    1917          24 :       banddos%dos = .FALSE.
    1918          24 :       banddos%band = .FALSE.
    1919          24 :       banddos%vacdos = .FALSE.
    1920          24 :       sliceplot%slice = .FALSE.
    1921          24 :       input%l_coreSpec = .FALSE.
    1922          24 :       input%l_wann = .FALSE.
    1923             : 
    1924          24 :       input%vchk = .FALSE.
    1925          24 :       input%cdinf = .FALSE.
    1926             : 
    1927          24 :       sliceplot%iplot = .FALSE.
    1928          24 :       input%score = .FALSE.
    1929          24 :       sliceplot%plpot = .FALSE.
    1930             : 
    1931          24 :       input%eonly = .FALSE.
    1932          24 :       input%l_bmt = .FALSE.
    1933             : 
    1934          24 :       xPathA = '/fleurInput/output'
    1935          24 :       numberNodes = xmlGetNumberOfNodes(xPathA)
    1936             : 
    1937          24 :       IF (numberNodes.EQ.1) THEN
    1938             : 
    1939             :          ! Read in general output switches
    1940          24 :          banddos%dos = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@dos'))
    1941          24 :          banddos%band = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@band'))
    1942          24 :          banddos%vacdos = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vacdos'))
    1943          24 :          sliceplot%slice = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@slice'))
    1944          24 :          input%l_coreSpec = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@coreSpec'))
    1945          24 :          input%l_wann = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@wannier'))
    1946          24 :          banddos%l_mcd = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@mcd'))
    1947             : 
    1948             :          ! Read in optional switches for checks
    1949             : 
    1950          24 :          xPathA = '/fleurInput/output/checks'
    1951          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1952             : 
    1953          24 :          IF (numberNodes.EQ.1) THEN
    1954          24 :             input%vchk = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@vchk'))
    1955          24 :             input%cdinf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@cdinf'))
    1956             : 
    1957             :          END IF
    1958             : 
    1959             :          ! Read in optional plotting parameters
    1960             : 
    1961          24 :          xPathA = '/fleurInput/output/plotting'
    1962          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1963             : 
    1964          24 :          IF (numberNodes.EQ.1) THEN
    1965          24 :             sliceplot%iplot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@iplot'))
    1966          24 :             input%score = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@score'))
    1967          24 :             sliceplot%plpot = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@plplot'))
    1968             :          END IF
    1969             : 
    1970             :          ! Read in optional specialOutput switches
    1971             : 
    1972          24 :          xPathA = '/fleurInput/output/specialOutput'
    1973          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1974             : 
    1975          24 :          IF (numberNodes.EQ.1) THEN
    1976          24 :             input%eonly = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eonly'))
    1977          24 :             input%l_bmt = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@bmt'))
    1978             :          END IF
    1979             : 
    1980             :          ! Read in optional densityOfStates output parameters
    1981             : 
    1982          24 :          xPathA = '/fleurInput/output/densityOfStates'
    1983          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    1984             : 
    1985          24 :          IF ((banddos%dos).AND.(numberNodes.EQ.0)) THEN
    1986           0 :             CALL juDFT_error("dos is true but densityOfStates parameters are not set!", calledby = "r_inpXML")
    1987             :          END IF
    1988             : 
    1989          24 :          IF (numberNodes.EQ.1) THEN
    1990          24 :             banddos%ndir = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ndir'))
    1991          24 :             banddos%e2_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minEnergy'))
    1992          24 :             banddos%e1_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxEnergy'))
    1993          24 :             banddos%sig_dos = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sigma'))
    1994             :          END IF
    1995             : 
    1996             :          ! Read in optional vacuumDOS parameters
    1997             : 
    1998          24 :          xPathA = '/fleurInput/output/vacuumDOS'
    1999          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2000             : 
    2001          24 :          IF ((banddos%vacdos).AND.(numberNodes.EQ.0)) THEN
    2002           0 :             CALL juDFT_error("vacdos is true but vacDOS parameters are not set!", calledby = "r_inpXML")
    2003             :          END IF
    2004             : 
    2005          24 :          vacuum%layers = 1
    2006          24 :          input%integ = .FALSE.
    2007          24 :          vacuum%starcoeff = .FALSE.
    2008          24 :          vacuum%nstars = 0
    2009          72 :          vacuum%locx = 0.0
    2010          72 :          vacuum%locy = 0.0
    2011          24 :          vacuum%nstm = 0
    2012          24 :          vacuum%tworkf = 0.0
    2013          24 :          IF ((banddos%vacdos).AND.(numberNodes.EQ.1)) THEN
    2014           0 :             vacuum%layers = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@layers'))
    2015           0 :             input%integ = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@integ'))
    2016           0 :             vacuum%starcoeff = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@star'))
    2017           0 :             vacuum%nstars = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nstars'))
    2018           0 :             vacuum%locx(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locx1'))
    2019           0 :             vacuum%locx(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locx2'))
    2020           0 :             vacuum%locy(1) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locy1'))
    2021           0 :             vacuum%locy(2) = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@locy2'))
    2022           0 :             vacuum%nstm = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nstm'))
    2023           0 :             vacuum%tworkf = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@tworkf'))
    2024             :          END IF
    2025          24 :          vacuum%layerd = vacuum%layers
    2026          24 :          ALLOCATE(vacuum%izlay(vacuum%layerd,2))
    2027             : 
    2028             :          ! Read in optional chargeDensitySlicing parameters
    2029             : 
    2030          24 :          xPathA = '/fleurInput/output/chargeDensitySlicing'
    2031          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2032             : 
    2033          24 :          IF ((sliceplot%slice).AND.(numberNodes.EQ.0)) THEN
    2034           0 :             CALL juDFT_error("slice is true but chargeDensitySlicing parameters are not set!", calledby = "r_inpXML")
    2035             :          END IF
    2036             : 
    2037          24 :          IF (numberNodes.EQ.1) THEN
    2038          24 :             sliceplot%kk = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@numkpt'))
    2039          24 :             sliceplot%e1s = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minEigenval'))
    2040          24 :             sliceplot%e2s = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxEigenval'))
    2041          24 :             sliceplot%nnne = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nnne'))
    2042          24 :             input%pallst = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@pallst'))
    2043             :          END IF
    2044             : 
    2045          24 :          IF (banddos%band) THEN
    2046           1 :             banddos%dos=.TRUE.
    2047           1 :             banddos%ndir = -4
    2048           1 :             WRITE(*,*) 'band="T" --> Overriding "dos" and "ndir"!'
    2049             :          ENDIF
    2050             : 
    2051             :          ! Read in optional core spectrum (EELS) input parameters
    2052             : 
    2053          24 :          xPathA = '/fleurInput/output/coreSpectrum'
    2054          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2055             : 
    2056          24 :          IF ((input%l_coreSpec).AND.(numberNodes.EQ.0)) THEN
    2057           0 :             CALL juDFT_error("coreSpec is true but coreSpectrum parameters are not set!", calledby = "r_inpXML")
    2058             :          END IF
    2059             : 
    2060          24 :          IF (numberNodes.EQ.1) THEN
    2061           0 :             coreSpecInput%verb = 0
    2062           0 :             tempBool = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@verbose'))
    2063           0 :             IF(tempBool) coreSpecInput%verb = 1
    2064           0 :             coreSpecInput%ek0 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eKin'))
    2065           0 :             coreSpecInput%atomType = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomType'))
    2066           0 :             coreSpecInput%lx = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@lmax'))
    2067           0 :             coreSpecInput%edge = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@edgeType')))
    2068           0 :             coreSpecInput%emn = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eMin'))
    2069           0 :             coreSpecInput%emx = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@eMax'))
    2070           0 :             tempInt = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@numPoints'))
    2071           0 :             coreSpecInput%ein = (coreSpecInput%emx - coreSpecInput%emn) / (tempInt - 1.0)
    2072           0 :             coreSpecInput%nqphi = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nqphi'))
    2073           0 :             coreSpecInput%nqr = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@nqr'))
    2074           0 :             coreSpecInput%alpha_ex = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@alpha_Ex'))
    2075           0 :             coreSpecInput%beta_ex = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@beta_Ex'))
    2076           0 :             coreSpecInput%I0 = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@I_initial'))
    2077           0 :             xPathB = TRIM(ADJUSTL(xPathA))//'/edgeIndices'
    2078           0 :             xPathB = TRIM(ADJUSTL(xPathB))//'/text()'
    2079           0 :             valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathB)))
    2080           0 :             numTokens = countStringTokens(valueString)
    2081           0 :             coreSpecInput%edgeidx(:) = 0
    2082           0 :             IF(numTokens.GT.SIZE(coreSpecInput%edgeidx)) THEN
    2083           0 :                CALL juDFT_error('More EELS edge indices provided than allowed.',calledby='r_inpXML')
    2084             :             END IF
    2085           0 :             DO i = 1, MAX(numTokens,SIZE(coreSpecInput%edgeidx))
    2086           0 :                coreSpecInput%edgeidx(i) = evaluateFirstIntOnly(popFirstStringToken(valueString))
    2087             :             END DO
    2088             :          END IF
    2089             : 
    2090             :          ! Read in optional Wannier functions parameters
    2091             : 
    2092          24 :          xPathA = '/fleurInput/output/wannier'
    2093          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2094             : 
    2095          24 :          IF ((input%l_wann).AND.(numberNodes.EQ.0)) THEN
    2096           0 :             CALL juDFT_error("wannier is true but Wannier parameters are not set!", calledby = "r_inpXML")
    2097             :          END IF
    2098             : 
    2099          24 :          IF (numberNodes.EQ.1) THEN
    2100           0 :             wann%l_ms = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@ms'))
    2101           0 :             wann%l_sgwf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@sgwf'))
    2102           0 :             wann%l_socgwf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@socgwf'))
    2103           0 :             wann%l_bs_comf = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@bsComf'))
    2104           0 :             wann%l_atomlist = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@atomList'))
    2105             :          END IF
    2106             : 
    2107          24 :          xPathA = '/fleurInput/output/wannier/bandSelection'
    2108          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2109             : 
    2110          24 :          IF (numberNodes.EQ.1) THEN
    2111           0 :             wann%l_byindex=.TRUE.
    2112           0 :             wann%band_min(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@minSpinUp'))
    2113           0 :             wann%band_max(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@maxSpinUp'))
    2114           0 :             xPathA = '/fleurInput/output/wannier/bandSelection/@minSpinDown'
    2115           0 :             numberNodes = xmlGetNumberOfNodes(xPathA)
    2116           0 :             IF (numberNodes.EQ.1) THEN
    2117           0 :                wann%band_min(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))
    2118             :             ELSE
    2119           0 :                wann%band_min(2) = wann%band_min(1)
    2120             :             END IF
    2121           0 :             xPathA = '/fleurInput/output/wannier/bandSelection/@maxSpinDown'
    2122           0 :             numberNodes = xmlGetNumberOfNodes(xPathA)
    2123           0 :             IF (numberNodes.EQ.1) THEN
    2124           0 :                wann%band_max(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))))
    2125             :             ELSE
    2126           0 :                wann%band_max(2) = wann%band_max(1)
    2127             :             END IF
    2128           0 :             wann%l_byindex = .TRUE.
    2129           0 :             IF(input%l_wann) THEN
    2130           0 :                IF (dimension%neigd.NE.-1) THEN
    2131           0 :                   IF (dimension%neigd.LT.MAX(wann%band_max(1),wann%band_max(2))) THEN
    2132           0 :                      dimension%neigd = MAX(wann%band_max(1),wann%band_max(2))
    2133             :                   END IF
    2134             :                END IF
    2135             :             END IF
    2136             :          END IF
    2137             : 
    2138          24 :          xPathA = '/fleurInput/output/wannier/jobList'
    2139          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2140             : 
    2141          24 :          IF (numberNodes.EQ.1) THEN
    2142           0 :             xPathA = '/fleurInput/output/wannier/jobList/text()'
    2143             : 
    2144             :             ! Note: At the moment only 255 characters for the text in this node. Maybe this is not enough.
    2145           0 :             valueString = xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))
    2146           0 :             numTokens = countStringTokens(valueString)
    2147           0 :             ALLOCATE(wann%jobList(numTokens))
    2148           0 :             DO i = 1, numTokens
    2149           0 :                wann%jobList(i) = popFirstStringToken(valueString)
    2150             :             END DO
    2151             :          END IF
    2152             : 
    2153             :          ! Read in optional magnetic circular dichroism parameters
    2154          24 :          xPathA = '/fleurInput/output/magneticCircularDichroism'
    2155          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2156             : 
    2157          24 :          IF ((banddos%l_mcd).AND.(numberNodes.EQ.0)) THEN
    2158           0 :             CALL juDFT_error("mcd is true but magneticCircularDichroism parameters are not set!", calledby = "r_inpXML")
    2159             :          END IF
    2160             : 
    2161          24 :          IF (numberNodes.EQ.1) THEN
    2162           3 :             banddos%e_mcd_lo = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@energyLo'))
    2163           3 :             banddos%e_mcd_up = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@energyUp'))
    2164             :          END IF
    2165             : 
    2166             :          ! Read in optional parameter for unfolding bandstructure of supercell
    2167          24 :          xPathA = '/fleurInput/output/unfoldingBand'
    2168          24 :          numberNodes = xmlGetNumberOfNodes(xPathA)
    2169             : 
    2170          24 :          IF (numberNodes.EQ.1) THEN
    2171           0 :             banddos%unfoldband = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@unfoldBand'))
    2172           0 :             banddos%s_cell_x = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellX'))
    2173           0 :             banddos%s_cell_y = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellY'))
    2174           0 :             banddos%s_cell_z = evaluateFirstOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@supercellZ'))
    2175             :          END IF
    2176             : 
    2177             :       END IF
    2178             : 
    2179             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2180             : !!! End of output section
    2181             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2182             : 
    2183             :       ! Generate / fill wann%atomlist(:) array
    2184          24 :       IF (wann%l_atomlist) THEN
    2185           0 :          absSum = 0
    2186           0 :          DO i = 1, atoms%nat
    2187           0 :             IF (wannAtomList(i)) absSum = absSum + 1
    2188             :          END DO
    2189           0 :          wann%atomlist_num = absSum
    2190           0 :          ALLOCATE(wann%atomlist(wann%atomlist_num))
    2191           0 :          j = 1
    2192           0 :          DO i = 1, atoms%nat
    2193           0 :             IF (wannAtomList(i)) THEN
    2194           0 :                wann%atomlist(j) = i
    2195           0 :                j = j + 1
    2196             :             END IF
    2197             :          END DO
    2198             :       ELSE
    2199          24 :          wann%atomlist_num = atoms%nat
    2200          24 :          ALLOCATE(wann%atomlist(wann%atomlist_num))
    2201          98 :          DO i = 1, atoms%nat
    2202          24 :             wann%atomlist(i) = i
    2203             :          END DO
    2204             :       END IF
    2205             : 
    2206          24 :       DEALLOCATE(wannAtomList)
    2207             : 
    2208             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2209             : !!! Start of non-XML input
    2210             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2211             : 
    2212             :       ! Read in enpara file iff available
    2213             : 
    2214          24 :       l_enpara = .FALSE.
    2215          24 :       INQUIRE (file ='enpara',exist= l_enpara)
    2216          24 :       IF (l_enpara) THEN
    2217           8 :          CALL enpara%READ(atoms,input%jspins,input%film,.FALSE.)
    2218             :       END IF
    2219             : 
    2220             :       ! Read in nocoinp file iff available
    2221          24 :       l_nocoinp = .FALSE.
    2222          24 :       INQUIRE (file ='nocoinp',exist= l_nocoinp)
    2223          24 :       IF (l_nocoinp) THEN
    2224           0 :          CALL inpnoco(atoms,input,vacuum,noco)
    2225             :       END IF
    2226             : 
    2227             : 
    2228             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2229             : !!! End of non-XML input
    2230             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    2231             : 
    2232             :       !CALL xmlFreeResources()
    2233             :       
    2234             :       !WRITE(*,*) 'Reading of inp.xml file finished'
    2235             : 
    2236          24 :       DEALLOCATE(speciesNLO)
    2237             : 
    2238          48 :    END SUBROUTINE r_inpXML
    2239             : 
    2240          24 :    SUBROUTINE setXCParameters(atoms,namex,relcor,jspins,vxc_id_x,vxc_id_c,exc_id_x,exc_id_c,xcpot)
    2241             :       USE m_juDFT
    2242             :       USE m_types
    2243             :       USE m_types_xcpot_inbuild
    2244             :       USE m_types_xcpot_libxc
    2245             : 
    2246             :       IMPLICIT NONE
    2247             :       TYPE(t_atoms),INTENT(IN)          :: atoms
    2248             :       CHARACTER(LEN=*),     INTENT(IN)  :: namex
    2249             :       LOGICAL,              INTENT(IN)  :: relcor
    2250             :       INTEGER,              INTENT(IN)  :: jspins,vxc_id_c,vxc_id_x,exc_id_x,exc_id_c
    2251             :       CLASS(t_xcpot),INTENT(OUT),ALLOCATABLE      :: xcpot
    2252             : 
    2253          24 :       IF (namex(1:5)=='LibXC') THEN
    2254           2 :          ALLOCATE(t_xcpot_libxc::xcpot)
    2255             :       ELSE
    2256          22 :          ALLOCATE(t_xcpot_inbuild::xcpot)
    2257             :       ENDIF
    2258             : 
    2259          24 :       xcpot%kinED%set = .FALSE.
    2260             :       SELECT TYPE(xcpot)
    2261             :       TYPE IS(t_xcpot_inbuild)
    2262          22 :          CALL xcpot%init(namex(1:4),relcor,atoms%ntype)
    2263             :       TYPE IS(t_xcpot_libxc)
    2264           2 :          CALL xcpot%init(jspins,vxc_id_x,vxc_id_c,exc_id_x,exc_id_c)
    2265             :       END SELECT
    2266             : 
    2267          24 :       CALL set_xcpot_usage(xcpot)
    2268          24 :    END SUBROUTINE setXCParameters
    2269             : 
    2270          24 :    SUBROUTINE set_xcpot_usage(xcpot)
    2271             :       use m_judft_usage
    2272             :       USE m_types
    2273             :       USE m_types_xcpot_inbuild
    2274             :       USE m_types_xcpot_libxc
    2275             :       implicit none
    2276             :       class(t_xcpot), intent(in)    :: xcpot
    2277             : 
    2278             :       ! give some information about XC functional to usage.json
    2279             :       ! 1 -> LDA
    2280             :       ! 2 -> GGA
    2281             :       ! 3 -> MetaGGA
    2282             :       ! 4 -> Hybrid functional
    2283          24 :       if(xcpot%vxc_is_lda()) then
    2284           4 :          call add_usage_data("XC-treatment", 1)
    2285           4 :          return
    2286             :       endif
    2287             : 
    2288          20 :       if(xcpot%exc_is_MetaGGA()) then
    2289           0 :          call add_usage_data("XC-treatment", 3)
    2290           0 :          return
    2291             :       endif
    2292             : 
    2293          20 :       if(xcpot%vxc_is_GGA()) then
    2294          20 :          call add_usage_data("XC-treatment", 2)
    2295          20 :          return
    2296             :       endif
    2297             : 
    2298           0 :       if(xcpot%is_hybrid()) then
    2299           0 :          call add_usage_data("XC-treatment", 4)
    2300           0 :          return
    2301             :       endif
    2302             : 
    2303             :    END SUBROUTINE set_xcpot_usage
    2304             : 
    2305          36 :    SUBROUTINE getIntegerSequenceFromString(string, sequence, count)
    2306             :       use m_juDFT_stop
    2307             :       IMPLICIT NONE
    2308             : 
    2309             :       CHARACTER(*),         INTENT(IN)  :: string
    2310             :       INTEGER, ALLOCATABLE, INTENT(OUT) :: sequence(:)
    2311             :       INTEGER,              INTENT(OUT) :: count
    2312             : 
    2313             :       INTEGER :: i, length, start, lastNumber, currentNumber, index
    2314             :       LOGICAL singleNumber, comma, dash
    2315             : 
    2316             :       ! 3 cases: 1. a single number
    2317             :       !          2. number - number
    2318             :       !          3. comma separated numbers
    2319             : 
    2320          36 :       length = LEN(string)
    2321          36 :       count = 0
    2322          36 :       start = 1
    2323          36 :       singleNumber = .TRUE.
    2324          36 :       comma = .FALSE.
    2325          36 :       dash = .FALSE.
    2326          36 :       lastNumber = 0
    2327             :       count = 0
    2328             : 
    2329             :       ! 1. Determine number count
    2330             : 
    2331          72 :       DO i = 1, length
    2332          36 :          SELECT CASE (string(i:i))
    2333             :          CASE ('0':'9')
    2334             :          CASE (',')
    2335           0 :             IF ((start.EQ.i).OR.(dash)) THEN
    2336           0 :                CALL judft_error('String has wrong syntax (in getIntegerSequenceFromString)')
    2337             :             END IF
    2338           0 :             singleNumber = .FALSE.
    2339           0 :             comma = .TRUE.
    2340           0 :             READ(string(start:i-1),*) lastNumber
    2341           0 :             count = count + 1
    2342           0 :             start = i+1
    2343             :          CASE ('-')
    2344           0 :             IF ((start.EQ.i).OR.(dash).OR.(comma)) THEN
    2345           0 :                CALL judft_error('String has wrong syntax (in getIntegerSequenceFromString)')
    2346             :             END IF
    2347           0 :             singleNumber = .FALSE.
    2348           0 :             dash = .TRUE.
    2349           0 :             READ(string(start:i-1),*) lastNumber
    2350           0 :             start = i+1
    2351             :          CASE DEFAULT
    2352          36 :             CALL judft_error('String has wrong syntax (in getIntegerSequenceFromString)')
    2353             :          END SELECT
    2354             :       END DO
    2355          36 :       IF(start.GT.length) THEN
    2356           0 :          CALL judft_error('String has wrong syntax (in getIntegerSequenceFromString)')
    2357             :       END IF
    2358          36 :       READ(string(start:length),*) currentNumber
    2359          36 :       IF (dash) THEN
    2360           0 :          count = currentNumber - lastNumber + 1
    2361             :       ELSE
    2362          36 :          count = count + 1
    2363             :       END IF
    2364             : 
    2365          36 :       IF (ALLOCATED(sequence)) THEN
    2366           0 :          DEALLOCATE(sequence)
    2367             :       END IF
    2368          36 :       ALLOCATE(sequence(count))
    2369             : 
    2370             :       ! 2. Read in numbers iff comma separation ...and store numbers in any case
    2371             : 
    2372          36 :       IF (singleNumber) THEN
    2373          36 :          sequence(1) = currentNumber
    2374           0 :       ELSE IF (dash) THEN
    2375           0 :          DO i = 1, count
    2376           0 :             sequence(i) = lastNumber + i - 1
    2377             :          END DO
    2378             :       ELSE
    2379             :          index = 1
    2380             :          start = 1
    2381           0 :          DO i = 1, length
    2382           0 :             SELECT CASE (string(i:i))
    2383             :             CASE (',')
    2384           0 :                comma = .TRUE.
    2385           0 :                READ(string(start:i-1),*) lastNumber
    2386           0 :                start = i+1
    2387           0 :                SEQUENCE(index) = lastNumber
    2388           0 :                index = index + 1
    2389             :             END SELECT
    2390             :          END DO
    2391           0 :          sequence(index) = currentNumber
    2392             :       END IF
    2393             : 
    2394          36 :    END SUBROUTINE getIntegerSequenceFromString
    2395             : 
    2396           0 :    FUNCTION popFirstStringToken(line) RESULT(firstToken)
    2397             : 
    2398             :       IMPLICIT NONE
    2399             : 
    2400             :       CHARACTER(*), INTENT(INOUT) :: line
    2401             :       CHARACTER(LEN = LEN(line)) :: firstToken
    2402             : 
    2403             :       INTEGER separatorIndex
    2404             : 
    2405           0 :       separatorIndex = 0
    2406           0 :       line = TRIM(ADJUSTL(line))
    2407             : 
    2408           0 :       separatorIndex = INDEX(line,' ')
    2409           0 :       IF (separatorIndex.LE.1) THEN
    2410           0 :          firstToken = ' '
    2411             :       ELSE
    2412           0 :          firstToken = line(1:separatorIndex-1)
    2413           0 :          line = line(separatorIndex+1:)
    2414             :       END IF
    2415             : 
    2416           0 :    END FUNCTION popFirstStringToken
    2417             : 
    2418           0 :    FUNCTION countStringTokens(line) RESULT(tokenCount)
    2419             : 
    2420             :       IMPLICIT NONE
    2421             : 
    2422             :       CHARACTER(*), INTENT(IN) :: line
    2423             :       INTEGER                  :: tokenCount
    2424             : 
    2425           0 :       CHARACTER(LEN=LEN(line)) :: tempLine
    2426             :       INTEGER separatorIndex
    2427             : 
    2428           0 :       tokenCount = 0
    2429             : 
    2430           0 :       tempLine = TRIM(ADJUSTL(line))
    2431             : 
    2432           0 :       DO WHILE (tempLine.NE.'')
    2433           0 :          separatorIndex = 0
    2434           0 :          separatorIndex = INDEX(tempLine,' ')
    2435           0 :          IF (separatorIndex.EQ.0) THEN
    2436           0 :             tempLine = ''
    2437             :          ELSE
    2438           0 :             tempLine = TRIM(ADJUSTL(tempLine(separatorIndex+1:)))
    2439             :          END IF
    2440           0 :          tokenCount = tokenCount + 1
    2441             :       END DO
    2442             : 
    2443           0 :    END FUNCTION countStringTokens
    2444             : 
    2445           0 :    FUNCTION priv_read_q_list(path)RESULT(q)
    2446             :       USE m_xmlIntWrapFort
    2447             :       USE m_calculator
    2448             :       IMPLICIT NONE
    2449             :       CHARACTER(len=*),INTENT(in):: path
    2450             :       REAL,ALLOCATABLE::q(:,:)
    2451             : 
    2452             :       INTEGER:: n,i
    2453             :       CHARACTER(len=256):: xpatha,valueString
    2454             : 
    2455           0 :       n=xmlGetNumberOfNodes(TRIM(ADJUSTL(path))//'/q')
    2456           0 :       ALLOCATE(q(3,n))
    2457           0 :       DO i = 1, n
    2458             :          !PRINT *, path,'/q[',i,']'
    2459           0 :          WRITE(xPathA,"(a,a,i0,a)") TRIM(ADJUSTL(path)),'/q[',i,']'
    2460             :          !PRINT *,xpatha
    2461           0 :          valueString = TRIM(ADJUSTL(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA)))))
    2462             :          !PRINT *,"Q:",valueString
    2463           0 :          READ(valueString,*) q(1,i),q(2,i),q(3,i)
    2464             :       END DO
    2465           0 :    END FUNCTION priv_read_q_list
    2466         154 : END MODULE m_rinpXML

Generated by: LCOV version 1.13