Start of new inpgen

libxc-git @ ca6f7114

-Subproject commit 3cb2231abf1d47fbd8b3e21c8478e9f26a73ce5f
+Subproject commit ca6f7114b9fffe0964ca9e8d24e09ef15b300316

inpgen/new/atompar.F90

+MODULE m_atompar
+  IMPLICIT NONE
+  USE m_judft
+  type t_atompar
+     integer :: id = -1
+     integer :: nucnumber = 0
+     real    :: rmt = 0.0
+     real    :: rmt_min=99.0
+     integer :: jri = 0
+     REAL    :: dx = 0.0
+     REAL    :: bmu = -9999.0
+     integer :: lmax = 0
+     integer :: lnonsph = 0
+     character(len=100)::lo=""
+     character(len=100)::econfig=""
+     character(len=100)::desc=""
+   contains
+     procedure :: add_defaults
+  end type t_atompar
+
+  type(t_atompar),allocatable :: atompar_list(:)
+  integer :: no_of_atompars
+  
+contains
+  subroutine add_defaults(ap)
+    class(t_atompar),intent(inout)::ap
+
+    TYPE(t_atompar):: ap_d
+
+
+    
+    if (ap%rmt>0) then
+       ap_d=find_atompar(ap%nucnumber,ap%rmt)
+    else
+       call judft_error("Defaults ...")
+    endif
+
+    if (ap%jri==0) ap%jri=ap_d%jri
+    if (ap%dx==0) ap%jri=ap_d%dx
+    if (ap%lmax==0) ap%jri=ap_d%lmax
+    if (ap%lnonsph==0) ap%jri=ap_d%lnonsph
+    if (ap%lo=="") ap%lo=ap_d%lo
+    IF (ap%econfig=="") ap%econfig=ap_d%econfig
+    
+    !now generate defaults for missing values
+    if(ap%jri==0) ap%jri = NINT(NINT(330*ap%rmt)*0.5)*2+1
+    if(ap%dx==0) ap%dx=LOG(3200*ap%nucnumber*ap%rmt)/(ap%jri-1)
+    if(ap%lmax==0) then
+       if (ap%rmt<1.8) then
+          ap%lmax=6
+       else if(ap%rmt<2.4) then
+          ap%lmax=8
+       else
+          ap%lmax=10
+       endif
+    endif
+    IF(ap%lnonsph==0) ap%lnonsph=MIN( MAX( (ap%lmax(:)-2),3 ), 8 )
+    IF (ap%lnonsph>ap%lmax) THEN
+       WRITE(*,*) "lnonsph had to be reduced for some atom"
+       ap%lnonsph=ap%lmax
+    ENDIF
+    !check of magnetism
+    IF (ap%bmu<-99.0) THEN
+       IF (ap_d%bmu>-99.0) THEN
+          ap%bmu=ap_d%bmu
+       ELSE
+          ap%bmu = 0.0
+          IF (ap%nucnumber.EQ.24) ap%bmu = 1.0  ! Cr - Ni
+          IF (ap%nucnumber.EQ.25) ap%bmu = 3.5
+          IF (ap%nucnumber.EQ.26) ap%bmu = 2.2
+          IF (ap%nucnumber.EQ.27) ap%bmu = 1.6
+          IF (ap%nucnumber.EQ.28) ap%bmu = 1.1
+          IF (ap%nucnumber.EQ.59) ap%bmu = 2.1  ! Pr - Tm
+          IF (ap%nucnumber.EQ.60) ap%bmu = 3.1
+          IF (ap%nucnumber.EQ.61) ap%bmu = 4.1
+          IF (ap%nucnumber.EQ.62) ap%bmu = 5.1
+          IF (ap%nucnumber.EQ.63) ap%bmu = 7.1
+          IF (ap%nucnumber.EQ.64) ap%bmu = 7.1 
+          IF (ap%nucnumber.EQ.65) ap%bmu = 6.1
+          IF (ap%nucnumber.EQ.66) ap%bmu = 5.1
+          IF (ap%nucnumber.EQ.67) ap%bmu = 4.1
+          IF (ap%nucnumber.EQ.68) ap%bmu = 3.1
+          IF (ap%nucnumber.EQ.69) ap%bmu = 2.1
+       END IF
+    END IF
+
+    !make sure there are no blanks in lo
+    DO n=1,len_TRIM(lo)
+       IF (ad%lo(n:n)=' ') THEN
+          ad%lo(n:LEN(ad%lo)-n)=ad%lo(n+1:)
+          ad%lo(LEN(ad%lo):LEN(ad%lo))=' '
+    ENDDO
+    
+    
+    
+  end subroutine add_defaults
+    
+       
+  subroutine add_atompar(ap)
+    TYPE(t_atompar),INTENT(in),OPTIONAL::ap
+    type(t_atompar),allocatable:: tmp_list(:)
+
+    
+    
+    IF (.NOT.ALLOCATED(atompar_list)) THEN
+       no_of_atompars=0
+       ALLOCATE(atompar_list(100))
+       !Try to read default parameter files
+       CALL read_params("default.econfig")
+       CALL read_params("fleur.econfig")
+       call read_params("my.econfig")
+    else
+       if (no_of_atompars==size(atompar_list)) then
+          !extend the list
+          call move_alloc(atompar_list,tmp_list)
+          allocate(atompar_list(no_of_atompars+50))
+          atompar_list(:no_of_atompars)=tmp_list
+          deallocate(tmp_list)
+       endif
+    END IF
+    IF (PRESENT(ap)) THEN
+       no_of_atompars=no_of_atompars+1
+       atompar_list(no_of_atompars)=ap
+    ENDIF
+  end subroutine add_atompar
+
+  function find_atompar(nucnumber,rmt_max,id)result(ap)
+    integer,intent(in)          :: nucnumber
+    real,intent(in)             :: rmt_max
+    integer,intent(in),optional :: id
+    type(t_atompar)    :: ap
+
+    integer :: n
+
+    !check if there is an id given
+    if (present(id)) then
+       DO n=no_of_atompars,1,-1
+          ap=atompar_list(n)
+          IF (ap%nucnumber==nucnumber.AND.ap%id==id) THEN
+             RETURN
+          ENDIF
+       end DO
+       call judft_error("You specified a specific id for an atom but never defined that")
+    end if
+
+    !Else we search if the MT has been given for this atom
+    DO n=no_of_atompars,1,-1
+       ap=atompar_list(n)
+       if (ap%nucnumber==nucnumber) then
+          IF (ap%rmt>0.AND.ap%rmt<rmt_max) THEN
+             RETURN
+          ENDIF 
+       endif
+    enddo
+
+    !Else we check if there is an atom definition available
+    DO n=no_of_atompars,1,-1
+       ap=atompar_list(n)
+       if (ap%nucnumber==nucnumber) then
+          if (ap%rmt_min<=rmt_max) then
+             ap%rmt=rmt_max
+             return
+          endif
+       endif
+    enddo
+    
+    call judft_error("No possible atomic parameter-set found")
+  end function find_atompar
+
+
+  SUBROUTINE read_params(filename)
+    CHARACTER(len=*),INTENT(in)::filename
+
+    CHARACTER(len=5)::str
+    INTEGER :: id,z,jri,lmax,lnonsph,io_stat
+    REAL    :: rmt,rmt_min,dx,bmu
+    CHARACTER(len=100)::desc,lo,econfig
+    LOGICAL :: l_exist
+    TYPE(t_atompar)::ap
+    NAMELIST /atom/desc,id,z,jri,lmax,lnonsph,rmt,rmt_min,dx,lo,econfig,bmu
+
+    INQUIRE(file=filename,exist=l_exist)
+    IF (.NOT.l_exist) RETURN
+    
+    OPEN(99,file=filename)
+    DO
+       READ(99,*,err=100,END=100) str
+       IF (str=="&atom") THEN
+          BACKSPACE(99)
+          id=-1;z=0;rmt=0.0;rmt_min=99.0;jri=0;dx=0.0;lmax=0;bmu=-9999.0;lnonsph=0;lo='';econfig=''
+          READ(99,atom,iostat=io_stat)
+          IF (io_stat==0) THEN
+             ap=t_atompar(id=id,nucnumber=z,rmt=rmt,rmt_min=rmt_min,jri=jri,dx=dx,bmu=bmu,lmax=lmax,lnonsph=lnonsph,lo=lo,econfig=econfig,desc=desc)
+          ELSE
+             !try old namelist
+             CALL read_params_old(ap)
+          END IF
+          CALL add_atompar(ap)
+       ENDIF
+    ENDDO
+100 CLOSE(99)
+  END SUBROUTINE read_params
+
+  SUBROUTINE read_params_old(ap)
+    !Try to read old namelist
+    TYPE(t_atompar),INTENT(out)::ap
+
+    REAL:: id,z,rmt,dx,bmu
+    INTEGER:: jri,lmax,lnonsph,ncst,nc,io_stat,nz
+    CHARACTER(len=100)::econfig,lo,element,name
+    
+    NAMELIST /atom/ id,z,rmt,dx,jri,lmax,lnonsph,ncst,econfig,bmu,lo,element,name
+
+    id=-9999.9;z=-1.0;rmt=0.0;dx=0.0;jri=0;lmax=0;lnonsph=0;ncst=-1;lo='';econfig='';name='';bmu=-9999.0
+    
+    BACKSPACE(99)
+    READ(99,atom,iostat=io_stat)
+    IF(io_stat.NE.0) THEN
+       BACKSPACE(99)
+       READ(99,*) name
+       WRITE(*,*) name
+       CALL judft_error("Found a &atom namelist in input that was incorrect")
+    END IF
+    
+    !determine nz and id
+    nz=-1
+    IF (element.NE."") THEN
+       nz=element_to_z(element)
+       IF (z.NE.0.AND.nz.NE.FLOOR(z)) CALL judft_error("z and z of specified element differ")
+    ELSE
+       nz=FLOOR(z)
+    ENDIF
+
+    IF (id.NE.-9999.9) THEN
+       IF (nz==-1) nz=FLOOR(id)
+       id=(id-nz)*100
+       IF (id>=100.OR.id<0) CALL judft_error("ID and element (or nuclear number do not fit")
+    ELSE
+       id=-1.0
+    ENDIF
+    IF (nz==-1) CALL judft_error("Neither z nor element specified")
+
+    IF (ncst>-1) CALL judft_warn("ncst is no longer supported as input")
+
+    ap=t_atompar(id=INT(id),nucnumber=nz,rmt=rmt,dx=dx,jri=jri,lmax=lmax,lnonsph=lnonsph,lo=lo,bmu=bmu,econfig=econfig,desc=name)
+    
+    
+  CONTAINS
+    INTEGER FUNCTION element_to_z(element)
+      USE m_constants,ONLY: namat_const
+      IMPLICIT NONE
+      CHARACTER(len=*),INTENT(in)::  element
+      
+      CHARACTER(len=2)   :: ele
+      INTEGER,parameter  :: adiff=ICHAR('A')-ICHAR('a')
+      INTEGER            :: n
+      
+      ele=ADJUSTL(element)
+      
+      IF (ele(1:1)>='A'.AND.ele(1:1)<='Z') ele(1:1)=CHAR(ICHAR(ele(1:1))-adiff)
+      IF (ele(2:2)>='A'.AND.ele(2:2)<='Z') ele(2:2)=CHAR(ICHAR(ele(2:2))-adiff)
+      
+      element_to_z = -1
+      DO n = 0, SIZE(namat_const)-1
+         IF (TRIM(ele) == TRIM(ADJUSTL(namat_const(n)))) THEN
+            element_to_z = n
+            EXIT
+         ENDIF
+      ENDDO
+      
+    END FUNCTION element_to_z
+  END SUBROUTINE read_params_old
+
+  SUBROUTINE dump_list()
+    INTEGER::n
+
+    INTEGER :: id,z,jri,lmax,lnonsph
+    REAL    :: rmt,rmt_min,dx
+    CHARACTER(len=100)::desc,lo,econfig
+    type(t_atompar)::ap
+    NAMELIST /atom/desc,id,z,jri,lmax,lnonsph,rmt,rmt_min,dx,lo,econfig
+
+    WRITE(6,*) "List of defined atomic parameters"
+    DO n=1,no_of_atompars
+       ap=atompar_list(n)
+       id=ap%id
+       z=ap%nucnumber
+       jri=ap%jri
+       lmax=ap%lmax
+       lnonsph=ap%lnonsph
+       rmt=ap%rmt
+       rmt_min=ap%rmt_min
+       dx=ap%dx
+       lo=ap%lo
+       econfig=ap%econfig
+       desc=ap%desc
+       WRITE(6,atom)
+    ENDDO
+    
+  END SUBROUTINE dump_list
+END MODULE m_atompar
+

