Merge lp:~albertog/siesta/4.0-semic-pol-fix into lp:siesta/4.0

[Alberto Garcia]

Fix some issues with polarization orbitals. Bessel orbs improvements.

A number of issues (mostly cosmetic, or appearing in rare occasions) with
polarization orbitals have been fixed.

* Fix 'n' quantum number for high-l polarization orbitals

  When a polarization orbital had angular momentum l greater than the
  rest of the the basis, its associated principal quantum number n was
  off (it was set to the polarized shell's n).

  This error affected the metadata in the .ion files, the wavefunction
  files, and the PDOS file, so some analysis operations might have
  failed to take into account those orbitals if they were requested
  explicitly (e.g., a COOP calculation with 'O_3d' symbols might show
  the bug; but not one with 'O_d' or 'O').

* Fix for the case of polarization orbitals having the same l as a
  semicore state.

  An example: For Ti, the electronic structure is 3s2 3p6 3d2 4s2 4p0*

  4p0* is a polarization orbital, and 3p6 can be a semicore orbital.
  The outer orbitals for l=0..3 reported by the 'periodic table'
  routines for Ti are:

      4s2 3p6 3d2 4f0

  When using a pseudopotential with 3p6 in the valence (a "semicore"
  state), if 4p did not appear explicitly in a PAO.Basis block (i.e.,
  with the entry for 4s marked as "polarized"), the program did not
  identify 3p as a "semicore state", and thus it did not generate
  *two* KB operators for l=1.

  This has been fixed, together with the wrong 'n' assignment to the
  polarization orbital.

* Improve output of basis specification

  Change misleading 'n=' prefix by 'i=' (it is simply an enumeration)
  and add a proper symbolic nl identifier at the end of the line.

Some cosmetic and usability improvements have also been made to the
handling of Bessel floating orbitals:

* Extend and clarify the handling of Bessel floating orbitals

  For Bessel floating orbitals different 'zetas' correspond to
  progressively excited states (more nodes) for a given l.

  A user might want to reproduce the basis coverage of a given 'nl'
  shell with a suite of bessel functions with increasing number of
  nodes, using the 'zeta' mechanism. Up to now, it was only possible to
  do this for a single 'n' for each l, and the 'n' quantum number
  specified by the user in the (compulsory for Bessel orbitals)
  PAO.Basis block was ignored.

  This patch enables the specification of multiple 'nl' shells for the
  same l, and keeps the user's choice for the 'n' quantum number in
  the output files of the program (.ion, .PDOS, etc).

* Enable short-hand for entering rc information in PAO.Basis block

  If the number of rc's for a given shell is less than the number of
  'zetas', the program will assign the last rc value to the remaining
  zetas, rather than stopping with an error. This is particularly
  useful for Bessel suites of orbitals.

  A new test (Tests/bessel-rich) exemplifies this, as well as the
  new Bessel extensions.

I had a new inspiration and I think I have fixed all the issues in a simpler and less intrusive way.
'n' quantum numbers are now assigned correctly, and the KB issue for Ti is solved.

The underlying problem is documented in an "archaeological note" at the end of the header of basis_specs.f.

The "basis specs" are now more clearly output, with symbolic names for shells (similar to the
output of the psml branch).

Thank you for your patience!

[Nick Papior (nickpapior)]

Looks good to me!

[lp:~albertog/siesta/4.0-semic-pol-fix]

599. By Alberto Garcia on 2019-01-30

Update manual. Extend 'assumed rc' feature to contraction factors

=== modified file 'Docs/siesta.tex'
--- Docs/siesta.tex	2018-07-19 19:31:44 +0000
+++ Docs/siesta.tex	2019-01-30 13:27:17 +0000
@@ -1669,7 +1669,11 @@
 \item[-] {\tt l}: Angular momentum of
 basis orbitals of this shell
 \item[-] {\tt nzls(lsh,is)}: Number of 'zetas' for this shell. For a filteret 
-basis this number is ignored since the number is controlled by the cutoff.
+basis this number is ignored since the number is controlled by the
+cutoff.
+For bessel-floating orbitals, the different 'zetas' map to
+increasingly excited states with the same angular momentum (with
+increasing number of nodes).
 \item[-] {\tt PolOrb(l+1)}: {\it Optional input}. If set equal to {\tt P}, a
 shell of
 polarization functions (with angular momentum $l+1$)  will be constructed
@@ -1718,9 +1722,16 @@
 each 'zeta' for this shell. For the second zeta onwards, if this value
 is negative, the actual rc used will be the given fraction of the
 first zeta's rc.
+If the number of rc's for a given shell is less than the number of
+'zetas', the program will assign the last rc value to the remaining
+zetas, rather than stopping with an error. This is particularly
+useful for Bessel suites of orbitals.
 \item[-] {\tt contrf(izeta,l,is)}: {\it Optional input}.
 Contraction factor\index{scale factor} of
 each 'zeta' for this shell.
+If the number of entries for a given shell is less than the number of
+'zetas', the program will assign the last contraction value to the remaining
+zetas, rather than stopping with an error.
 {\it Default value}: {\tt 1.0}
 \end{itemize}
 
@@ -1751,7 +1762,7 @@
 a constant-pseudopotential atom, {\it i.e.}, the basis functions
 generated will be spherical Bessel functions\index{Bessel functions}
 with the specified $r_c$. In this case, $r_c$ has to be given, as
-{\bf EnergyShift} will not calculate it.\index{basis!Bessel functions}
+{\bf PAO.EnergyShift} will not calculate it.\index{basis!Bessel functions}
 
 Other negative atomic numbers will be interpreted by {\sc Siesta} as
 {\it ghosts}\index{ghost atoms}\index{basis!ghost atoms}

=== modified file 'Src/atom.F'
--- Src/atom.F	2018-04-26 13:04:08 +0000
+++ Src/atom.F	2019-01-30 13:27:17 +0000
@@ -289,8 +289,9 @@
      .      ' ( Floating basis ) ' 
 
         elseif(iz.eq.-100) then
-          write(6,'(a,i4,a)')
-     .      'atom: Called for Z=',iz,'( Floating Bessel functions)'  
+          write(6,'(3a,i4,a,a)')
+     .      'atom: Called for ', atm_label, '  (Z =', iz,')',
+     .      ' ( Floating Bessel functions ) ' 
 
         endif
 ! 
@@ -702,7 +703,7 @@
           slfe(is) = 0.0d0
 ! Calculate Bessel functions
           call Bessel(is,a,b,rofi,drdi,s,
-     .                lmxo,nzeta,rco,lambda,filtercut,no)
+     .                lmxo,nzeta,rco,lambda,filtercut,no,nsemic)
  
           write(6,'(/a,i3)')
      .      'atom: Total number of floating Bessel orbitals:', no
@@ -2578,7 +2579,8 @@
 C
 C       Internal and common variables
 C
-        integer  l, lmax, nzetamax, nkblmx,nsm, nsm_max
+        integer  l, lmax, nzetamax, nkblmx,nsm, nsm_max, n_user
+        integer  npols_lm1
 C
 C Adding a new species to the list
 C          
@@ -2635,27 +2637,42 @@
             
         cnfigtb(:,:,is) = 0
 
-        if (iz.ne.-100) then  
-          do l=0,lmxo
-            do nsm=1, nsemic(l)+1
-              cnfigtb(l,nsm,is)=cnfigmx(l)-(nsemic(l)+1)+nsm
-            enddo
-          enddo 
-          do l=min(lmaxd,lmxo+1),lmaxd
-            cnfigtb(l,1,is)=l+1         ! AG*** Why??
-                                           ! To deal with polarization
-                                           ! orbitals beyond lmxo?
-! BUT, what happens to "in-between" polarization orbitals?
-! 
-!
-          enddo
-        else
-!
-!          Arbitrary "n" for "semicore" bessel states...
-!
+        ! Generate the array with 'n' quantum number information
+        if (iz.ne.-100) then  ! Normal case
+           
+          ! Note that cnfigmx and nsemic are now properly extended
+          ! all the way to lmaxd (in basis_types) to give a proper cnfigtb
+
           do l=0,lmaxd
