Merge branch 'develop' of fleur-git:fleur into develop

core/cored.F90

@@ -95,7 +95,7 @@ CONTAINS
        !         rn = rmt(jatom)
        dxx = atoms%dx(jatom)
        bmu = 0.0
-       CALL setcor(jatom,DIMENSION%jspd,atoms,bmu,nst,kappa,nprnc,occ_h)
+       CALL setcor(jatom,DIMENSION%jspd,atoms,input,bmu,nst,kappa,nprnc,occ_h)
        IF ((bmu > 99.)) THEN
           occ(1:nst) = input%jspins *  occ_h(1:nst,jspin)
        ELSE

core/coredr.F90

@@ -73,7 +73,7 @@ CONTAINS
        END DO
     ELSE
        OPEN (58,file='core.dat',form='formatted',status='new')
-       CALL etabinit(atoms,DIMENSION, vr, etab,ntab,ltab,nkmust)
+       CALL etabinit(atoms,DIMENSION,input, vr, etab,ntab,ltab,nkmust)
     END IF
     !
     ncmsh = DIMENSION%msh

core/etabinit.F90

@@ -7,7 +7,7 @@ MODULE m_etabinit
   !     ntab & ltab transport this info to core.F        gb`02
   !------------------------------------------------------------
 CONTAINS
-  SUBROUTINE etabinit(atoms,DIMENSION, vr,&
+  SUBROUTINE etabinit(atoms,DIMENSION,input, vr,&
        etab,ntab,ltab,nkmust)
 
     USE m_constants, ONLY : c_light
@@ -17,6 +17,7 @@ CONTAINS
     IMPLICIT NONE
     TYPE(t_dimension),INTENT(IN)   :: DIMENSION
     TYPE(t_atoms),INTENT(IN)       :: atoms
+    TYPE(t_input),INTENT(IN)       :: input
     !
     !     .. Scalar Arguments ..
     !     ..
@@ -47,7 +48,7 @@ CONTAINS
        rn = atoms%rmt(jatom)
        dxx = atoms%dx(jatom)
        bmu = 0.0
-       CALL setcor(jatom,1,atoms,bmu,nst,kappa,nprnc,occ)
+       CALL setcor(jatom,1,atoms,input,bmu,nst,kappa,nprnc,occ)
        rnot = atoms%rmsh(1,jatom)
        d = EXP(atoms%dx(jatom))
        rn = rnot* (d** (ncmsh-1))

core/setcor.f90

 MODULE m_setcor
   USE m_juDFT
 CONTAINS
-  SUBROUTINE setcor(itype,jspins,atoms,bmu,nst,kappa,nprnc,occ)
+  SUBROUTINE setcor(itype,jspins,atoms,input,bmu,nst,kappa,nprnc,occ)
     !
     !     *****************************************************
     !     sets the values of kappa and occupation numbers of
@@ -9,11 +9,10 @@ CONTAINS
     !         following code by m. weinert  february 1982
     !     *****************************************************
 
-    USE m_types
-    USE m_types
     USE m_types
     IMPLICIT NONE
     TYPE(t_atoms),INTENT(IN)   :: atoms
+    TYPE(t_input),INTENT(IN)   :: input
     !
     !     .. Scalar Arguments ..
     INTEGER,INTENT (IN) :: itype,jspins
@@ -24,9 +23,9 @@ CONTAINS
     REAL,INTENT (OUT)    :: occ(:,:)
     !     ..
     !     .. Local Scalars ..
-    INTEGER iz,jz,jz0,k,n,jspin
+    INTEGER iz,jz,jz0,k,n,m,i,jspin,tempInt
     INTEGER k_h(2),n_h(2)
-    REAL fj,l,bmu_l,o_h(2), fac(2)
+    REAL fj,l,bmu_l,o_h(2), fac(2),tempReal
     LOGICAL l_clf
     CHARACTER(len=13) :: fname
     !     ..
@@ -35,7 +34,7 @@ CONTAINS
     WRITE(fname,"('corelevels.',i2.2)") NINT(atoms%zatom(itype))
     INQUIRE (file=fname, exist=l_clf)
 
-    IF (l_clf) THEN
+    IF (l_clf.AND..NOT.input%l_inpXML) THEN
        OPEN (61,file=fname,form='formatted')
        READ (61,'(i3)') nst
        IF (bmu.LT.0.001) bmu = 999.
@@ -212,5 +211,33 @@ CONTAINS
        ENDIF
     ENDIF
 
+    ! modify default electron configuration according to explicitely provided setting in inp.xml
+    IF(input%l_inpXML) THEN
+       nst = max(nst,atoms%numStatesProvided(itype))
+       DO n = 1, atoms%numStatesProvided(itype)
+          IF((nprnc(n).NE.atoms%coreStateNprnc(n,itype)).OR.(kappa(n).NE.atoms%coreStateKappa(n,itype))) THEN
+             m = 0
+             DO m = n, nst
+                IF((nprnc(m).EQ.atoms%coreStateNprnc(n,itype)).AND.(kappa(m).EQ.atoms%coreStateKappa(n,itype))) THEN
+                   EXIT
+                END IF
+             END DO
+             DO i = m-1, n, -1
+                nprnc(i+1) = nprnc(i)
+                kappa(i+1) = kappa(i)
+                occ(i+1,:) = occ(i,:)
+             END DO
+          END IF
+          nprnc(n) = atoms%coreStateNprnc(n,itype)
+          kappa(n) = atoms%coreStateKappa(n,itype)
+          IF (jspins.EQ.1) THEN
+             occ(n,1) = atoms%coreStateOccs(n,1,itype) + atoms%coreStateOccs(n,2,itype)
+          ELSE
+             occ(n,1) = atoms%coreStateOccs(n,1,itype)
+             occ(n,2) = atoms%coreStateOccs(n,2,itype)
+          END IF
+       END DO
+    END IF
+
   END SUBROUTINE setcor
       END MODULE m_setcor

global/constants.f

       MODULE m_constants
       IMPLICIT NONE
 
+      INTEGER,PARAMETER::noState_const = 0
+      INTEGER,PARAMETER::coreState_const = 1
+      INTEGER,PARAMETER::valenceState_const = 2
       REAL,PARAMETER:: pi_const=3.1415926535897932
       REAL,PARAMETER:: tpi_const=2.*3.1415926535897932
       REAL,PARAMETER:: fpi_const=4.*3.1415926535897932

global/types.F90

@@ -191,8 +191,8 @@
        INTEGER,ALLOCATABLE::jri(:)
        !core states
        INTEGER,ALLOCATABLE::ncst(:)
-       !Are core states explicitely provided?
-       LOGICAL,ALLOCATABLE::coreStatesProvided(:)
+       !How many states are explicitely provided?
+       INTEGER,ALLOCATABLE::numStatesProvided(:)
        !core state occupations
        REAL,ALLOCATABLE::coreStateOccs(:,:,:)
        !core state nprnc

init/efield.f90

@@ -58,7 +58,7 @@
       DO n = 1,atoms%ntype
          IF (atoms%zatom(n).GE.1.0) THEN
             bmu = 0.0
-            CALL setcor(n,1,atoms,bmu, nst,kappa,nprnc,occ)
+            CALL setcor(n,1,atoms,input,bmu, nst,kappa,nprnc,occ)
             DO nc=1,atoms%ncst(n)
               qe=qe+atoms%neq(n)*occ(nc,1)
             ENDDO

inpgen/atom_input.f

@@ -12,7 +12,7 @@
       SUBROUTINE atom_input(
      >                      infh,xl_buffer,buffer,
      >                      jspins,film,idlist,xmlCoreRefOccs,
-     X                      nline,xmlCoreStates,
+     X                      nline,xmlElectronStates,
      X                      xmlPrintCoreStates,xmlCoreOccs,
      <                      nel,atoms,enpara )
 
@@ -20,7 +20,7 @@
       USE m_juDFT_init
       USE m_readrecord
       USE m_setatomcore, ONLY : setatom_bystr, setcore_bystr
-      USE m_constants,   ONLY : namat_const
+      USE m_constants
       USE m_enpara,      ONLY : w_enpara,default_enpara
 
       IMPLICIT NONE
@@ -38,7 +38,7 @@
       REAL   , INTENT (IN)     :: idlist(atoms%ntype)
       REAL   , INTENT (IN)     :: xmlCoreRefOccs(29)
       REAL, INTENT (INOUT)     :: xmlCoreOccs(2,29,atoms%ntype)
-      LOGICAL, INTENT (INOUT)  :: xmlCoreStates(29,atoms%ntype)
+      INTEGER, INTENT (INOUT)  :: xmlElectronStates(29,atoms%ntype)
       LOGICAL, INTENT (INOUT)  :: xmlPrintCoreStates(29,atoms%ntype)
       CHARACTER(len=xl_buffer) :: buffer
 
@@ -415,7 +415,7 @@
                   IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 25  !(7s1/2)
             END SELECT
             IF(xmlCoreStateNumber.EQ.0) STOP 'Invalid core state!'
-            xmlCoreStates(xmlCoreStateNumber,n) = .TRUE.
+            xmlElectronStates(xmlCoreStateNumber,n) = coreState_const
             xmlPrintCoreStates(xmlCoreStateNumber,n) = 
      +         coreocc(i,n).NE.xmlCoreRefOccs(xmlCoreStateNumber)
             SELECT CASE(xmlCoreStateNumber)
@@ -465,6 +465,60 @@ c           in s and p states equal occupation of up and down states
             END IF
             WRITE(27,'(4i3,i4,a1)') coreqn(1,i,n),coreqn(2,i,n),up,dn,
      &                      coreqn(1,i,n),lotype(lval(i,n))
+            xmlCoreStateNumber = 0
+            SELECT CASE(coreqn(1,i,n))
+               CASE (1)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 1   !(1s1/2)
+               CASE (2)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 2   !(2s1/2)
+                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 3    !(2p1/2)
+                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 4   !(2p3/2)
+               CASE (3)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 5   !(3s1/2)
+                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 6    !(3p1/2)
+                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 7   !(3p3/2)
+                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 9    !(3d3/2)
+                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 10  !(3d5/2)
+               CASE (4)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 8   !(4s1/2)
+                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 11   !(4p1/2)
+                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 12  !(4p3/2)
+                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 14   !(4d3/2)
+                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 15  !(4d5/2)
+                  IF(coreqn(2,i,n).EQ.3) xmlCoreStateNumber = 19   !(4f5/2)
+                  IF(coreqn(2,i,n).EQ.-4) xmlCoreStateNumber = 20  !(4f7/2)
+               CASE (5)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 13  !(5s1/2)
+                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 16   !(5p1/2)
+                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 17  !(5p3/2)
+                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 21   !(5d3/2)
+                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 22  !(5d5/2)
+                  IF(coreqn(2,i,n).EQ.3) xmlCoreStateNumber = 26   !(5f5/2)
+                  IF(coreqn(2,i,n).EQ.-4) xmlCoreStateNumber = 27  !(5f7/2)
+               CASE (6)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 18  !(6s1/2)
+                  IF(coreqn(2,i,n).EQ.1) xmlCoreStateNumber = 23   !(6p1/2)
+                  IF(coreqn(2,i,n).EQ.-2) xmlCoreStateNumber = 24  !(6p3/2)
+                  IF(coreqn(2,i,n).EQ.2) xmlCoreStateNumber = 28   !(6d3/2)
+                  IF(coreqn(2,i,n).EQ.-3) xmlCoreStateNumber = 29  !(6d5/2)
+               CASE (7)
+                  IF(coreqn(2,i,n).EQ.-1) xmlCoreStateNumber = 25  !(7s1/2)
+            END SELECT
+            IF(xmlCoreStateNumber.EQ.0) STOP 'Invalid valence state!'
+            xmlElectronStates(xmlCoreStateNumber,n) = valenceState_const
+            xmlPrintCoreStates(xmlCoreStateNumber,n) = 
+     +         coreocc(i,n).NE.xmlCoreRefOccs(xmlCoreStateNumber)
+!            SELECT CASE(xmlCoreStateNumber)
+!               CASE (9:10,14:15,19:22,26:29)
+!                  up = MIN((xmlCoreRefOccs(xmlCoreStateNumber)/2),
+!     +                          coreocc(i,n))
+!                  dn = MAX(0.0,coreocc(i,n)-up)
+!               CASE DEFAULT
+!                  up = CEILING(coreocc(i,n)/2)
+!                  dn = FLOOR(coreocc(i,n)/2)
+!            END SELECT
+            xmlCoreOccs(1,xmlCoreStateNumber,n) = up
+            xmlCoreOccs(2,xmlCoreStateNumber,n) = dn
           ENDDO
           WRITE (6,*) '----------'
 

inpgen/atom_sym.f

@@ -5,7 +5,7 @@
      >                    dispfh,outfh,errfh,dispfn,natmax,
      X                    natin,atomid,atompos,
      X                    ngen,mmrot,ttr,
-     >                    cartesian,symor,as,bs,nop48,
+     >                    cartesian,i_c,symor,as,bs,nop48,
      <                    ntype,nat,nops,mrot,tau,
      <                    neq,ntyrep,zatom,natype,natrep,natmap,pos)
 
@@ -53,6 +53,7 @@
       LOGICAL, INTENT (IN)    :: symor             ! whether to reduce to symmorphic subgroup
       INTEGER, INTENT (INOUT) :: natin             ! formerly 'ntype0'
       INTEGER, INTENT (IN)    :: ngen              ! Number of generators
+      INTEGER, INTENT (IN)    :: i_c               ! centering of lattice
       INTEGER, INTENT (IN)    :: nop48, natmax     ! dimensioning
       INTEGER, INTENT (INOUT) :: mmrot(3,3,nop48)
       REAL,    INTENT (INOUT) :: ttr(3,nop48)
@@ -76,7 +77,7 @@
       REAL    tr(3),tt(3),disp(3,natmax)
 
       INTEGER mp(3,3),mtmp(3,3)
-      REAL    ttau(3),orth(3,3)
+      REAL    ttau(3),orth(3,3),tc(3,3),td(3,3)
 
       INTEGER mmrot2(3,3,ngen)
       REAL    ttr2(3,ngen)
@@ -84,20 +85,49 @@
       LOGICAL l_exist,lclose,l_inipos
       INTEGER n,na,ng,ncyl,nc,no,nop0,nn,nt,i,j,mops
       INTEGER ios,istep0
-      REAL    eps7
       CHARACTER(len=30) :: filen
       REAL,    ALLOCATABLE :: inipos(:,:)
 
-      eps7= 1.0e-7 ; istep0 = 0
+      REAL, PARAMETER :: eps = 1.0e-7, isqrt3 = 1.0/sqrt(3.0), 
+     &                   thrd = 1.0/3.0, mtthrd = -2.0/3.0
+
+      REAL :: lmat(3,3,8)
+      DATA  lmat /  1.0,  0.0,  0.0,      ! 1: primitive     : P
+     &              0.0,  1.0,  0.0,
+     &              0.0,  0.0,  1.0,
+     +             -1.0,  1.0,  1.0,      ! 2: Inverse (F) 
+     &              1.0, -1.0,  1.0,
+     &              1.0,  1.0, -1.0,
+     +              0.0,  1.0,  1.0,      ! 3: Inverse (I)
+     &              1.0,  0.0,  1.0,
+     &              1.0,  1.0,  0.0,
+     +              1.0,  1.0,  0.0,      ! 4: Inverse (hP)
+     &         -isqrt3, isqrt3, 0.0,
+     &              0.0,  0.0,  1.0,
+     +         0.0, isqrt3, -isqrt3,      ! 5: Inverse (hR)
+     &           mtthrd, thrd, thrd,
+     &             thrd, thrd, thrd,
+     +              1.0,  1.0,  0.0,      ! 6: Inverse ( S (C) )
+     &             -1.0,  1.0,  0.0,
+     &              0.0,  0.0,  1.0,
+     +              1.0,  0.0,  1.0,      ! 7: Inverse (B)
+     &              0.0,  1.0,  0.0,
+     &             -1.0,  0.0,  1.0, 
+     +              1.0,  0.0,  0.0,      ! 8: Inverse (A)
+     &              0.0,  1.0, -1.0,
+     &              0.0,  1.0,  1.0/
+
+       istep0 = 0
 !
 !---> take atomic positions and shift to (-1/2,1/2] in lattice coords.
 !
       natin = abs(natin)
       DO n=1,natin
          IF (cartesian) THEN  ! convert to lattice coords. if necessary
-            atompos(:,n) = matmul( bs, atompos(:,n) )
+!            atompos(:,n) = matmul( bs, atompos(:,n) )
+            atompos(:,n) = matmul( lmat(:,:,i_c), atompos(:,n) )
          ENDIF
-         atompos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps7 )
+         atompos(:,n) = atompos(:,n) - anint( atompos(:,n) - eps )
       ENDDO
 
 !--->  store the positions (in lattice coord.s) given in the input file
@@ -124,7 +154,7 @@
             tr = matmul( bs, tr )
           ENDIF
           atompos(:,n) = atompos(:,n) + tr(:)
-          atompos(:,n) = atompos(:,n) - anint( atompos(:,n)- eps7 )
+          atompos(:,n) = atompos(:,n) - anint( atompos(:,n)- eps )
         ENDDO
         CLOSE (dispfh)
         IF ( ios==0 ) THEN
@@ -142,10 +172,16 @@
 !--->   save generators
         IF (cartesian) THEN    ! convert to lattice coords. if necessary
           DO ng = 2, ngen+1
-            mmrot2(:,:,1) = matmul( bs, mmrot(:,:,ng) )
-            mmrot(:,:,ng) = matmul( mmrot2(:,:,1), as )
-            ttr2(:,1) = matmul( bs, ttr(:,ng) )
+!            mmrot2(:,:,1) = matmul( bs, mmrot(:,:,ng) )
+!            mmrot(:,:,ng) = matmul( mmrot2(:,:,1), as )
+            tc = mmrot(:,:,ng)
+            td = matmul( bs, tc ) 
+            tc = matmul( td, as ) 
+            mmrot(:,:,ng) = NINT(tc)
+            write(*,*) i_c, ttr(:,ng)
+            ttr2(:,1) = matmul( lmat(:,:,i_c), ttr(:,ng) )
             ttr(:,ng) = ttr2(:,1)
+            write(*,*) mmrot(:,:,ng),ttr(:,ng)
           ENDDO
         ENDIF
         mmrot2(:,:,1:ngen) = mmrot(:,:,2:ngen+1) 
@@ -214,7 +250,7 @@
       ENDIF
 
 !---> rewrite all the non-primitive translations so in (-1/2,1/2]
-      ttr(:,1:nops) = ttr(:,1:nops) - anint( ttr(:,1:nops)- eps7 )
+      ttr(:,1:nops) = ttr(:,1:nops) - anint( ttr(:,1:nops)- eps )
 
 !--->  allocate arrays for space group information (mod_spgsym)
 !      if( nopd < nops )then
@@ -226,6 +262,7 @@
          mrot(:,:,n) = mmrot(:,:,n)
          tau(:,n) = ttr(:,n)
          index_op(n) = n
+         write(*,*) n,mrot(:,:,n),tau(:,n)
       ENDDO
 
 !--->    check that the group is closed, etc.
@@ -245,7 +282,7 @@
 !--->   reduce symmetry to the largest symmorphic subgroup
         j = 1
         DO i = 1, nops
-          IF ( all ( abs( tau(:,i) ) < eps7 ) ) THEN
+          IF ( all ( abs( tau(:,i) ) < eps ) ) THEN
             IF ( j<i ) then
               mrot(:,:,j) = mrot(:,:,i)
             ENDIF
@@ -276,10 +313,10 @@
          DO n = 1, nat
             tt = ( atompos(:,nt) - tpos(:,n) )
      &           - anint( atompos(:,nt) - tpos(:,n) )
-            IF ( all( abs(tt) < eps7 ) ) THEN
+            IF ( all( abs(tt) < eps ) ) THEN
                icount(n) = icount(n) + 1
                imap(nt) = n
-               IF ( abs( atomid(nt)-atomid(ity(n)) ) < eps7 )
+               IF ( abs( atomid(nt)-atomid(ity(n)) ) < eps )
      &             CYCLE repres_atoms
                CALL juDFT_error("ERROR! mismatch between atoms."
      +              ,calledby ="atom_sym")
@@ -297,14 +334,14 @@
 !--->    loop over operations
          opts: DO no = 2, nops
             tr = matmul( mrot(:,:,no) , atompos(:,nt) ) + tau(:,no)
-            tr = tr - anint( tr - eps7 )
+            tr = tr - anint( tr - eps )
 !--->       check whether this is a new atom
             DO n = 1, nneq(ntype)
                tt = ( tr-tpos(:,nat+n) ) - anint( tr-tpos(:,nat+n) )
-               IF ( all( abs(tt) < eps7 ) ) THEN
+               IF ( all( abs(tt) < eps ) ) THEN
 
                   nn = ity(nat+n)
-                  IF ( abs( atomid(nt)-atomid(nn) ) < eps7 ) CYCLE opts
+                  IF ( abs( atomid(nt)-atomid(nn) ) < eps ) CYCLE opts
                   WRITE (6,'(" Mismatch between atoms and",
      &                       " symmetry input")')
                   CALL juDFT_error("atom_sym: mismatch rotated",calledby

inpgen/crystal.f

@@ -8,7 +8,7 @@
      >                   dbgfh,errfh,outfh,dispfh,dispfn,
      >                   cal_symm,cartesian,symor,oldfleur,
      >                   natin,natmax,nop48,
-     >                   atomid,atompos,a1,a2,a3,aa,scale,noangles,
+     >                   atomid,atompos,a1,a2,a3,aa,scale,noangles,i_c,
      <                   invs,zrfs,invs2,nop,nop2,
      <                   ngen,mmrot,ttr,ntype,nat,nops,
      <                   neq,ntyrep,zatom,natype,natrep,natmap,
@@ -26,7 +26,7 @@
 
 !===> Arguments
       LOGICAL, INTENT(IN)    :: cal_symm,cartesian,oldfleur,noangles
-      INTEGER, INTENT(IN)    :: ngen,natmax,nop48
+      INTEGER, INTENT(IN)    :: ngen,natmax,nop48,i_c
       INTEGER, INTENT(IN)    :: dbgfh,errfh,outfh,dispfh ! file handles, mainly 6
       REAL,    INTENT(IN)    :: aa
       LOGICAL, INTENT(INOUT) :: symor                    ! on input: if true, reduce symmetry if oldfleur