inpgen/new/bravais_symm.f90

+MODULE m_bravaissymm
+  use m_juDFT
+  !********************************************************************
+  !     determines the point group of the bravais lattice given the
+  !     lattice vectors. the idea is to determine all the lattice
+  !     vectors that have the same length as a_{1,2,3}, and then use
+  !     these to determine the possible rotation matrices.
+  !     these rotation matrices are in lattice coordinates.    mw 12-99
+  !********************************************************************
+CONTAINS
+  SUBROUTINE bravais_symm(cell,nops,mrot)
+
+    IMPLICIT NONE
+
+    !==> Arguments
+    TYPE(t_cell),INTENT(in) :: cell
+    INTEGER, INTENT(OUT) :: nops, mrot(:,:,:)    ! point group operations
+
+    !==> Locals
+    REAL    amet(3,3),b1,b2,b3,d1,d2,d3,dmax,dt
+    INTEGER i,k,k1,k2,k3,m1,m2,m3,n1,n2,n3
+    INTEGER irot(3,3)
+    
+    INTEGER,PARAMETER::neig12=12! max. number of lattice vectors with same length
+    ! (max occurs for close-packed fcc: 12)
+    INTEGER lv1(3,neig12),lv2(3,neig12),lv3(3,neig12)
+
+    REAL, PARAMETER :: eps=1.0e-9
+ 
+
+    !---> distances for the lattice vectors
+    d1 = cell%aamat(1,1)
+    d2 = cell%aamat(2,2)
+    d3 = cell%aamat(3,3)
+    b1 = ( cell%bmat(1,1)/scale(1) )**2 + ( cell%bmat(1,2)/scale(2) )**2 + ( cell%bmat(1,3)/scale(3) )**2
+    b2 = ( cell%bmat(2,1)/scale(1) )**2 + ( cell%bmat(2,2)/scale(2) )**2 + ( cell%bmat(2,3)/scale(3) )**2
+    b3 = ( cell%bmat(3,1)/scale(1) )**2 + ( cell%bmat(3,2)/scale(2) )**2 + ( cell%bmat(3,3)/scale(3) )**2
+
+    !---> determine the cutoffs along each direction a_i:
+    dmax = max( d1,d2,d3)
+
+    m1 = nint( dmax * b1 )
+    m2 = nint( dmax * b2 )
+    m3 = nint( dmax * b3 )
+
+    !---->loop over all possible lattice vectors to find those with the
+    !---->length, i.e., ones that could be rotations
+    n1 = 1
+    n2 = 1
+    n3 = 1
+
+    lv1(1:3,1) = (/ 1,0,0 /)
+    lv2(1:3,1) = (/ 0,1,0 /)
+    lv3(1:3,1) = (/ 0,0,1 /)
+
+    DO k3=-m3,m3
+       DO k2=-m2,m2
+          DO k1=-m1,m1
+
+             dt = distance2(k1,k2,k3)
+
+             !---->    check if the same length
+             IF ( abs( dt - d1 ) < eps ) THEN
+                IF (.not.( k1==1 .and. k2==0 .and. k3==0 ) ) THEN
+                   n1 = n1+1
+                   IF(n1>neig12)  CALL juDFT_error("n1>neig12", calledby ="bravais_symm")
+                   lv1(1,n1) = k1
+                   lv1(2,n1) = k2
+                   lv1(3,n1) = k3
+                ENDIF
+             ENDIF
+
+             IF ( abs( dt - d2 ) < eps ) THEN
+                IF (.not.( k1==0 .and. k2==1 .and. k3==0 ) ) THEN
+                   n2 = n2+1
+                   IF(n2>neig12)  CALL juDFT_error("n2>neig12",calledby="bravais_symm")
+                   lv2(1,n2) = k1
+                   lv2(2,n2) = k2
+                   lv2(3,n2) = k3
+                ENDIF
+             ENDIF
+
+             IF ( abs( dt - d3 ) < eps ) THEN
+                IF (.not.( k1==0 .and. k2==0 .and. k3==1 ) ) THEN
+                   n3 = n3+1
+                   IF(n3>neig12)  CALL juDFT_error("n3>neig12",calledby="bravais_symm")
+                   lv3(1,n3) = k1
+                   lv3(2,n3) = k2
+                   lv3(3,n3) = k3
+                ENDIF
+             ENDIF
+
+          ENDDO
+       ENDDO
+    ENDDO
+
+    !---> the possible rotation matrices are given by the matrix of
+    !---> column vectors of lv_{1,2,3}
+    nops = 0
+    DO k3 = 1,n3
+       DO k2 = 1,n2
+          DO k1 = 1,n1
+
+             !--->          check whether determinant is +/-1 (needs to be for rotation)
+             IF ( abs(mdet(k1,k2,k3)) .NE. 1 ) CYCLE
+
+             !--->          check whether this maintains lengths correctly
+             !--->          if M is the metric, then must have R^T M R = M 
+             irot = reshape( (/ lv1(:,k1),lv2(:,k2),lv3(:,k3) /) , (/ 3,3 /) )
+             IF ( any( abs(matmul( transpose(irot), matmul(cell%aamat,irot) ) - cell%aamat) > eps ) ) CYCLE
+
+             nops = nops + 1
+             IF ( nops > SIZE(mrot,3) )  CALL juDFT_error("nop > size(mrot)", calledby="bravais_symm")
+             mrot(:,:,nops) = irot
+
+          ENDDO
+       ENDDO
+    ENDDO
+
+    WRITE (6,'(//," Point group of the Bravais lattice has ",i2," operations")') nops
+
+    RETURN
+
+  CONTAINS   ! INTERNAL routines
+
+    REAL FUNCTION distance2(l1,l2,l3)
+      !*********************************************************************
+      !     calculates the magnitude square for a vector (l1,l2,l3) given in
+      !     lattice units
+      !*********************************************************************
+      IMPLICIT NONE
+
+      INTEGER, INTENT(IN) :: l1,l2,l3
+
+      distance2 = l1*(l1*cell%aamat(1,1) + 2*l2*cell%aamat(2,1)) + l2*(l2*cell%aamat(2,2) + 2*l3*cell%aamat(3,2)) + l3*(l3*cell%aamat(3,3) + 2*l1*cell%aamat(1,3))
+
+      RETURN
+    END FUNCTION distance2
+
+    INTEGER FUNCTION mdet(k1,k2,k3)
+      !*********************************************************************
+      !     determines the determinant for possible rotation matrix
+      !     ( lv1(:,k1) ; lv2(:,k2) ; lv3(:,k3) )
+      !*********************************************************************
+      IMPLICIT NONE
+
+      INTEGER, INTENT(IN) :: k1,k2,k3
+
+      mdet = lv1(1,k1)*( lv2(2,k2)*lv3(3,k3) - lv2(3,k2)*lv3(2,k3) ) + lv1(2,k1)*( lv2(3,k2)*lv3(1,k3) - lv2(1,k2)*lv3(3,k3) ) + lv1(3,k1)*( lv2(1,k2)*lv3(2,k3) - lv2(2,k2)*lv3(1,k3) )
+
+      RETURN
+    END FUNCTION mdet
+
+  END SUBROUTINE bravais_symm
+END MODULE m_bravaissymm