+             do nsm=1, nsemic(l)+1
+                cnfigtb(l,nsm,is)=cnfigmx(l)-(nsemic(l)+1)+nsm
+             enddo
+
+             if (l==0) cycle
+
+             ! Special case where a polarization orbital has the same l
+             ! as a lower-lying state (e.g. Ti: 3s2 3p6 4s2 3d2 4p0*) 4p0 is a polarization
+             ! orbital, with the same l as 3p.
+             ! (See 'archaeological note' in the header of 'basis_specs.f')
+             npols_lm1 = sum(npolorb(l-1,:))
+             if (npols_lm1 > 0) then
+                nsm = nsemic(l)+2
+                cnfigtb(l,nsm,is)= cnfigmx(l) + 1  ! in the above case, 4p 
+             endif
+            
+          enddo 
+
+        else
+!
+!          For bessel-function states, the 'n' quantum number
+!          does not follow the heuristics for real 'atomic'
+!          states, but since the user has to specify them in the
+!          PAO.Basis block, we use them.           
+
+           do l=0,lmxo
             do nsm=1,nsemic(l)+1
-              cnfigtb(l,nsm,is)=nsm
+              n_user = cnfigmx(l)-(nsemic(l)+1)+nsm
+              cnfigtb(l,nsm,is) = n_user
             enddo
           enddo 
         endif   ! iz.ne.-100
@@ -6727,7 +6744,7 @@
 !
         subroutine BESSEL(is,a,b,rofi,drdi,s,lmxo,
      .                    nzeta,rco,lambda,filtercut,
-     .                    norb) 
+     .                    norb,nsemic) 
 
         use basis_specs, only: restricted_grid
 C
@@ -6744,14 +6761,13 @@
      .    lambda(nzetmx, 0:lmaxd,nsemx),
      .    filtercut(0:lmaxd,nsemx)
 
-
         integer
-     .    lmxo, is, nzeta(0:lmaxd,nsemx), norb
+     .    lmxo, is, nzeta(0:lmaxd,nsemx), norb, nsemic(0:lmaxd)
 C
 C Internal variables
 C
         integer
-     .    l,nprin, nnodes, nodd, nrc, ir, indx, izeta
+     .    l,nprin, nnodes, nodd, nrc, ir, indx, izeta, nsm
 
         real(dp)
      .    rc, dnrm, phi,
@@ -6763,17 +6779,28 @@
         indx=0
         do l=0,lmxo 
 
-          if (nzeta(l,1).gt.0) then
-
-            write(6,'(/2A,I2)')
-     .        'Bessel: floating Bessel functions ',
-     .        'with angular momentum L=',l
-
-            do izeta=1, nzeta(l,1) 
+           do nsm=1,nsemic(l)+1 
+              if (nzeta(l,nsm).gt.0) then
+                 write(6,'(/2A,I2)')
+     .                'Bessel: floating Bessel functions ',
+     .                'with angular momentum L=',l
+                 exit
+              endif
+           enddo 
+
+           do nsm=1,nsemic(l)+1     ! Note extra shell loop
+           
+             if (nzeta(l,nsm).eq.0) CYCLE
+
+             write(6,'(/A,I1,A)')
+     .            'Bessel: Basis orbitals for "state" ',
+     .            cnfigtb(l,nsm,is), sym(l)
+
+             do izeta=1, nzeta(l,nsm) 
 C
 C Cut-off radius for Bessel functions must be an explicit input
 C
-              if (rco(izeta,l,1).lt.1.0d-5) then 
+              if (rco(izeta,l,nsm).lt.1.0d-5) then 
                 write(6,'(a)')
      .     'Bessel: ERROR Zero cut-off radius with Z=-100 option'
                 write(6,'(a)')
@@ -6783,15 +6810,15 @@
                 call die
               endif
 C
-              if (abs(lambda(izeta,l,1)).lt.1.0d-3) 
-     .          lambda(izeta,l,1)=1.0d0
-              if (abs(lambda(izeta,l,1)-1.0d0).gt.1.0d-3) then
+              if (abs(lambda(izeta,l,nsm)).lt.1.0d-3) 
+     .          lambda(izeta,l,nsm)=1.0d0
+              if (abs(lambda(izeta,l,nsm)-1.0d0).gt.1.0d-3) then
                 write(6,'(/,a)')
      . 'Bessel: WARNING Scale factor is not active with Z=-100 option' 
               endif 
-              lambda(izeta,l,1)=1.0d0
+              lambda(izeta,l,nsm)=1.0d0
 
-              rc=rco(izeta,l,1)
+              rc=rco(izeta,l,nsm)
               nrc=nint(log(rc/b+1.0d0)/a)+1
               if (restricted_grid) then
                 nodd=mod(nrc,2)
@@ -6800,10 +6827,17 @@
                 endif
               endif
               rc=b*(exp(a*(nrc-1))-1.0d0)
-              rco(izeta,l,1)=rc
-
-              nnodes=izeta
-              nprin=l+1
+              rco(izeta,l,nsm)=rc
+
+              ! Instead of generating directly j_l(k*r), where
+              ! k is chosen to have a zero at rc, the effective Schrodinger
+              ! equation is integrated with a zero potential.
+
+              ! Note that the the different zetas here mean successively
+              ! excited states for this l
+
+              nnodes=izeta      ! Including the one at r=rc...
+              nprin=l+1       
               call schro_eq(1.0d0,rofi,v,v,s,drdi,
      .                      nrc,l,a,b,nnodes,nprin,
      .                      eorb,g) 
@@ -6826,19 +6860,17 @@
 
               write(6,'(/,(3x,a,i2),2(/,a25,f12.6))')
      .          'izeta =',izeta,
-     .          'rc =',rco(izeta,l,1),
+     .          'rc =',rco(izeta,l,nsm),
      .          'energy =',eorb  
 
               norb=norb+(2*l+1)
               indx=indx+1
               call comBasis(is,a,b,rofi,g,l,
-     .           rco(izeta,l,1),lambda(izeta,l,1),
-     .           filtercut(l,1),izeta,1,nrc,indx)
+     .           rco(izeta,l,nsm),lambda(izeta,l,nsm),
+     .           filtercut(l,nsm),izeta,nsm,nrc,indx)
 
             enddo 
-        
-          endif  
-
+           enddo
         enddo 
 
         end subroutine bessel

=== modified file 'Src/basis_specs.f'
--- Src/basis_specs.f	2017-12-15 10:24:23 +0000
+++ Src/basis_specs.f	2019-01-30 13:27:17 +0000
@@ -47,9 +47,11 @@
 ! 
 !          rc_1  rc_2 .... rc_nzeta  
 ! 
-!   are the cutoff radii in bohrs. This line is mandatory and there must
-!   be at least nzeta values (the extra ones are discarded)
-! 
+!   are the cutoff radii in bohrs. This line is mandatory
+!   If the number of rc's for a given shell is less than the number of
+!   'zetas', the program will assign the last rc value to the remaining
+!   zetas, rather than stopping with an error. This is particularly
+!   useful for Bessel suites of orbitals.
 !   A line containing contraction (scale) factors is optional, but if it
 !   appears, the values *must be* real numbers, and there must be at
 !   least nzeta of them.
@@ -76,7 +78,7 @@
 !   There are no 'per l-shell' polarization orbitals, except if the
 !   third character of 'basis_size' is 'p' (as in 'dzp'), in which
 !   case polarization orbitals are defined so they have the minimum      
-!   angular momentum l such that there are not occupied orbitals 
+!   angular momentum l such that there are no occupied orbitals 
 !   with the same l in the valence shell of the ground-state 
 !   atomic configuration. They polarize the corresponding l-1 shell.
 ! 
