Add k point density specification option to inpgen

The k point density can now be defined by specifying denX, denY, denZ
parameters in the &kpts line.

global/types.F90

@@ -282,6 +282,7 @@ MODULE m_types
      REAL,ALLOCATABLE ::wtkpt(:)
      INTEGER               ::  nkptf !<k-vectors in full BZ
      INTEGER               ::  nkpt3(3)
+     REAL                  ::  kPointDensity(3) ! only used if k point set is defined as density
      REAL   ,ALLOCATABLE   ::  bkf(:,:)
      INTEGER,ALLOCATABLE   ::  bkp(:)
      INTEGER,ALLOCATABLE   ::  bksym(:)

inpgen/lapw_input.f

@@ -8,7 +8,7 @@
       CONTAINS
       SUBROUTINE lapw_input(
      >                      infh,nline,xl_buffer,bfh,buffer,
-     <                      jspins,kcrel,ndvgrd,nkpt,div,
+     <                      jspins,kcrel,ndvgrd,nkpt,div,kPointDensity,
      <                      frcor,ctail,chng,tria,kmax,gmax,gmaxxc,
      <                      dvac,dtild,tkb,namex,relcor)
      
@@ -19,12 +19,14 @@
       INTEGER, INTENT (OUT) :: jspins,kcrel,ndvgrd,nkpt,div(3)
       LOGICAL, INTENT (OUT) :: frcor,ctail,tria
       REAL,    INTENT (OUT) :: kmax,gmax,gmaxxc,tkb,chng
+      REAL,    INTENT (OUT) :: kPointDensity(3)
       REAL,    INTENT (INOUT) :: dvac,dtild
       CHARACTER(len=4), INTENT (OUT) :: namex
       CHARACTER(len=12),INTENT (OUT) :: relcor
       CHARACTER(len=xl_buffer)       :: buffer
       
       INTEGER iflag,div1,div2,div3,nline,nbuffer,ios
+      REAL denX, denY, denZ
       LOGICAL h_film,h_comp,h_exco,h_kpt,fatalerror,relxc
       CHARACTER(len=4) :: xctyp
 
@@ -33,7 +35,7 @@
      &                  gmax, gmaxxc, kmax
       NAMELIST /exco/   xctyp, relxc 
       NAMELIST /film/   dvac, dtild
-      NAMELIST /kpt/    nkpt, div1, div2, div3, tkb, tria
+      NAMELIST /kpt/    nkpt,div1,div2,div3,tkb,tria,denX,denY,denZ
 
 
       h_film=.false. ; h_comp=.false.
@@ -120,11 +122,15 @@
       ELSEIF (buffer(1:4)=='&kpt') THEN
         IF (h_kpt) CALL err(1)
 
-        div1 = 1 ; div2 = 1 ; div3 = 1
+        div1 = 0 ; div2 = 0 ; div3 = 0
+        denX = 0.0 ; denY = 0.0 ; denZ = 0.0
 
         READ (bfh,kpt,err=912, end=912, iostat=ios)
         h_kpt=.true.
-        div(1) = div1 ; div(2) = div2 ; div(3) = div3 
+        div(1) = div1 ; div(2) = div2 ; div(3) = div3
+        kPointDensity(1) = denX
+        kPointDensity(2) = denY
+        kPointDensity(3) = denZ
 
 !===> end
       ELSEIF (buffer(1:4)=='&end') THEN

inpgen/set_inp.f90

@@ -52,6 +52,7 @@
  
       INTEGER nel,i,j, nkptOld
       REAL    kmax,dtild,dvac1,n1,n2,gam,kmax0,dtild0,dvac0,sumWeight
+      REAL    recVecLength
       LOGICAL l_test,l_gga,l_exists, l_explicit
       REAL     dx0(atoms%ntype), rmtTemp(atoms%ntype)
       REAL     a1Temp(3),a2Temp(3),a3Temp(3) 
@@ -336,7 +337,7 @@
 !
       CALL lapw_input(&
      &                infh,nline,xl_buffer,bfh,buffer,&
-     &                input%jspins,input%kcrel,obsolete%ndvgrd,kpts%nkpt,div,&
+     &                input%jspins,input%kcrel,obsolete%ndvgrd,kpts%nkpt,div,kpts%kPointDensity,&
      &                input%frcor,input%ctail,obsolete%chng,input%tria,input%rkmax,stars%gmax,xcpot%gmaxxc,&
      &                vacuum%dvac,dtild,input%tkb,namex,relcor)
 !
@@ -363,10 +364,27 @@
           kpts%nkpt = MAX(nint((216000/cell%omtil)/sym%nop),1)
         ENDIF
       ENDIF
-      IF((div(1).EQ.0).OR.(div(2).EQ.0)) THEN
-         kpts%specificationType = 1
-      ELSE
+
+      IF((ANY(div(:).NE.0)).AND.(ANY(kpts%kPointDensity(:).NE.0.0))) THEN
+         CALL juDFT_error('Double specification of k point set', calledby = 'set_inp')
+      END IF
+      IF (ANY(div(:).NE.0)) THEN
          kpts%specificationType = 2