inpgen/new/check_mt_radii.f90

+!--------------------------------------------------------------------------------
+! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
+! This file is part of FLEUR and available as free software under the conditions
+! of the MIT license as expressed in the LICENSE file in more detail.
+!--------------------------------------------------------------------------------
+
+MODULE m_check_mt_radii
+  USE m_juDFT
+  !---------------------------------------------------------------------
+  !  Check muffin tin radii and determine a reasonable choice for MTRs.
+  !---------------------------------------------------------------------
+CONTAINS
+  SUBROUTINE check_mt_radii(atoms,input,vacuum,cell,oneD,l_test,rmt1,overlap)
+
+    USE m_types
+    USE m_sort
+    USE m_inv3
+    USE m_juDFT
+    IMPLICIT NONE
+    !     ..
+    !     .. Scalar Arguments ..
+    TYPE(t_atoms),INTENT(IN) :: atoms
+    TYPE(t_input),INTENT(IN) :: input
+    TYPE(t_vacuum),INTENT(IN):: vacuum
+    TYPE(t_cell),INTENT(IN)  :: cell
+    TYPE(t_oneD),INTENT(IN)  :: oneD
+    LOGICAL, INTENT (IN)     :: l_test
+    !     ..
+    !     .. Array Arguments ..
+    REAL,    INTENT (OUT),OPTIONAL    :: rmt1(atoms%ntype)
+    REAL,OPTIONAL,INTENT(OUT):: overlap(0:atoms%ntype,atoms%ntype)
+    !     ..
+    !     .. Local Scalars ..
+    INTEGER na,n
+    INTEGER i,j,k,l
+    INTEGER maxCubeAtoms, iAtom, numAtoms, iNeighborAtom, identicalAtoms
+    INTEGER typeA, typeB
+    REAL    sum_r,facA,facB
+    REAL    rmtMax, rmtMin, rmtMaxDefault, rmtDelta
+    REAL    rmtFac, cubeLength, amatAuxDet
+    REAL    maxSqrDist, dist, currentDist
+    LOGICAL error, outOfBounds
+    !     ..
+    !     .. Local Arrays ..
+    REAL    t_rmt(0:103), minRmts(0:103)
+    REAL    amatAux(3,3), invAmatAux(3,3)
+    REAL    taualAux(3,atoms%nat), posAux(3,atoms%nat)
+    REAL    minPos(3), maxPos(3), pos(3), point(3), realCellPos(3)
+    REAL    offsetPos(3)
+    REAL    nearestAtomDists(atoms%ntype)
+    INTEGER nearestAtoms(atoms%ntype)
+    INTEGER sortedDistList(atoms%ntype)
+    INTEGER minCubeIndex(3), maxCubeIndex(3), cubeIndex(3)
+    INTEGER minCellIndices(3), maxCellIndices(3)
+
+    INTEGER, ALLOCATABLE :: numAtomsInCubes(:,:,:)
+    INTEGER, ALLOCATABLE :: atomRefsInCubes(:,:,:,:)
+    INTEGER, ALLOCATABLE :: refCubes(:,:)
+    INTEGER, ALLOCATABLE :: nearestNeighbors(:,:)
+    INTEGER, ALLOCATABLE :: numNearestNeighbors(:)
+    INTEGER, ALLOCATABLE :: neighborAtoms(:)
+    INTEGER, ALLOCATABLE :: distIndexList(:)
+    REAL,    ALLOCATABLE :: posInCubes(:,:,:,:,:)
+    REAL,    ALLOCATABLE :: refPos(:,:)
+    REAL,    ALLOCATABLE :: nearestNeighborDists(:,:)
+    REAL,    ALLOCATABLE :: sqrDistances(:)
+
+    !     Plan for this routine:
+    !     0. Do initializations and set some constants
+    !     1. Locally replace unit cell by an auxiliary unit cell with:
+    !        a) all atoms within the unit cell
+    !        b) basis vectors obtained by lattice reduction of the original cell. 
+    !           [not in 1st (this) version of routine. Can be implemented with the LLL algorithm when needed.]
+    !     2. Get minimal and maximal coordinates within auxiliary unit cell
+    !     3. Construct mesh of cubes covering the auxiliary unit cell and a boundary of width 2*rmtMax + rmtDelta
+    !     4. Fill mesh of cubes with atoms
+    !        a) Store atoms in cubes and representative cube for each atom type
+    !     5. For each atom in auxiliary unit cell select cube and collect shortest distances to other atoms in neighborhood
+    !        a) Sort distances and set MT radii for the atoms
+    !     6. Correct bad choices and set missing MT radii, vacuum distances, and other parameters
+    !     7. Test old MT radii
+
+
+    !     0. Do initializations and set some constants
+    if (present(overlap)) overlap=0.0
+
+    rmtMaxDefault = 2.8
+    rmtMax = rmtMaxDefault
+    rmtMin = 1.0
+    IF (l_test) THEN
+       rmtMax = MAX(rmtMax,MAXVAL(atoms%rmt(:)))
+       rmtMin = MIN(rmtMin,MINVAL(atoms%rmt(:)))
+    END IF
+    rmtDelta = 0.3
+    IF (input%film) THEN
+       rmtFac = 0.95
+    ELSE
+       rmtFac = 0.975
+    ENDIF
+    t_rmt(0:103) = 2.3 ! default value
+    t_rmt(1) = 1.0 ; t_rmt(5:9) = 1.3 ; t_rmt(16:17) = 1.8
+    cubeLength = 2*rmtMax+rmtDelta
+    maxCubeAtoms = (FLOOR(cubeLength / (0.7*2.0*rmtMin)) + 1)**3
+    error  = .FALSE.
+
+    !     1. For the 1st version the auxiliary unit cell is just a copy of the original unit cell with
+    !        all atoms within the cell.
+
+    DO i = 1, 3
+       DO j = 1, 3
+          amatAux(i,j) = cell%amat(i,j)