@@ -86,7 +88,7 @@
 ! 
 !   rc(1:nzeta) is set to 0.0
 !   lambda(1:nzeta) is set to 1.0  (this is a change from old practice)
-! 
+!
 !  ----------------------------------
 !  
 !   Next come the blocks associated to the KB projectors:
@@ -115,7 +117,28 @@
 !       n-shells in the corresponding PAO shell with the same l. For l
 !       greater than lmxo, it is set to 1. The reference energies are
 !       in all cases set to huge(1.d0) to mark the default.
-!       
+!
+!       Archaeological note: The first implementation of the basis-set generation
+!       module had only non-polarization orbitals. Polarization orbitals were added
+!       later as "second-class" companions. This shows in details like "lmxo" (the
+!       maximum l of the basis set) not taking into account polarization orbitals.
+!       Polarization orbitals were tagged at the end, without maintaining l-shell
+!       ordering.
+!       Even later, support for "semicore" orbitals was added. A new ('nsm') index was
+!       used to distinguish the different orbitals in a l-shell. Polarization orbitals
+!       were not brought into this classification. 
+!       The code is strained when semicore and polarization orbitals coexist for the same l,
+!       as in Ti, whose electronic structure is []3s2 3p6 3d2 4s2 (4p0*)  with 4p as the
+!       polarization orbital.  
+!       Here 'nsemic' (the number of semicore states) is kept at zero for l=1 
+!       (the polarization orbital is not counted), with the side-effect that only one KB
+!       projector is generated for l=1.
+!       Special checks have been implemented to cover these cases.
+! 
+!       Future work should probably remove the separate treatment of polarization orbitals.
+!       (Note that if *all* orbitals are specified in a PAO.Basis block, in effect turning
+!       perturbative polarization orbitals into 'normal' orbitals, this problem is not present.)
+!     
 ! =======================================================================
 !
       use precision
@@ -489,8 +512,21 @@
             k%l = l
             if (l.gt.basp%lmxo) then
               k%nkbl = 1
-            else           ! Set equal to the number of PAO shells
-              k%nkbl = basp%lshell(l)%nn
+            else
+              ! Set equal to the number of PAO shells with this l
+              k%nkbl = basp%lshell(l)%nn     
+              ! Should include polarization orbs (as in Ti case: 3p..4p*)
+              ! (See 'archaeological note' in the header of this file)
+              if (l>0) then
+                 do i = 1, basp%lshell(l-1)%nn
+                    if (basp%lshell(l-1)%shell(i)%polarized) then
+                       k%nkbl = k%nkbl + 1
+                       write(6,"(a,i1,a)") trim(basp%label) //
+     $                  ': nkbl increased for l=',l,
+     $                  ' due to the presence of a polarization orbital'
+                    endif
+                 enddo
+              endif
               if (k%nkbl.eq.0) then
                 write(6,*) 'Warning: Empty PAO shell. l =', l
                 write(6,*) 'Will have a KB projector anyway...'
@@ -592,6 +628,7 @@
       subroutine repaobasis()
 
       integer isp, ish, nn, i, ind, l, indexp, index_splnorm
+      integer nrcs_zetas
 
       type(block_fdf)            :: bfdf
       type(parsed_line), pointer :: pline
@@ -742,11 +779,21 @@
           s%rc(:) = 0.d0
           s%lambda(:) = 1.d0
           if (.not. fdf_bline(bfdf,pline)) call die("No rc's")
-          if (fdf_bnvalues(pline) .ne. s%nzeta)
-     .      call die("Wrong number of rc's")
+          
+          ! Use the last rc entered for the successive zetas
+          ! if there are not enough values (useful for Bessel)
+          nrcs_zetas = fdf_bnvalues(pline)
+          if (nrcs_zetas < 1) then
+           call die("Need at least one rc per shell in PAO.Basis block")
+          endif
           do i= 1, s%nzeta
-            s%rc(i) = fdf_bvalues(pline,i)
+             if (i <= nrcs_zetas) then
+                s%rc(i) = fdf_bvalues(pline,i)
+             else
+                s%rc(i) = s%rc(nrcs_zetas)
+             endif
           enddo
+
           if (s%split_norm_specified) then
             do i = 2,s%nzeta
               if (s%rc(i) /= 0.0_dp) then
@@ -770,11 +817,20 @@
      .          call die('repaobasis: ERROR in PAO.Basis block')
               cycle shells
             else
-              if (fdf_bnreals(pline) .ne. s%nzeta)
-     .          call die("Wrong number of lambda's")
-              do i=1,s%nzeta
-                s%lambda(i) = fdf_breals(pline,i)
-              enddo
+              ! Read scale factors
+              ! Use the last scale factor entered for the successive zetas
+              ! if there are not enough values 
+               nrcs_zetas = fdf_bnreals(pline)
+               if (nrcs_zetas < 1) then
+                 call die("Need at least one scale factor in PAO.Basis")
+               endif
+               do i= 1, s%nzeta
+                  if (i <= nrcs_zetas) then
+                     s%lambda(i) = fdf_breals(pline,i)
+                  else
+                     s%lambda(i) = s%lambda(nrcs_zetas)
+                  endif
+               enddo
             endif
           endif
 

=== modified file 'Src/basis_types.f'
--- Src/basis_types.f	2017-12-15 10:24:23 +0000
+++ Src/basis_types.f	2019-01-30 13:27:17 +0000
@@ -578,6 +578,24 @@
      $           cnfigmx(l,isp) = basp%ground_state%n(l)
 
          enddo
+
+         ! NOTE: cnfigmx and nsemic are only initialized for l up to lmxo
+         !       in the above loop
+         !       Extend them so that we can deal properly with outer polarization states
+
+         do l=basp%lmxo+1, lmaxd
+            !     gs is only setup up to l=3 (f)
+            if (l <= 3) then
+               cnfigmx(l,isp) = basp%ground_state%n(l)
+            else
+               ! g orbitals. Use the "l+1" heuristic
+               ! For example, 5g pol orb associated to a 4f orb.
+               cnfigmx(l,isp) = l + 1
+            endif
+            nsemic(l,isp) = 0
+         enddo
+         
+
          do l=0,basp%lmxkb
             k=>basp%kbshell(l)
             nkbl(l,isp) = k%nkbl
@@ -594,6 +612,10 @@
       integer, intent(in)  :: is
 
       integer l,  n, i
+      integer :: nprin
+      character(len=4) :: orb_id
+      character(len=1), parameter   ::
+     $                           sym(0:4) = (/ 's','p','d','f','g' /)
 
       write(lun,'(/a/79("="))') '<basis_specs>'
          write(lun,'(a20,1x,a2,i4,4x,a5,g12.5,4x,a7,g12.5)')
@@ -609,9 +631,11 @@
      $           'Cnfigmx=', cnfigmx(l,is)
             do n=1,nsemic(l,is)+1
                if (nzeta(l,n,is) == 0) exit
-               write(lun,'(10x,a2,i1,2x,a6,i1,2x,a7,i1)')
-     $                         'n=', n, 'nzeta=',nzeta(l,n,is),
-     $                         'polorb=', polorb(l,n,is)
+               nprin = cnfigmx(l,is) - nsemic(l,is) + n - 1
+               write(orb_id,"(a1,i1,a1,a1)") "(",nprin, sym(l), ")"
+               write(lun,'(10x,a2,i1,2x,a6,i1,2x,a7,i1,2x,a4)')
+     $                         'i=', n, 'nzeta=',nzeta(l,n,is),
+     $                         'polorb=', polorb(l,n,is), orb_id
                if (basistype(is).eq.'filteret') then
                  write(lun,'(10x,a10,2x,g12.5)') 
      $                         'fcutoff:', filtercut(l,n,is)