+      ELSE IF (ANY(kpts%kPointDensity(:).NE.0.0)) THEN
+         kpts%specificationType = 4
+      ELSE
+         kpts%specificationType = 1
+      END IF
+      IF (kpts%specificationType.EQ.4) THEN
+         DO i = 1, 3
+            IF (kpts%kPointDensity(i).LE.0.0) THEN
+               CALL juDFT_error('Error: Nonpositive kpointDensity provided', calledby = 'set_inp')
+            END IF
+            recVecLength = SQRT(cell%bmat(i,1)**2 + cell%bmat(i,2)**2 + cell%bmat(i,3)**2)
+            kpts%nkpt3(i) = CEILING(kpts%kPointDensity(i) * recVecLength)
+         END DO
+         kpts%nkpt = kpts%nkpt3(1) * kpts%nkpt3(2) * kpts%nkpt3(3)
+         div(:) = kpts%nkpt3(:)
       END IF
 
       IF(TRIM(ADJUSTL(sym%namgrp)).EQ.'any') THEN

io/r_inpXML.F90

@@ -116,7 +116,7 @@ SUBROUTINE r_inpXML(&
   INTEGER            :: speciesEParams(0:3)
   INTEGER            :: mrotTemp(3,3,48)
   REAL               :: tauTemp(3,48)
-  REAL               :: bk(3), kPointDensity(3)
+  REAL               :: bk(3)
   LOGICAL            :: flipSpin, l_eV, invSym, l_qfix, relaxX, relaxY, relaxZ
   LOGICAL            :: l_vca, coreConfigPresent, l_enpara, l_orbcomp
   REAL               :: magMom, radius, logIncrement, qsc(3), latticeScale, dr
@@ -406,13 +406,14 @@ SUBROUTINE r_inpXML(&
   END IF
 
   ! Option kPointDensity
+  kpts%kPointDensity(:) = 0.0
   xPathA = '/fleurInput/calculationSetup/bzIntegration/kPointDensity'
   numberNodes = xmlGetNumberOfNodes(xPathA)
   IF (numberNodes.EQ.1) THEN
      l_kpts = .FALSE.
-     kPointDensity(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denX'))
-     kPointDensity(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denY'))
-     kPointDensity(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denZ'))
+     kpts%kPointDensity(1) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denX'))
+     kpts%kPointDensity(2) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denY'))
+     kpts%kPointDensity(3) = evaluateFirstIntOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@denZ'))
      kpts%l_gamma = evaluateFirstBoolOnly(xmlGetAttributeValue(TRIM(ADJUSTL(xPathA))//'/@gamma'))
      kpts%specificationType = 4
   END IF
@@ -1075,11 +1076,11 @@ SUBROUTINE r_inpXML(&
   ! Construction of k point mesh if kPointDensity is provided
   IF (kpts%specificationType.EQ.4) THEN
      DO i = 1, 3
-        IF (kPointDensity(i).LE.0.0) THEN
+        IF (kpts%kPointDensity(i).LE.0.0) THEN
            CALL juDFT_error('Error: Nonpositive kpointDensity provided', calledby = 'r_inpXML')
         END IF
         recVecLength = SQRT(cell%bmat(i,1)**2 + cell%bmat(i,2)**2 + cell%bmat(i,3)**2)
-        kpts%nkpt3(i) = CEILING(kPointDensity(i) * recVecLength)
+        kpts%nkpt3(i) = CEILING(kpts%kPointDensity(i) * recVecLength)
      END DO
      kpts%nkpt = kpts%nkpt3(1) * kpts%nkpt3(2) * kpts%nkpt3(3)
      kpts%specificationType = 2

io/w_inpXML.f90

@@ -265,10 +265,13 @@ SUBROUTINE w_inpXML(&
          WRITE (fileNum,208) kpts%nkpt,kptGamma
       END IF
 
-   ELSE !(kpts%specificationType.EQ.2)
+   ELSE IF (kpts%specificationType.EQ.2) THEN
 !            <kPointMesh nx="10" ny="10" nz="10" gamma="F"/>
       210 FORMAT('         <kPointMesh nx="',i0,'" ny="',i0,'" nz="',i0,'" gamma="',l1,'"/>')
       WRITE (fileNum,210) div(1),div(2),div(3),kptGamma
+   ELSE !(kpts%specificationType.EQ.4)
+      212 FORMAT('         <kPointDensity denX="',f0.6,'" denY="',f0.6,'" denZ="',f0.6,'" gamma="',l1,'"/>')
+      WRITE (fileNum,212) kpts%kPointDensity(1),kpts%kPointDensity(2),kpts%kPointDensity(3),kptGamma
    END IF
    WRITE (fileNum,'(a)') '      </bzIntegration>'