=== modified file 'Src/old_atmfuncs.f'
--- Src/old_atmfuncs.f	2016-01-25 16:00:16 +0000
+++ Src/old_atmfuncs.f	2019-01-30 13:27:17 +0000
@@ -399,18 +399,17 @@
       integer, intent(in) :: is    ! Species index
       integer, intent(in) :: io    ! Orbital index (within atom)
 
-C Returns the valence-shell configuration in the atomic ground state
-C (i.e. the principal quatum number for orbitals of angular momentum l)
-
-C   INTEGER CNFIGFIO: Principal quantum number of the shell to what 
-C                     the orbital belongs ( for polarization orbitals
-C                     the quantum number corresponds to the shell which
-C                     is polarized by the orbital io) 
-
-      integer l, norb, lorb, izeta, ipol,nsm
-      integer  indx, nsmorb
-
-C
+C   INTEGER CNFIGFIO: Principal quantum number of the shell to which
+C                     the orbital belongs
+
+C               (Formerly, for polarization orbitals
+C               the quantum number corresponded to the shell which
+C               is polarized by the orbital io.
+C               This behavior for polarization orbitals is meaningless,
+C               and has been replaced.)
+      
+      integer l, norb, izeta, ipol, nsm
+
       call check_is('cnfigfio',is)
       if ((io.gt.nomax(is)).or.(io.lt.1)) then
             write(6,*) 'CNFIGFIO: THERE ARE NO DATA FOR IO=',IO
@@ -419,25 +418,32 @@
          call die()
       endif
 
+      ! This routine assumes that polarization orbitals are the last ones
+      ! in the list indexed by 'io'
+      
         norb=0
-        indx=0
         do 10 l=0,lmxosave(is)
          do 8 nsm=1,nsemicsave(l,is)+1
           do 5 izeta=1,nzetasave(l,nsm,is)
             norb=norb+(2*l+1)
-            indx=indx+1
-            if(norb.ge.io) goto 30
+            if(norb.ge.io) then
+               cnfigfio=cnfigtb(l,nsm,is)
+               return
+            endif
  5        continue
  8       continue
 10      continue
 
-        indx=0
         do  20 l=0,lmxosave(is)
           do 18 nsm=1,nsemicsave(l,is)+1
             do 15 ipol=1, npolorbsave(l,nsm,is)
               norb=norb+(2*(l+1)+1)
-              indx=indx+1
-              if(norb.ge.io) goto 40
+              if(norb.ge.io) then
+                 ! Deal properly with polarization orbitals
+                 ! Return the highest n for l+1
+                 cnfigfio=maxval(cnfigtb(l+1,:,is))
+                 return
+              endif
 15          continue
 18        continue
 20      continue
@@ -445,16 +451,6 @@
            write(6,*) 'CNFIGFIO: ERROR: NOT FOUND'
         call die()
 
-30      lorb=l
-        nsmorb=nsm
-        cnfigfio=cnfigtb(lorb,nsmorb,is)
-        return
-
-40      lorb=l
-        nsmorb=nsm
-        cnfigfio=cnfigtb(lorb,nsmorb,is)  
-        return
-
       end function cnfigfio
 !
 !

=== modified file 'Tests/Makefile'
--- Tests/Makefile	2018-04-25 11:14:24 +0000
+++ Tests/Makefile	2019-01-30 13:27:17 +0000
@@ -36,7 +36,7 @@
 #
 label = work
 
-tests = 32_h2o ag anneal-cont ar2_vdw batio3 benzene bessel born \
+tests = 32_h2o ag anneal-cont ar2_vdw batio3 benzene bessel bessel-rich born \
 	born_spin carbon_nanoscroll ch4 chargeconf-h2o constant_volume \
 	dipole_correction fe fe_broyden fe_clust_noncollinear fe_cohp \
 	fen fe_noncol_kp fire_benzene floating force_2 force_constants \

=== added file 'Tests/Reference/bessel-rich.out'
--- Tests/Reference/bessel-rich.out	1970-01-01 00:00:00 +0000
+++ Tests/Reference/bessel-rich.out	2019-01-30 13:27:17 +0000
@@ -0,0 +1,817 @@
+Siesta Version  : siesta-4.0--589--n-fix-4
+Architecture    : gfortran-macosx64-openmpi
+Compiler version: GNU Fortran (Homebrew GCC 7.2.0) 7.2.0
+Compiler flags  : mpif90 -g -O2 
+PP flags        : -DCDF -DMPI -DF2003  
+PARALLEL version
+NetCDF support
+
+* Running in serial mode with MPI
+>> Start of run:  17-JAN-2019  15:47:45
+
+                           ***********************       
+                           *  WELCOME TO SIESTA  *       
+                           ***********************       
+
+reinit: Reading from standard input
+************************** Dump of input data file ****************************
+#
+# Needs new feature: handling of fewer rc's than nzetas in PAO.Basis block
+#
+write-ion-plot-files T
+#
+SystemName          Water molecule with various Bessel Orbitals
+SystemLabel         bessel-rich
+NumberofAtoms   7
+NumberOfSpecies 4
+%block ChemicalSpeciesLabel
+ 1  8  O      # Species index, atomic number, species label
+ 2  1  H
+ 3  -100 Bessel
+ 4  -100 J
+%endblock ChemicalSpeciesLabel
+AtomicCoordinatesFormat  Ang
+%block AtomicCoordinatesAndAtomicSpecies
+ 0.000  0.000  0.000  1
+ 0.757  0.586  0.000  2
+-0.757  0.586  0.000  2
+ 0.3785 0.293  0.000  3
+-0.3785 0.293  0.000  3
+ 0.3785 0.293  0.000  4
+-0.3785 0.293  0.000  4
+%endblock AtomicCoordinatesAndAtomicSpecies
+%block PAO.Basis
+Bessel    3
+ n=1    0    1
+   2.0
+   1.0
+ n=2    0    1
+   2.5
+   1.0
+ n=3    1    1
+   3.5
+   1.0
+J         2  # l-shells
+n=2   0   7    # Note new feature: fewer rc's than zetas
+  4.5
+n=2   1   3
+  4.5 4.5 5.0
+%endblock PAO.Basis
+************************** End of input data file *****************************
+
+reinit: -----------------------------------------------------------------------
+reinit: System Name: Water molecule with various Bessel Orbitals
+reinit: -----------------------------------------------------------------------
+reinit: System Label: bessel-rich                                                 
+reinit: -----------------------------------------------------------------------
+
+initatom: Reading input for the pseudopotentials and atomic orbitals ----------
+ Species number:            1  Label: O Atomic number:           8
+ Species number:            2  Label: H Atomic number:           1
+ Species number:            3  Label: Bessel (floating Bessel functions)
+ Species number:            4  Label: J (floating Bessel functions)
+Ground state valence configuration:   2s02  2p04
+Reading pseudopotential information in formatted form from O.psf
+
+Valence configuration for pseudopotential generation:
+2s( 2.00) rc: 1.14
+2p( 4.00) rc: 1.14
+3d( 0.00) rc: 1.14
+4f( 0.00) rc: 1.14
+Dumping pseudopotential information in formatted form in O.psdump
+Ground state valence configuration:   1s01
+Reading pseudopotential information in formatted form from H.psf
+
+Valence configuration for pseudopotential generation:
+1s( 1.00) rc: 1.25
+2p( 0.00) rc: 1.25
+3d( 0.00) rc: 1.25
+4f( 0.00) rc: 1.25
+Dumping pseudopotential information in formatted form in H.psdump
+For O, standard SIESTA heuristics set lmxkb to 3
+ (one more than the basis l, including polarization orbitals).
+Use PS.lmax or PS.KBprojectors blocks to override.
+For H, standard SIESTA heuristics set lmxkb to 2
+ (one more than the basis l, including polarization orbitals).
+Use PS.lmax or PS.KBprojectors blocks to override.
+
+<basis_specs>
+===============================================================================
+O                    Z=   8    Mass=  16.000        Charge= 0.17977+309
+Lmxo=1 Lmxkb= 3    BasisType=split      Semic=F
+L=0  Nsemic=0  Cnfigmx=2
+          n=1  nzeta=2  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    0.0000      0.0000    
+            lambdas:    1.0000      1.0000    
+L=1  Nsemic=0  Cnfigmx=2
+          n=1  nzeta=2  polorb=1
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    0.0000      0.0000    
+            lambdas:    1.0000      1.0000    
+-------------------------------------------------------------------------------
+L=0  Nkbl=1  erefs: 0.17977+309
+L=1  Nkbl=1  erefs: 0.17977+309
+L=2  Nkbl=1  erefs: 0.17977+309
+L=3  Nkbl=1  erefs: 0.17977+309
+===============================================================================
+</basis_specs>
+
+atom: Called for O                     (Z =   8)
+
+read_vps: Pseudopotential generation method:
+read_vps: ATM3      Troullier-Martins                       
+Total valence charge:    6.00000
+
+xc_check: Exchange-correlation functional:
+xc_check: Ceperley-Alder
+V l=0 = -2*Zval/r beyond r=  1.1278
+V l=1 = -2*Zval/r beyond r=  1.1278
+V l=2 = -2*Zval/r beyond r=  1.1278
+V l=3 = -2*Zval/r beyond r=  1.1138
+All V_l potentials equal beyond r=  1.1278
+This should be close to max(r_c) in ps generation
+All pots = -2*Zval/r beyond r=  1.1278
+
+VLOCAL1: 99.0% of the norm of Vloc inside     34.126 Ry
+VLOCAL1: 99.9% of the norm of Vloc inside     77.774 Ry
+atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge    1.37759
+atom: Maximum radius for r*vlocal+2*Zval:    1.18566
+GHOST: No ghost state for L =  0
+GHOST: No ghost state for L =  1
+GHOST: No ghost state for L =  2
+GHOST: No ghost state for L =  3
+
+KBgen: Kleinman-Bylander projectors: 
+   l= 0   rc=  1.294105   el= -1.742414   Ekb=  9.135903   kbcos=  0.326910
+   l= 1   rc=  1.294105   el= -0.676589   Ekb= -8.124878   kbcos= -0.395047
+   l= 2   rc=  1.448233   el=  0.002386   Ekb= -2.039267   kbcos= -0.003484
+   l= 3   rc=  1.561052   el=  0.003508   Ekb= -0.799141   kbcos= -0.000344
+
+KBgen: Total number of  Kleinman-Bylander projectors:   16
+atom: -------------------------------------------------------------------------
+
+atom: SANKEY-TYPE ORBITALS:
+atom: Selected multiple-zeta basis: split     
+
+SPLIT: Orbitals with angular momentum L= 0
+
+SPLIT: Basis orbitals for state 2s
+
+SPLIT: PAO cut-off radius determined from an
+SPLIT: energy shift=  0.020000 Ry
+
+   izeta = 1
+                 lambda =    1.000000
+                     rc =    3.305093
+                 energy =   -1.723766
+                kinetic =    1.614911
+    potential(screened) =   -3.338677
+       potential(ionic) =  -11.304675
+
+   izeta = 2
+                 rmatch =    2.510382
+              splitnorm =    0.150000
+                 energy =   -1.471299
+                kinetic =    2.446434
+    potential(screened) =   -3.917732
+       potential(ionic) =  -12.476133
+
+SPLIT: Orbitals with angular momentum L= 1
+
+SPLIT: Basis orbitals for state 2p
+
+SPLIT: PAO cut-off radius determined from an
+SPLIT: energy shift=  0.020000 Ry
+
+   izeta = 1
+                 lambda =    1.000000
+                     rc =    3.937239
+                 energy =   -0.658841
+                kinetic =    5.005986
+    potential(screened) =   -5.664827
+       potential(ionic) =  -13.452360
+
+   izeta = 2
+                 rmatch =    2.541963
+              splitnorm =    0.150000
+                 energy =   -0.367441
+                kinetic =    7.530509
+    potential(screened) =   -7.897949
+       potential(ionic) =  -16.611953
+
+POLgen: Perturbative polarization orbital with L=  2
+
+POLgen: Polarization orbital for state 2p
+
+   izeta = 1
+                     rc =    3.937239
+                 energy =    2.398520
+                kinetic =    4.716729
+    potential(screened) =   -2.318209
+       potential(ionic) =   -8.603170
+atom: Total number of Sankey-type orbitals: 13
+
+atm_pop: Valence configuration (for local Pseudopot. screening):
+ 2s( 2.00)                                                            
+ 2p( 4.00)                                                            
+Vna: chval, zval:    6.00000   6.00000
+
+Vna:  Cut-off radius for the neutral-atom potential:   3.937239
+
+atom: _________________________________________________________________________
+
+<basis_specs>
+===============================================================================
+H                    Z=   1    Mass=  1.0100        Charge= 0.17977+309
+Lmxo=0 Lmxkb= 2    BasisType=split      Semic=F
+L=0  Nsemic=0  Cnfigmx=1
+          n=1  nzeta=2  polorb=1
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    0.0000      0.0000    
+            lambdas:    1.0000      1.0000    
+-------------------------------------------------------------------------------
+L=0  Nkbl=1  erefs: 0.17977+309
+L=1  Nkbl=1  erefs: 0.17977+309
+L=2  Nkbl=1  erefs: 0.17977+309
+===============================================================================
+</basis_specs>
+
+atom: Called for H                     (Z =   1)
+
+read_vps: Pseudopotential generation method:
+read_vps: ATM3      Troullier-Martins                       
+Total valence charge:    1.00000
+
+xc_check: Exchange-correlation functional:
+xc_check: Ceperley-Alder
+V l=0 = -2*Zval/r beyond r=  1.2343
+V l=1 = -2*Zval/r beyond r=  1.2189
+V l=2 = -2*Zval/r beyond r=  1.2189
+All V_l potentials equal beyond r=  1.2343
+This should be close to max(r_c) in ps generation
+All pots = -2*Zval/r beyond r=  1.2343
+
+VLOCAL1: 99.0% of the norm of Vloc inside     28.493 Ry
+VLOCAL1: 99.9% of the norm of Vloc inside     64.935 Ry
+atom: Maximum radius for 4*pi*r*r*local-pseudopot. charge    1.45251
+atom: Maximum radius for r*vlocal+2*Zval:    1.21892
+GHOST: No ghost state for L =  0
+GHOST: No ghost state for L =  1
+GHOST: No ghost state for L =  2
+
+KBgen: Kleinman-Bylander projectors: 
+   l= 0   rc=  1.364359   el= -0.467325   Ekb= -2.005361   kbcos= -0.336422
+   l= 1   rc=  1.434438   el=  0.001430   Ekb= -0.501708   kbcos= -0.021697
+   l= 2   rc=  1.470814   el=  0.002365   Ekb= -0.190555   kbcos= -0.002281
+
+KBgen: Total number of  Kleinman-Bylander projectors:    9
+atom: -------------------------------------------------------------------------
+
+atom: SANKEY-TYPE ORBITALS:
+atom: Selected multiple-zeta basis: split     
+
+SPLIT: Orbitals with angular momentum L= 0
+
+SPLIT: Basis orbitals for state 1s
+
+SPLIT: PAO cut-off radius determined from an
+SPLIT: energy shift=  0.020000 Ry
+
+   izeta = 1
+                 lambda =    1.000000
+                     rc =    4.828263
+                 energy =   -0.449375
+                kinetic =    0.929372
+    potential(screened) =   -1.378747
+       potential(ionic) =   -1.915047
+
+   izeta = 2
+                 rmatch =    3.854947
+              splitnorm =    0.150000
+                 energy =   -0.336153
+                kinetic =    1.505294
+    potential(screened) =   -1.841447
+       potential(ionic) =   -2.413582
+
+POLgen: Perturbative polarization orbital with L=  1
+
+POLgen: Polarization orbital for state 1s
+
+   izeta = 1
+                     rc =    4.828263
+                 energy =    0.706972
+                kinetic =    1.396397
+    potential(screened) =   -0.689424
+       potential(ionic) =   -1.169792
+atom: Total number of Sankey-type orbitals:  5
+
+atm_pop: Valence configuration (for local Pseudopot. screening):
+ 1s( 1.00)                                                            
+Vna: chval, zval:    1.00000   1.00000
+
+Vna:  Cut-off radius for the neutral-atom potential:   4.828263
+
+atom: _________________________________________________________________________
+
+<basis_specs>
+===============================================================================
+Bessel               Z=-100    Mass= 0.10000E+41    Charge= 0.17977+309
+Lmxo=1 Lmxkb=-1    BasisType=split      Semic=F
+L=0  Nsemic=1  Cnfigmx=2
+          n=1  nzeta=1  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    2.0000    
+            lambdas:    1.0000    
+          n=2  nzeta=1  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    2.5000    
+            lambdas:    1.0000    
+L=1  Nsemic=0  Cnfigmx=3
+          n=1  nzeta=1  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    3.5000    
+            lambdas:    1.0000    
+-------------------------------------------------------------------------------
+===============================================================================
+</basis_specs>
+
+atom: Called for Bessel                (Z =-100) ( Floating Bessel functions ) 
+
+Bessel: floating Bessel functions with angular momentum L= 0
+
+Bessel: Basis orbitals for "state" 1s
+
+   izeta = 1
+                     rc =    2.011274
+                 energy =    2.439817
+
+Bessel: Basis orbitals for "state" 2s
+
+   izeta = 1
+                     rc =    2.519390
+                 energy =    1.554924
+
+Bessel: floating Bessel functions with angular momentum L= 1
+
+Bessel: Basis orbitals for "state" 3p
+
+   izeta = 1
+                     rc =    3.487864
+                 energy =    1.659713
+
+atom: Total number of floating Bessel orbitals:  5
+
+atom: _________________________________________________________________________
+
+<basis_specs>
+===============================================================================
+J                    Z=-100    Mass= 0.10000E+41    Charge= 0.17977+309
+Lmxo=1 Lmxkb=-1    BasisType=split      Semic=F
+L=0  Nsemic=0  Cnfigmx=2
+          n=1  nzeta=7  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    4.5000      4.5000      4.5000      4.5000    
+            lambdas:    1.0000      1.0000      1.0000      1.0000    
+L=1  Nsemic=0  Cnfigmx=2
+          n=1  nzeta=3  polorb=0
+            splnorm:   0.15000    
+               vcte:    0.0000    
+               rinn:    0.0000    
+               qcoe:    0.0000    
+               qyuk:    0.0000    
+               qwid:   0.10000E-01
+                rcs:    4.5000      4.5000      5.0000    
+            lambdas:    1.0000      1.0000      1.0000    
+-------------------------------------------------------------------------------
+===============================================================================
+</basis_specs>
+
+atom: Called for J                     (Z =-100) ( Floating Bessel functions ) 
+
+Bessel: floating Bessel functions with angular momentum L= 0
+
+Bessel: Basis orbitals for "state" 2s
+
+   izeta = 1
+                     rc =    4.479210
+                 energy =    0.491923
+
+   izeta = 2
+                     rc =    4.479210
+                 energy =    1.967691
+
+   izeta = 3
+                     rc =    4.479210
+                 energy =    4.427304
+
+   izeta = 4
+                     rc =    4.479210
+                 energy =    7.870760
+
+   izeta = 5
+                     rc =    4.479210
+                 energy =   12.298054
+
+   izeta = 6
+                     rc =    4.479210
+                 energy =   17.709175
+
+   izeta = 7
+                     rc =    4.479210
+                 energy =   24.104105
+
+Bessel: floating Bessel functions with angular momentum L= 1
+
+Bessel: Basis orbitals for "state" 2p
+
+   izeta = 1
+                     rc =    4.479210
+                 energy =    1.006350
+
+   izeta = 2
+                     rc =    4.479210
+                 energy =    2.974558
+
+   izeta = 3
+                     rc =    5.075940
+                 energy =    4.614751
+
+atom: Total number of floating Bessel orbitals: 16
+
+atom: _________________________________________________________________________
+
+prinput: Basis input ----------------------------------------------------------
+
+PAO.BasisType split     
+
+%block ChemicalSpeciesLabel
+    1    8 O                       # Species index, atomic number, species label
+    2    1 H                       # Species index, atomic number, species label
+    3 -100 Bessel                  # Species index, atomic number, species label
+    4 -100 J                       # Species index, atomic number, species label
+%endblock ChemicalSpeciesLabel
+
+%block PAO.Basis                 # Define Basis set
+O                     2                    # Species label, number of l-shells
+ n=2   0   2                         # n, l, Nzeta 
+   3.305      2.510   
+   1.000      1.000   
+ n=2   1   2 P   1                   # n, l, Nzeta, Polarization, NzetaPol
+   3.937      2.542   
+   1.000      1.000   
+H                     1                    # Species label, number of l-shells
+ n=1   0   2 P   1                   # n, l, Nzeta, Polarization, NzetaPol
+   4.828      3.855   
+   1.000      1.000   
+Bessel                3                    # Species label, number of l-shells
+ n=1   0   1                         # n, l, Nzeta 
+   2.011   
+   1.000   
+ n=2   0   1                         # n, l, Nzeta 
+   2.519   
+   1.000   
+ n=3   1   1                         # n, l, Nzeta 
+   3.488   
+   1.000   
+J                     2                    # Species label, number of l-shells
+ n=2   0   7                         # n, l, Nzeta 
+   4.479      4.479      4.479      4.479      4.479      4.479      4.479   
+   1.000      1.000      1.000      1.000      1.000      1.000      1.000   
+ n=2   1   3                         # n, l, Nzeta 
+   4.479      4.479      5.076   
+   1.000      1.000      1.000   
+%endblock PAO.Basis
+
+prinput: ----------------------------------------------------------------------
+
+Dumping basis to NetCDF file O.ion.nc
+Dumping basis to NetCDF file H.ion.nc
+Dumping basis to NetCDF file Bessel.ion.nc
+Dumping basis to NetCDF file J.ion.nc
+coor:   Atomic-coordinates input format  =     Cartesian coordinates
+coor:                                          (in Angstroms)
+
+siesta: Atomic coordinates (Bohr) and species
+siesta:      0.00000   0.00000   0.00000  1        1
+siesta:      1.43052   1.10738   0.00000  2        2
+siesta:     -1.43052   1.10738   0.00000  2        3
+siesta:      0.71526   0.55369   0.00000  3        4
+siesta:     -0.71526   0.55369   0.00000  3        5
+siesta:      0.71526   0.55369   0.00000  4        6
+siesta:     -0.71526   0.55369   0.00000  4        7
+
+siesta: Automatic unit cell vectors (Ang):
+siesta:    7.286412    0.000000    0.000000
+siesta:    0.000000    6.087484    0.000000
+siesta:    0.000000    0.000000    5.909356
+
+siesta: System type = molecule  
+
+initatomlists: Number of atoms, orbitals, and projectors:      7    65    34
+
+siesta: ******************** Simulation parameters ****************************
+siesta:
+siesta: The following are some of the parameters of the simulation.
+siesta: A complete list of the parameters used, including default values,
+siesta: can be found in file out.fdf
+siesta:
+redata: Non-Collinear-spin run           =     F
+redata: SpinPolarized (Up/Down) run      =     F
+redata: Number of spin components        =     1
+redata: Long output                      =     F
+redata: Number of Atomic Species         =        4
+redata: Charge density info will appear in .RHO file
+redata: Write Mulliken Pop.              =     NO
+redata: Matel table size (NRTAB)         =     1024
+redata: Mesh Cutoff                      =   100.0000 Ry
+redata: Net charge of the system         =     0.0000 |e|
+redata: Min. number of SCF Iter          =        0
+redata: Max. number of SCF Iter          =       50
+redata: Mix DM or H after convergence    =     F
+redata: Recompute H after scf cycle      =     F
+redata: Mixing is linear
+redata: Mix DM in first SCF step ?       =     F
+redata: Write Pulay info on disk?        =     F
+redata: Discard 1st Pulay DM after  kick =     F
+redata: New DM Mixing Weight             =     0.2500
+redata: New DM Occupancy tolerance       = 0.000000000001
+redata: No kicks to SCF
+redata: DM Mixing Weight for Kicks       =     0.5000
+redata: DM Tolerance for SCF             =     0.000100
+redata: Require (free) Energy convergence in SCF =     F
+redata: DM (free)Energy tolerance for SCF =     0.000010 eV
+redata: Require Harris convergence for SCF =     F
+redata: DM Harris energy tolerance for SCF =     0.000010 eV
+redata: Using Saved Data (generic)   =     F
+redata: Use continuation files for DM    =     F
+redata: Neglect nonoverlap interactions  =     F
+redata: Method of Calculation            =     Diagonalization
+redata: Divide and Conquer               =     T
+redata: Electronic Temperature           =     0.0019  Ry
+redata: Fix the spin of the system       =     F
+redata: Dynamics option                  =     Single-point calculation
+redata: ***********************************************************************
+Total number of electrons:     8.000000
+Total ionic charge:     8.000000
+
+* ProcessorY, Blocksize:    1  24
+
+
+* Orbital distribution balance (max,min):    65    65
+
+ Kpoints in:            1 . Kpoints trimmed:            1
+
+siesta: k-grid: Number of k-points =     1
+siesta: k-grid: Cutoff (effective) =     2.955 Ang
+siesta: k-grid: Supercell and displacements
+siesta: k-grid:    1   0   0      0.000
+siesta: k-grid:    0   1   0      0.000
+siesta: k-grid:    0   0   1      0.000
+
+                     ====================================
+                        Single-point calculation
+                     ====================================
+
+outcell: Unit cell vectors (Ang):
+        7.286412    0.000000    0.000000
+        0.000000    6.087484    0.000000
+        0.000000    0.000000    5.909356
+
+outcell: Cell vector modules (Ang)   :    7.286412    6.087484    5.909356
+outcell: Cell angles (23,13,12) (deg):     90.0000     90.0000     90.0000
+outcell: Cell volume (Ang**3)        :    262.1149
+New_DM. Step:     1
+Initializing Density Matrix...
+New grid distribution:   1
+           1       1:   24    1:   20    1:   18
+
+InitMesh: MESH =    48 x    40 x    36 =       69120
+InitMesh: (bp) =    24 x    20 x    18 =        8640
+InitMesh: Mesh cutoff (required, used) =   100.000   102.571 Ry
+ExtMesh (bp) on 0 =    60 x    56 x    54 =      181440
+PhiOnMesh: Number of (b)points on node 0 =                 8640
+PhiOnMesh: nlist on node 0 =               119152
+
+stepf: Fermi-Dirac step function
+
+siesta: Program's energy decomposition (eV):
+siesta: Ebs     =      -124.892071
+siesta: Eions   =       815.854478
+siesta: Ena     =       175.155695
+siesta: Ekin    =       341.667406
+siesta: Enl     =       -52.736859
+siesta: DEna    =        -0.000003
+siesta: DUscf   =         0.000000
+siesta: DUext   =         0.000000
+siesta: Exc     =      -109.898170
+siesta: eta*DQ  =         0.000000
+siesta: Emadel  =         0.000000
+siesta: Emeta   =         0.000000
+siesta: Emolmec =         0.000000
+siesta: Ekinion =         0.000000
+siesta: Eharris =      -467.297614
+siesta: Etot    =      -461.666410
+siesta: FreeEng =      -461.666410
+
+   scf: iscf   Eharris(eV)      E_KS(eV)   FreeEng(eV)    dDmax  Ef(eV)
+   scf:    1     -467.2976     -461.6664     -461.6664  1.77421 -8.4794
+timer: Routine,Calls,Time,% =    IterSCF     1       0.073   0.70
+   scf:    2     -467.9568     -465.9473     -465.9473  0.17421  0.7795
+   scf:    3     -466.7456     -466.2075     -466.2075  0.03742 -1.0131
+   scf:    4     -466.6648     -466.3262     -466.3262  0.02266 -1.4679
+   scf:    5     -466.6584     -466.4096     -466.4096  0.02027 -1.5893
+   scf:    6     -466.6578     -466.4717     -466.4717  0.01648 -1.6214
+   scf:    7     -466.6578     -466.5182     -466.5182  0.01273 -1.6293
+   scf:    8     -466.6577     -466.5531     -466.5531  0.00965 -1.6311
+   scf:    9     -466.6577     -466.5793     -466.5793  0.00727 -1.6314
+   scf:   10     -466.6577     -466.5989     -466.5989  0.00545 -1.6314
+   scf:   11     -466.6577     -466.6136     -466.6136  0.00409 -1.6315
+   scf:   12     -466.6577     -466.6246     -466.6246  0.00306 -1.6316
+   scf:   13     -466.6577     -466.6329     -466.6329  0.00229 -1.6317
+   scf:   14     -466.6577     -466.6391     -466.6391  0.00171 -1.6318
+   scf:   15     -466.6577     -466.6438     -466.6438  0.00128 -1.6319
+   scf:   16     -466.6577     -466.6472     -466.6472  0.00096 -1.6319
+   scf:   17     -466.6577     -466.6499     -466.6499  0.00072 -1.6320
+   scf:   18     -466.6577     -466.6518     -466.6518  0.00054 -1.6320
+   scf:   19     -466.6577     -466.6533     -466.6533  0.00040 -1.6321
+   scf:   20     -466.6577     -466.6544     -466.6544  0.00030 -1.6321
+   scf:   21     -466.6577     -466.6552     -466.6552  0.00022 -1.6321
+   scf:   22     -466.6577     -466.6559     -466.6559  0.00017 -1.6321
+   scf:   23     -466.6577     -466.6563     -466.6563  0.00012 -1.6322
+   scf:   24     -466.6577     -466.6567     -466.6567  0.00009 -1.6322
+
+SCF Convergence by dMax criterion
+max |DM_out - DM_in|:     0.00009343
+SCF cycle converged after   24 iterations
+
+Using DM_out to compute the final energy and forces
+
+siesta: E_KS(eV) =             -466.6577
+
+siesta: E_KS - E_eggbox =      -466.6577
+
+siesta: Atomic forces (eV/Ang):
+----------------------------------------
+   Tot   -0.000000   -0.032764   -0.000000
+----------------------------------------
+   Max    0.376555
+   Res    0.167768    sqrt( Sum f_i^2 / 3N )
+----------------------------------------
+   Max    0.376555    constrained
+
+Stress-tensor-Voigt (kbar):       -3.42       -0.89       -1.16        0.00        0.00        0.00
+(Free)E + p*V (eV/cell)     -466.3595
+Target enthalpy (eV/cell)     -466.6577
+
+siesta: Program's energy decomposition (eV):
+siesta: Ebs     =      -107.092904
+siesta: Eions   =       815.854478
+siesta: Ena     =       175.155695
+siesta: Ekin    =       351.559238
+siesta: Enl     =       -63.040421
+siesta: DEna    =        -2.430290
+siesta: DUscf   =         0.775212
+siesta: DUext   =         0.000000
+siesta: Exc     =      -112.822678
+siesta: eta*DQ  =         0.000000
+siesta: Emadel  =         0.000000
+siesta: Emeta   =         0.000000
+siesta: Emolmec =         0.000000
+siesta: Ekinion =         0.000000
+siesta: Eharris =      -466.657724
+siesta: Etot    =      -466.657724
+siesta: FreeEng =      -466.657724
+
+siesta: Final energy (eV):
+siesta:  Band Struct. =    -107.092904
+siesta:       Kinetic =     351.559238
+siesta:       Hartree =     388.214304
+siesta:    Ext. field =       0.000000
+siesta:   Exch.-corr. =    -112.822678
+siesta:  Ion-electron =   -1087.255665
+siesta:       Ion-ion =      -6.352922
+siesta:       Ekinion =       0.000000
+siesta:         Total =    -466.657724
+
+siesta: Atomic forces (eV/Ang):
+siesta:      1    0.000000   -0.018600    0.000000
+siesta:      2    0.376555    0.275107   -0.000000
+siesta:      3   -0.376555    0.275107   -0.000000
+siesta:      4   -0.001135   -0.003186    0.000000
+siesta:      5    0.001135   -0.003186    0.000000
+siesta:      6   -0.005249   -0.279003    0.000000
+siesta:      7    0.005249   -0.279003    0.000000
+siesta: ----------------------------------------
+siesta:    Tot   -0.000000   -0.032764   -0.000000
+
+siesta: Stress tensor (static) (eV/Ang**3):
+siesta:    -0.002134    0.000000   -0.000000
+siesta:     0.000000   -0.000556   -0.000000
+siesta:     0.000000    0.000000   -0.000724
+
+siesta: Cell volume =        262.114919 Ang**3
+
+siesta: Pressure (static):
+siesta:                Solid            Molecule  Units
+siesta:           0.00001239          0.00000224  Ry/Bohr**3
+siesta:           0.00113783          0.00020528  eV/Ang**3
+siesta:           1.82302327          0.32889750  kBar
+(Free)E+ p_basis*V_orbitals  =        -464.513364
+(Free)Eharris+ p_basis*V_orbitals  =        -464.513364
+
+siesta: Electric dipole (a.u.)  =    0.000000    0.657893    0.000000
+siesta: Electric dipole (Debye) =    0.000000    1.672200    0.000000
+
+timer: Elapsed wall time (sec) =      11.961
+timer: CPU execution times (sec):
+
+Routine            Calls   Time/call    Tot.time        %
+siesta                 1      11.909      11.909   100.00
+Setup                  1       0.365       0.365     3.06
+bands                  1       0.000       0.000     0.00
+KSV_init               1       0.000       0.000     0.00
+IterGeom               1      11.539      11.539    96.90
+state_init             1       2.127       2.127    17.86
+hsparse                1       0.002       0.002     0.02
+overlap                1       2.123       2.123    17.83
+Setup_H0               1       7.766       7.766    65.22
+naefs                  2       0.000       0.001     0.00
+MolMec                 2       0.000       0.000     0.00
+kinefsm                2       0.990       1.979    16.62
+nlefsm                 2       2.848       5.695    47.82
+DHSCF_Init             1       0.097       0.097     0.81
+DHSCF1                 1       0.007       0.007     0.06
+INITMESH               1       0.000       0.000     0.00
+DHSCF2                 1       0.090       0.090     0.75
+REMESH                 1       0.008       0.008     0.07
+REORD                 58       0.000       0.006     0.05
+PHION                  1       0.067       0.067     0.56
+COMM_BSC              53       0.000       0.006     0.05
+POISON                26       0.007       0.174     1.46
+fft                   52       0.003       0.152     1.27
+IterSCF               24       0.061       1.464    12.29
+setup_H               24       0.059       1.420    11.92
+DHSCF                 25       0.064       1.592    13.37
+DHSCF3                25       0.059       1.475    12.38
+rhoofd                25       0.029       0.714     6.00
+cellXC                25       0.004       0.108     0.90
+vmat                  25       0.018       0.449     3.77
+compute_dm            24       0.001       0.028     0.24
+diagon                24       0.001       0.027     0.22
+r-eigvec              24       0.001       0.025     0.21
+r-buildHS             24       0.000       0.000     0.00
+rdiag                 24       0.001       0.025     0.21
+rdiag1                24       0.000       0.001     0.01
+rdiag2                24       0.000       0.002     0.02
+rdiag3                24       0.001       0.018     0.15
+rdiag4                24       0.000       0.003     0.02
+r-buildD              24       0.000       0.001     0.01
+MIXER                 23       0.000       0.000     0.00
+WriteDM               24       0.000       0.010     0.09
+PostSCF                1       0.180       0.180     1.51
+DHSCF4                 1       0.118       0.118     0.99
+dfscf                  1       0.112       0.112     0.94
+overfsm                1       0.001       0.001     0.01
+state_analysis         1       0.002       0.002     0.02
+siesta_move            1       0.000       0.000     0.00
+siesta_analysis        1       0.003       0.003     0.02
+optical                1       0.000       0.000     0.00
+  
+>> End of run:  17-JAN-2019  15:47:57
+Job completed

=== added directory 'Tests/bessel-rich'
=== added file 'Tests/bessel-rich/bessel-rich.fdf'
--- Tests/bessel-rich/bessel-rich.fdf	1970-01-01 00:00:00 +0000
+++ Tests/bessel-rich/bessel-rich.fdf	2019-01-30 13:27:17 +0000
@@ -0,0 +1,45 @@
+#
+# Needs new feature: handling of fewer rc's than nzetas in PAO.Basis block
+#
+write-ion-plot-files T
+#
+SystemName          Water molecule with various Bessel Orbitals
+SystemLabel         bessel-rich
+NumberofAtoms   7
+NumberOfSpecies 4
+%block ChemicalSpeciesLabel
+ 1  8  O      # Species index, atomic number, species label
+ 2  1  H
+ 3  -100 Bessel
+ 4  -100 J
+%endblock ChemicalSpeciesLabel
+
+AtomicCoordinatesFormat  Ang
+%block AtomicCoordinatesAndAtomicSpecies
+ 0.000  0.000  0.000  1
+ 0.757  0.586  0.000  2
+-0.757  0.586  0.000  2
+ 0.3785 0.293  0.000  3
+-0.3785 0.293  0.000  3
+ 0.3785 0.293  0.000  4
+-0.3785 0.293  0.000  4
+%endblock AtomicCoordinatesAndAtomicSpecies
+
+%block PAO.Basis
+Bessel    3
+ n=1    0    1   
+   2.0
+   1.0
+ n=2    0    1   
+   2.5
+   1.0
+ n=3    1    1   
+   3.5
+   1.0
+J         2  # l-shells
+n=2   0   7    # Note new feature: fewer rc's than zetas
+  4.5 
+n=2   1   3
+  4.5 4.5 5.0
+%endblock PAO.Basis
+  
\ No newline at end of file

=== added file 'Tests/bessel-rich/bessel-rich.pseudos'
--- Tests/bessel-rich/bessel-rich.pseudos	1970-01-01 00:00:00 +0000
+++ Tests/bessel-rich/bessel-rich.pseudos	2019-01-30 13:27:17 +0000
@@ -0,0 +1,1 @@
+H O

=== added file 'Tests/bessel-rich/makefile'
--- Tests/bessel-rich/makefile	1970-01-01 00:00:00 +0000
+++ Tests/bessel-rich/makefile	2019-01-30 13:27:17 +0000
@@ -0,0 +1,6 @@
+#
+# Single-test makefile
+#
+name=bessel-rich
+#
+include ../test.mk

=== modified file 'version.info'
--- version.info	2018-12-19 13:52:34 +0000
+++ version.info	2019-01-30 13:27:17 +0000
@@ -1,1 +1,4 @@
-siesta-4.0--589
+siesta-4.0--589--n-fix-5--semic-5
+
+
+