From cb63df029cc8bb00da2884235fd374d3a635e4db Mon Sep 17 00:00:00 2001
From: Jerome Jackson <jerome@neuroplug.net>
Date: Wed, 22 May 2024 16:29:57 +0100
Subject: [PATCH] cleanup set_atoms and move projector read to _read_special()

(and allow projector definition to be skipped)
---
 src/c_interface.F90               |   2 +-
 src/library_interface.F90         |  41 ++-
 src/postw90/postw90_readwrite.F90 |   6 +-
 src/readwrite.F90                 | 550 ++++++++++++++++++------------
 src/types.F90                     |   1 +
 src/wannier90_readwrite.F90       |  91 ++---
 6 files changed, 389 insertions(+), 302 deletions(-)

diff --git a/src/c_interface.F90 b/src/c_interface.F90
index f217d887c..b3dc44f70 100644
--- a/src/c_interface.F90
+++ b/src/c_interface.F90
@@ -75,7 +75,7 @@ subroutine cinput_setopt(common_cptr, seedname, ierr, comm) bind(c)
     call w90_get_fortran_stderr(istderr)
     call w90_get_fortran_stdout(istdout)
     call c_f_pointer(common_cptr, common_fptr)
-    call w90_input_setopt(common_fptr, seedname, comm, istdout, istderr, ierr)
+    call w90_input_setopt(common_fptr, seedname, istdout, istderr, ierr)
   end subroutine cinput_setopt
 
   subroutine cinput_reader(common_cptr, ierr) bind(c)
diff --git a/src/library_interface.F90 b/src/library_interface.F90
index 758df1d81..5b3b532d2 100644
--- a/src/library_interface.F90
+++ b/src/library_interface.F90
@@ -183,6 +183,7 @@ module w90_library
     module procedure w90_set_option_i2d
     module procedure w90_set_option_r1d
     module procedure w90_set_option_r2d
+    module procedure w90_set_option_c2d
     module procedure w90_set_option_real
   end interface w90_set_option
 
@@ -245,6 +246,8 @@ subroutine w90_input_setopt(common_data, seedname, istdout, istderr, ierr)
       return
     endif
 
+    ! input_setopt() processes options stored in a container, common_data%settings
+    ! meanwhile the .win file contents are stored in common_data%in_data, which should be empty here
     if (allocated(common_data%settings%in_data)) then
       call set_error_fatal(error, ' readinput and setopt clash at input_setopt call', common_data%comm)
       call prterr(error, ierr, istdout, istderr, common_data%comm)
@@ -338,7 +341,7 @@ subroutine w90_input_setopt(common_data, seedname, istdout, istderr, ierr)
       return
     endif
 
-    ! clear any settings (from settings interface not .win file)
+    ! clear settings container (from settings interface not .win file)
     deallocate (common_data%settings%entries, stat=ierr)
     if (ierr /= 0) then
       call set_error_alloc(error, 'Error in deallocating entries data in input_setopt() library call', common_data%comm)
@@ -437,6 +440,7 @@ subroutine w90_disentangle(common_data, istdout, istderr, ierr)
 
     ! local variables
     type(w90_error_type), allocatable :: error
+    integer :: ioff
 
     ierr = 0
 
@@ -463,8 +467,15 @@ subroutine w90_disentangle(common_data, istdout, istderr, ierr)
     if (common_data%dis_manifold%win_min == -huge(0.0_dp)) common_data%dis_manifold%win_min = minval(common_data%eigval)
     if (common_data%dis_manifold%win_max == huge(0.0_dp)) common_data%dis_manifold%win_max = maxval(common_data%eigval)
     if (common_data%dis_manifold%frozen_states) then
-      if (common_data%dis_manifold%froz_min == -huge(0.0_dp)) &
-        common_data%dis_manifold%froz_min = common_data%dis_manifold%win_min
+      !write(*,*) shape(common_data%eigval)
+      !write(*,*) shape(common_data%dis_manifold%lwindow)
+      if (common_data%dis_manifold%froz_min == -huge(0.0_dp)) then
+        ioff = maxval(common_data%exclude_bands)
+        common_data%dis_manifold%froz_min = minval(common_data%eigval(ioff + 1, :))
+      endif
+      !common_data%dis_manifold%froz_min = minval(common_data%eigval, mask=common_data%dis_manifold%lwindow)
+      !write(*,*)common_data%dis_manifold%froz_min,allocated(common_data%dis_manifold%lwindow),&
+      !  count(common_data%dis_manifold%lwindow)
     endif
 
     if (common_data%num_bands > common_data%num_wann) then
@@ -879,6 +890,24 @@ subroutine w90_set_option_r2d(common_data, keyword, arr)
     if (common_data%settings%num_entries == common_data%settings%num_entries_max) call expand_settings(common_data%settings)
   endsubroutine w90_set_option_r2d
 
+  subroutine w90_set_option_c2d(common_data, keyword, arr)
+    use w90_readwrite, only: init_settings, expand_settings
+
+    implicit none
+
+    character(*), intent(in) :: keyword
+    character(len=*), intent(in) :: arr(:)
+    type(lib_common_type), intent(inout) :: common_data
+    integer :: i
+
+    if (.not. allocated(common_data%settings%entries)) call init_settings(common_data%settings)
+    i = common_data%settings%num_entries + 1
+    common_data%settings%entries(i)%keyword = keyword
+    common_data%settings%entries(i)%c2d = arr
+    common_data%settings%num_entries = i + 1
+    if (common_data%settings%num_entries == common_data%settings%num_entries_max) call expand_settings(common_data%settings)
+  endsubroutine w90_set_option_c2d
+
   subroutine w90_set_option_real(common_data, keyword, rval)
     use w90_readwrite, only: init_settings, expand_settings
 
@@ -1009,7 +1038,7 @@ subroutine w90_get_proj(common_data, n, site, l, m, s, rad, x, z, sqa, istdout,
     integer, intent(in) :: istdout, istderr
 
     ! probably the remaining variables find limited use, allow them to be absent
-    integer, intent(inout), optional :: rad(:)
+    integer, intent(inout) :: rad(:)
     real(kind=dp), intent(inout) :: sqa(:, :), z(:, :), x(:, :) !, zona(:)
 
     ! local variables
@@ -1027,6 +1056,8 @@ subroutine w90_get_proj(common_data, n, site, l, m, s, rad, x, z, sqa, istdout,
 
     n = size(common_data%proj_input)
 
+    write (*, *) "expecting: ", n
+
     ! check allocation of main output arrays
     if (size(l) < n) then
       call set_error_fatal(error, 'Array argument l in get_proj() call is insufficiently sized', common_data%comm)
@@ -1089,6 +1120,8 @@ subroutine w90_get_proj(common_data, n, site, l, m, s, rad, x, z, sqa, istdout,
       z(:, ip) = proj%z(:)
       x(:, ip) = proj%x(:)
       rad(ip) = proj%radial
+
+      write (*, *) " projector   : ", l(ip), m(ip), s(ip)
       !if (present(zona)) zona(ip) = proj%zona
     enddo
   end subroutine w90_get_proj
diff --git a/src/postw90/postw90_readwrite.F90 b/src/postw90/postw90_readwrite.F90
index 474c0f1ac..70e58730c 100644
--- a/src/postw90/postw90_readwrite.F90
+++ b/src/postw90/postw90_readwrite.F90
@@ -168,8 +168,7 @@ subroutine w90_postw90_readwrite_read(settings, ws_region, w90_system, exclude_b
     if (allocated(error)) return
     call w90_readwrite_read_exclude_bands(settings, exclude_bands, num_exclude_bands, error, comm)
     if (allocated(error)) return
-    call w90_readwrite_read_num_bands(settings, effective_model, num_bands, num_wann, stdout, &
-                                      error, comm)
+    call w90_readwrite_read_num_bands(settings, effective_model, num_bands, num_wann, error, comm)
     if (allocated(error)) return
     disentanglement = (num_bands > num_wann)
     call w90_readwrite_read_mp_grid(settings, effective_model, mp_grid, num_kpts, error, comm)
@@ -224,7 +223,8 @@ subroutine w90_postw90_readwrite_read(settings, ws_region, w90_system, exclude_b
     dis_manifold%win_max = 0.0_dp
     if (eig_found) dis_manifold%win_min = minval(eigval)
     if (eig_found) dis_manifold%win_max = maxval(eigval)
-    call w90_readwrite_read_dis_manifold(settings, eig_found, dis_manifold, error, comm)
+    call w90_readwrite_read_dis_manifold(settings, dis_manifold, error, comm)
+
     if (allocated(error)) return
     call w90_wannier90_readwrite_read_geninterp(settings, pw90_geninterp, error, comm)
     if (allocated(error)) return
diff --git a/src/readwrite.F90 b/src/readwrite.F90
index 2a2d793e2..4a143b3f4 100644
--- a/src/readwrite.F90
+++ b/src/readwrite.F90
@@ -19,7 +19,6 @@
 !------------------------------------------------------------!
 
 module w90_readwrite
-
   !! Common read/write routines for data needed by both
   !! wannier90.x and postw90.x executables
 
@@ -31,29 +30,26 @@ module w90_readwrite
 
   private
 
-  public :: w90_readwrite_read_distk
   public :: w90_readwrite_chkpt_dist
+  public :: w90_readwrite_clean_infile
+  public :: w90_readwrite_clear_keywords
   public :: w90_readwrite_dealloc
-  public :: w90_readwrite_get_convention_type
-  public :: w90_readwrite_get_smearing_type
-  public :: w90_readwrite_lib_set_atoms
-  public :: w90_readwrite_read_chkpt
-  public :: w90_readwrite_read_chkpt_header
-  public :: w90_readwrite_read_chkpt_matrices
-  public :: w90_readwrite_write_header
   public :: w90_readwrite_get_block_length
   public :: w90_readwrite_get_centre_constraints
+  public :: w90_readwrite_get_convention_type
   public :: w90_readwrite_get_projections
   public :: w90_readwrite_get_range_vector
   public :: w90_readwrite_get_smearing_index
+  public :: w90_readwrite_get_smearing_type
   public :: w90_readwrite_get_vector_length
   public :: w90_readwrite_in_file
-  public :: w90_readwrite_set_kmesh
-  public :: w90_readwrite_clean_infile
-  public :: w90_readwrite_clear_keywords
   public :: w90_readwrite_read_algorithm_control
   public :: w90_readwrite_read_atoms
+  public :: w90_readwrite_read_chkpt
+  public :: w90_readwrite_read_chkpt_header
+  public :: w90_readwrite_read_chkpt_matrices
   public :: w90_readwrite_read_dis_manifold
+  public :: w90_readwrite_read_distk
   public :: w90_readwrite_read_eigvals
   public :: w90_readwrite_read_exclude_bands
   public :: w90_readwrite_read_fermi_energy
@@ -67,11 +63,15 @@ module w90_readwrite
   public :: w90_readwrite_read_mp_grid
   public :: w90_readwrite_read_num_bands
   public :: w90_readwrite_read_num_wann
-  public :: w90_readwrite_read_total_bands
   public :: w90_readwrite_read_system
+  public :: w90_readwrite_read_total_bands
   public :: w90_readwrite_read_units
   public :: w90_readwrite_read_verbosity
   public :: w90_readwrite_read_ws_data
+  public :: w90_readwrite_set_kmesh
+  public :: w90_readwrite_write_header
+
+  private :: w90_readwrite_set_atoms
 
   public :: w90_readwrite_get_keyword
   public :: w90_readwrite_get_keyword_block
@@ -83,17 +83,21 @@ module w90_readwrite
 contains
   !================================================!
   subroutine w90_readwrite_read_verbosity(settings, print_output, svd_omega, error, comm)
+    !! read verbosity "iprint" and timing "timing_level" variables
+    !! if iprint>2 svd_omega printing is enabled
+    !! printing is supressed on all non-root MPI ranks
     use w90_error, only: w90_error_type
     use w90_comms, only: mpirank
     implicit none
-    type(print_output_type), intent(inout) :: print_output
+
     logical, intent(inout) :: svd_omega
-    logical :: found
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+    type(print_output_type), intent(inout) :: print_output
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: unlucky_rank
+    logical :: found
 
     call w90_readwrite_get_keyword(settings, 'timing_level', found, error, comm, &
                                    i_value=print_output%timing_level)
@@ -118,6 +122,7 @@ subroutine w90_readwrite_read_verbosity(settings, print_output, svd_omega, error
   end subroutine w90_readwrite_read_verbosity
 
   subroutine w90_readwrite_read_algorithm_control(settings, optimisation, error, comm)
+    !! reads the "optimisation" flag
     use w90_error, only: w90_error_type
     implicit none
     integer, intent(inout) :: optimisation
@@ -134,15 +139,16 @@ end subroutine w90_readwrite_read_algorithm_control
 
   subroutine w90_readwrite_read_units(settings, lenconfac, length_unit, energy_unit, bohr, error, &
                                       comm)
+    !! reads the "energy_unit" and "length_unit" (valid: "ang" or "bohr") variables
     use w90_error, only: w90_error_type, set_error_input
     implicit none
-    real(kind=dp), intent(inout) :: lenconfac
-    character(len=*), intent(inout) :: length_unit
     character(len=*), intent(inout) :: energy_unit
+    character(len=*), intent(inout) :: length_unit
     real(kind=dp), intent(in) :: bohr
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+    real(kind=dp), intent(inout) :: lenconfac
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: ic
     logical :: found
@@ -152,24 +158,24 @@ subroutine w90_readwrite_read_units(settings, lenconfac, length_unit, energy_uni
 
     call w90_readwrite_get_keyword(settings, 'length_unit', found, error, comm, c_value=length_unit)
     if (allocated(error)) return
-    if (found) then ! if not specified, it may take already have capitalised values, Ang or Bohr
+    if (found) then
       if (length_unit .ne. 'ang' .and. length_unit .ne. 'bohr') then
-        call set_error_input(error, 'Error: value of length_unit not recognised in w90_readwrite_read_units', comm)
+        call set_error_input(error, &
+                             'Error: value of length_unit not recognised in w90_readwrite_read_units', comm)
         return
       else if (length_unit .eq. 'bohr') then
         lenconfac = 1.0_dp/bohr
       endif
     endif
 
-    ! Length unit (ang --> Ang, bohr --> Bohr)
-    ! this is set to uppercase only for printout...
+    ! Length unit (ang --> Ang, bohr --> Bohr) set to uppercase for printout
     ic = ichar(length_unit(1:1))
     if ((ic .ge. ichar('a')) .and. (ic .le. ichar('z'))) &
       length_unit(1:1) = char(ic + ichar('Z') - ichar('z'))
-
   end subroutine w90_readwrite_read_units
 
   subroutine w90_readwrite_read_num_wann(settings, num_wann, error, comm)
+    !! reads the number of wannier functions "num_wann" (mandatory input)
     use w90_error, only: w90_error_type, set_error_input
     implicit none
     integer, intent(inout) :: num_wann
@@ -185,48 +191,41 @@ subroutine w90_readwrite_read_num_wann(settings, num_wann, error, comm)
     if (.not. found) then
       call set_error_input(error, 'Error: You must specify num_wann', comm)
       return
-    endif
-    if (num_wann <= 0) then
+    else if (num_wann <= 0) then
       call set_error_input(error, 'Error: num_wann must be greater than zero', comm)
       return
     endif
   end subroutine w90_readwrite_read_num_wann
 
-  ! fixme
-  subroutine w90_readwrite_read_total_bands(settings, num_wann, error, comm)
+  subroutine w90_readwrite_read_total_bands(settings, total_bands, error, comm)
+    !! read the "total_bands" variable
+    !! this is a convenience for combination with an "exclude_bands" to evaluate num_bands
     use w90_error, only: w90_error_type, set_error_input
     implicit none
-    integer, intent(inout) :: num_wann
+    integer, intent(inout) :: total_bands
     type(w90_error_type), allocatable, intent(out) :: error
     type(w90_comm_type), intent(in) :: comm
     type(settings_type), intent(inout) :: settings
 
     logical :: found
-    num_wann = 0
+    total_bands = 0
 
-    call w90_readwrite_get_keyword(settings, 'total_bands', found, error, comm, i_value=num_wann)
+    call w90_readwrite_get_keyword(settings, 'total_bands', found, error, comm, i_value=total_bands)
     if (allocated(error)) return
-
-    !if (.not. found) then
-    !  call set_error_input(error, 'Error: You must specify num_wann', comm)
-    !  return
-    !endif
-    if (num_wann <= 0) then
-      call set_error_input(error, 'Error: num_wann must be greater than zero', comm)
-      return
-    endif
   end subroutine w90_readwrite_read_total_bands
 
   subroutine w90_readwrite_read_distk(settings, distk, nkin, error, comm)
-    ! read distribution of kpoints
+    !! Read MPI distribution of k-points
+    !! The array to be read must have num_kpt entries, with each entry being
+    !! the MPI rank to which each k-point is assigned
     use w90_error, only: w90_error_type, set_error_input, set_error_alloc
     implicit none
 
     integer, allocatable, intent(inout) :: distk(:)
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
-    type(settings_type), intent(inout) :: settings
     integer, intent(in) :: nkin
+    type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: nk
     logical :: found
@@ -238,31 +237,30 @@ subroutine w90_readwrite_read_distk(settings, distk, nkin, error, comm)
 
     if (found) then
       if (nk /= nkin) then
-        call set_error_input(error, 'incorrect length of k-distribution', comm)
+        call set_error_input(error, 'Error: incorrect length of k-distribution (distk)', comm)
         return
       endif
-
       allocate (distk(nkin))
-
       call w90_readwrite_get_range_vector(settings, 'distk', found, nk, .false., error, comm, distk)
       if (allocated(error)) return
     else
+      ! fixme JJ, some output might be helpful here
       allocate (distk(nkin))
-
       distk = 0 ! default to no distribution if not specified
     end if
   endsubroutine w90_readwrite_read_distk
 
   subroutine w90_readwrite_read_exclude_bands(settings, exclude_bands, num_exclude_bands, error, &
                                               comm)
+    !! Read (and allocate) excluded_bands list "exclude_bands"
     use w90_error, only: w90_error_type, set_error_input, set_error_alloc
     implicit none
 
     integer, allocatable, intent(inout) :: exclude_bands(:)
     integer, intent(out) :: num_exclude_bands
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: ierr
     logical :: found = .false.
@@ -277,10 +275,10 @@ subroutine w90_readwrite_read_exclude_bands(settings, exclude_bands, num_exclude
         call set_error_input(error, 'Error: problem reading exclude_bands', comm)
         return
       endif
-      if (allocated(exclude_bands)) deallocate (exclude_bands)
       allocate (exclude_bands(num_exclude_bands), stat=ierr)
       if (ierr /= 0) then
-        call set_error_alloc(error, 'Error allocating exclude_bands in w90_readwrite_read_exclude_bands', comm)
+        call set_error_alloc(error, &
+                             'Error allocating exclude_bands in w90_readwrite_read_exclude_bands', comm)
         return
       endif
       call w90_readwrite_get_range_vector(settings, 'exclude_bands', found, num_exclude_bands, &
@@ -294,16 +292,18 @@ subroutine w90_readwrite_read_exclude_bands(settings, exclude_bands, num_exclude
   end subroutine w90_readwrite_read_exclude_bands
 
   subroutine w90_readwrite_read_num_bands(settings, pw90_effective_model, num_bands, num_wann, &
-                                          stdout, error, comm)
+                                          error, comm)
+    !! Read the number of bands ("num_bands")
+    !! If not specified (and exclude_bands and total_bands are not both provided), defaults to num_wann
     use w90_error, only: w90_error_type, set_error_input
     implicit none
-    logical, intent(in) :: pw90_effective_model
-    integer, intent(inout) :: num_bands
+
     integer, intent(in) :: num_wann
-    integer, intent(in) :: stdout
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+    integer, intent(inout) :: num_bands
+    logical, intent(in) :: pw90_effective_model
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: i_temp
     logical :: found
@@ -316,9 +316,8 @@ subroutine w90_readwrite_read_num_bands(settings, pw90_effective_model, num_band
       if (found) then
         num_bands = i_temp
         if (num_bands < num_wann) then
-          write (stdout, *) 'num_bands', num_bands
-          write (stdout, *) 'num_wann', num_wann
-          call set_error_input(error, 'Error: num_bands must be greater than or equal to num_wann', comm)
+          call set_error_input(error, &
+                               'Error: num_bands must be greater than or equal to num_wann', comm)
           return
         endif
       else
@@ -328,20 +327,26 @@ subroutine w90_readwrite_read_num_bands(settings, pw90_effective_model, num_band
   end subroutine w90_readwrite_read_num_bands
 
   subroutine w90_readwrite_read_gamma_only(settings, gamma_only, num_kpts, error, comm)
+    !! Reads the flag for Gamma-only mode ("gamma_only")
     use w90_error, only: w90_error_type, set_error_input
     implicit none
-    logical, intent(inout) :: gamma_only
+
     integer, intent(in) :: num_kpts
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+    logical, intent(inout) :: gamma_only
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     logical :: found, ltmp
 
     ltmp = .false.
+    gamma_only = .false.
+
     call w90_readwrite_get_keyword(settings, 'gamma_only', found, error, comm, l_value=ltmp)
     if (allocated(error)) return
-    gamma_only = ltmp
+
+    if (found) gamma_only = ltmp
+
     if (gamma_only .and. (num_kpts .ne. 1)) then
       call set_error_input(error, 'Error: gamma_only is true, but num_kpts > 1', comm)
       return
@@ -350,13 +355,15 @@ end subroutine w90_readwrite_read_gamma_only
 
   subroutine w90_readwrite_read_mp_grid(settings, pw90_effective_model, mp_grid, num_kpts, error, &
                                         comm)
+    !! Read the mandatory k-point mesh input ("mp_grid")
     use w90_error, only: w90_error_type, set_error_input
     implicit none
-    logical, intent(in) :: pw90_effective_model
+
     integer, intent(inout) :: mp_grid(3), num_kpts
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+    logical, intent(in) :: pw90_effective_model
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: iv_temp(3)
     logical :: found
@@ -369,7 +376,8 @@ subroutine w90_readwrite_read_mp_grid(settings, pw90_effective_model, mp_grid, n
     if (.not. pw90_effective_model) then
       if (found) mp_grid = iv_temp
       if (.not. found) then
-        call set_error_input(error, 'Error: You must specify dimensions of the Monkhorst-Pack grid by setting mp_grid', comm)
+        call set_error_input(error, &
+                             'Error: You must specify dimensions of the Monkhorst-Pack grid by setting mp_grid', comm)
         return
       elseif (any(mp_grid < 1)) then
         call set_error_input(error, 'Error: mp_grid must be greater than zero', comm)
@@ -380,12 +388,17 @@ subroutine w90_readwrite_read_mp_grid(settings, pw90_effective_model, mp_grid, n
   end subroutine w90_readwrite_read_mp_grid
 
   subroutine w90_readwrite_read_system(settings, w90_system, error, comm)
+    !! Read a group of variables defining the system
+    !! "spinors" -- coupled spins
+    !! "num_elec_per_state" -- spin degeneracy
+    !! "num_valence_bands"
     use w90_error, only: w90_error_type, set_error_input
     implicit none
-    type(w90_system_type), intent(inout) :: w90_system
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
+    type(w90_system_type), intent(inout) :: w90_system
 
     logical :: found, ltmp
     integer :: itmp
@@ -396,10 +409,9 @@ subroutine w90_readwrite_read_system(settings, w90_system, error, comm)
     if (found) then
       w90_system%spinors = ltmp
     else
-      w90_system%spinors = .false. ! specify a default behaviour
+      w90_system%spinors = .false.
     endif
-    ! We need to know if the bands are double degenerate due to spin, e.g. when
-    ! calculating the DOS
+    ! We need to know if the bands are double degenerate due to spin, e.g. when calculating DOS
     if (w90_system%spinors) then
       w90_system%num_elec_per_state = 1
     else
@@ -432,24 +444,27 @@ subroutine w90_readwrite_read_system(settings, w90_system, error, comm)
     endif
   end subroutine w90_readwrite_read_system
 
-  subroutine w90_readwrite_read_kpath(settings, kpoint_path, ok, bands_plot, error, comm)
+  subroutine w90_readwrite_read_kpath(settings, kpoint_path, path_found, bands_plot, error, comm)
+    !! Read band plotting path variables: "kpoint_path" and "bands_num_points"
     use w90_error, only: w90_error_type, set_error_input, set_error_alloc
     implicit none
+
     logical, intent(in) :: bands_plot
+    logical, intent(out) :: path_found
     type(kpoint_path_type), intent(inout) :: kpoint_path
-    logical, intent(out) :: ok
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
     integer :: i_temp, ierr, bands_num_spec_points
     logical :: found
 
+    path_found = .false.
     bands_num_spec_points = 0
-    call w90_readwrite_get_block_length(settings, 'kpoint_path', found, i_temp, error, comm)
+
+    call w90_readwrite_get_block_length(settings, 'kpoint_path', path_found, i_temp, error, comm)
     if (allocated(error)) return
-    if (found) then
-      ok = .true.
+    if (path_found) then
       bands_num_spec_points = i_temp*2
       if (allocated(kpoint_path%labels)) deallocate (kpoint_path%labels)
       allocate (kpoint_path%labels(bands_num_spec_points), stat=ierr)
@@ -465,9 +480,8 @@ subroutine w90_readwrite_read_kpath(settings, kpoint_path, ok, bands_plot, error
       endif
       call w90_readwrite_get_keyword_kpath(settings, kpoint_path, error, comm)
       if (allocated(error)) return
-    else
-      ok = .false.
     end if
+
     call w90_readwrite_get_keyword(settings, 'bands_num_points', found, error, comm, &
                                    i_value=kpoint_path%num_points_first_segment)
     if (allocated(error)) return
@@ -531,91 +545,91 @@ end subroutine w90_readwrite_read_explicit_kpath
 
   subroutine w90_readwrite_read_fermi_energy(settings, found_fermi_energy, fermi_energy_list, &
                                              error, comm)
+    !! Read Fermi energy ("fermi_energy") and/or ranges ("fermi_energy_min", "fermi_energy_max" and
+    !! "fermi_energy_step") used to setup fermi_energy_list tabulation
+    !! _max and _step are only sought if _min found and are optional
     use w90_error, only: w90_error_type, set_error_input, set_error_alloc
     implicit none
-    logical, intent(out) :: found_fermi_energy
+
+    ! arguments
+    logical, intent(out) :: found_fermi_energy ! flags that E_F provided
     real(kind=dp), allocatable, intent(out) :: fermi_energy_list(:)
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
+    ! local variables
+    integer :: i, ierr, n
+    logical :: found, fermi_energy_scan
     real(kind=dp) :: fermi_energy
-    logical :: fermi_energy_scan
-    real(kind=dp) :: fermi_energy_min
     real(kind=dp) :: fermi_energy_max
+    real(kind=dp) :: fermi_energy_min
     real(kind=dp) :: fermi_energy_step
-    integer :: i, ierr, n
-    logical :: found
 
-    n = 0
     found_fermi_energy = .false.
-    call w90_readwrite_get_keyword(settings, 'fermi_energy', found, error, comm, &
+    fermi_energy_scan = .false.
+    n = 1
+    fermi_energy = 0.0_dp
+    fermi_energy_min = fermi_energy
+    fermi_energy_step = 0.0_dp
+
+    call w90_readwrite_get_keyword(settings, 'fermi_energy', found_fermi_energy, error, comm, &
                                    r_value=fermi_energy)
     if (allocated(error)) return
-    if (found) then
-      found_fermi_energy = .true.
+    if (found_fermi_energy) then
       n = 1
+      fermi_energy_step = 0.0_dp
+      fermi_energy_min = fermi_energy
     endif
 
-    fermi_energy_scan = .false.
-    call w90_readwrite_get_keyword(settings, 'fermi_energy_min', found, error, comm, &
+    call w90_readwrite_get_keyword(settings, 'fermi_energy_min', fermi_energy_scan, error, comm, &
                                    r_value=fermi_energy_min)
     if (allocated(error)) return
-    if (found) then
+    if (fermi_energy_scan) then
       if (found_fermi_energy) then
-        call set_error_input(error, 'Error: Cannot specify both fermi_energy and fermi_energy_min', comm)
+        call set_error_input(error, &
+                             'Error: Cannot specify both fermi_energy and fermi_energy_min', comm)
         return
       endif
-      fermi_energy_scan = .true.
-      fermi_energy_max = fermi_energy_min + 1.0_dp
+
       call w90_readwrite_get_keyword(settings, 'fermi_energy_max', found, error, comm, &
                                      r_value=fermi_energy_max)
       if (allocated(error)) return
-      if (found .and. fermi_energy_max <= fermi_energy_min) then
-        call set_error_input(error, 'Error: fermi_energy_max must be larger than fermi_energy_min', comm)
+      if (.not. found) then
+        fermi_energy_max = fermi_energy_min + 1.0_dp !default
+      else if (found .and. fermi_energy_max <= fermi_energy_min) then
+        call set_error_input(error, &
+                             'Error: fermi_energy_max must be larger than fermi_energy_min', comm)
         return
       endif
-      fermi_energy_step = 0.01_dp
+
       call w90_readwrite_get_keyword(settings, 'fermi_energy_step', found, error, comm, &
                                      r_value=fermi_energy_step)
       if (allocated(error)) return
-      if (found .and. fermi_energy_step <= 0.0_dp) then
+      if (.not. found) then
+        fermi_energy_step = 0.01_dp !default
+      else if (found .and. fermi_energy_step <= 0.0_dp) then
         call set_error_input(error, 'Error: fermi_energy_step must be positive', comm)
         return
       endif
+
       n = nint(abs((fermi_energy_max - fermi_energy_min)/fermi_energy_step)) + 1
+      fermi_energy_step = (fermi_energy_max - fermi_energy_min)/real(n - 1, dp)
     endif
 
-    if (found_fermi_energy) then
-      if (allocated(fermi_energy_list)) deallocate (fermi_energy_list)
-      allocate (fermi_energy_list(1), stat=ierr)
-      fermi_energy_list(1) = fermi_energy
-    elseif (fermi_energy_scan) then
-      if (n .eq. 1) then
-        fermi_energy_step = 0.0_dp
-      else
-        fermi_energy_step = (fermi_energy_max - fermi_energy_min)/real(n - 1, dp)
-      endif
-      if (allocated(fermi_energy_list)) deallocate (fermi_energy_list)
-      allocate (fermi_energy_list(n), stat=ierr)
-      do i = 1, n
-        fermi_energy_list(i) = fermi_energy_min + (i - 1)*fermi_energy_step
-      enddo
-!! AAM_2017-03-27: if fermi_energy* parameters are not set in input file
-!! then allocate fermi_energy_list with length 1 and set to zero as default.
-    else
-      if (allocated(fermi_energy_list)) deallocate (fermi_energy_list)
-      allocate (fermi_energy_list(1), stat=ierr)
-      fermi_energy_list(1) = 0.0_dp
-    endif
+    allocate (fermi_energy_list(n), stat=ierr)
     if (ierr /= 0) then
       call set_error_alloc(error, &
                            'Error allocating fermi_energy_list in w90_readwrite_read_fermi_energy', comm)
       return
     endif
+    do i = 1, n
+      fermi_energy_list(i) = fermi_energy_min + (i - 1)*fermi_energy_step
+    enddo
   end subroutine w90_readwrite_read_fermi_energy
 
   subroutine w90_readwrite_read_ws_data(settings, ws_region, error, comm)
+    !! Reads "use_ws_distance", "ws_distance_tol", "ws_search_size"
     use w90_error, only: w90_error_type, set_error_input
     implicit none
     type(ws_region_type), intent(inout) :: ws_region
@@ -662,11 +676,12 @@ end subroutine w90_readwrite_read_ws_data
   subroutine w90_readwrite_read_eigvals(eig_found, eigval, num_bands, num_kpts, stdout, &
                                         seedname, error, comm)
     !! Read the eigenvalues from wannier.eig
-
+    ! fixme (jj) consider relocating to library_extra
     use w90_error, only: w90_error_type, set_error_file, set_error_file, set_error_alloc
 
     implicit none
 
+    ! arguments
     character(len=*), intent(in)  :: seedname
     integer, intent(in) :: num_bands, num_kpts
     integer, intent(in) :: stdout
@@ -713,12 +728,13 @@ subroutine w90_readwrite_read_eigvals(eig_found, eigval, num_bands, num_kpts, st
     return
   end subroutine w90_readwrite_read_eigvals
 
-  subroutine w90_readwrite_read_dis_manifold(settings, eig_found, dis_manifold, error, comm)
+  subroutine w90_readwrite_read_dis_manifold(settings, dis_manifold, error, comm)
+    !! Reads disentanglement windows "dis_win_min" and "dis_win_max" (both are optional)
+    !! Reads frozen window "dis_froz_min" and "dis_froz_max" (either neither or both to be supplied)
     use w90_error, only: w90_error_type, set_error_input
     implicit none
 
     ! arguments
-    logical, intent(in) :: eig_found
     type(dis_manifold_type), intent(inout) :: dis_manifold
     type(w90_error_type), allocatable, intent(out) :: error
     type(w90_comm_type), intent(in) :: comm
@@ -734,8 +750,8 @@ subroutine w90_readwrite_read_dis_manifold(settings, eig_found, dis_manifold, er
     call w90_readwrite_get_keyword(settings, 'dis_win_max', found, error, comm, &
                                    r_value=dis_manifold%win_max)
     if (allocated(error)) return
-    ! eig_found check because eigenvalue reading may reset win_min,max limits
-    if (eig_found .and. (dis_manifold%win_max .lt. dis_manifold%win_min)) then
+
+    if (dis_manifold%win_max .lt. dis_manifold%win_min) then
       call set_error_input(error, &
                            'Error: w90_readwrite_read_dis_manifold: check disentanglement windows (win_max < win_min !)', comm)
       return
@@ -744,21 +760,20 @@ subroutine w90_readwrite_read_dis_manifold(settings, eig_found, dis_manifold, er
     call w90_readwrite_get_keyword(settings, 'dis_froz_max', found, error, comm, &
                                    r_value=dis_manifold%froz_max)
     if (allocated(error)) return
-    if (found) then
-      dis_manifold%frozen_states = .true.
-    end if
+
+    if (found) dis_manifold%frozen_states = .true.
+
     call w90_readwrite_get_keyword(settings, 'dis_froz_min', found2, error, comm, &
                                    r_value=dis_manifold%froz_min)
     if (allocated(error)) return
-    if (eig_found) then
-      if (dis_manifold%froz_max .lt. dis_manifold%froz_min) then
-        call set_error_input(error, 'Error: w90_readwrite_read_dis_manifold: check disentanglement frozen windows', comm)
-        return
-      endif
-      if (found2 .and. .not. found) then
-        call set_error_input(error, 'Error: w90_readwrite_read_dis_manifold: found dis_froz_min but not dis_froz_max', comm)
-        return
-      endif
+
+    if (dis_manifold%froz_max .lt. dis_manifold%froz_min) then
+      call set_error_input(error, 'Error: w90_readwrite_read_dis_manifold: check disentanglement frozen windows', comm)
+      return
+    endif
+    if (found2 .and. .not. found) then
+      call set_error_input(error, 'Error: w90_readwrite_read_dis_manifold: found dis_froz_min but not dis_froz_max', comm)
+      return
     endif
 
     ! ndimwin/lwindow are not read
@@ -799,6 +814,12 @@ subroutine w90_readwrite_read_dis_manifold(settings, eig_found, dis_manifold, er
   end subroutine w90_readwrite_read_dis_manifold
 
   subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
+    !! Reads finite-difference input variables:
+    !!   "search_shells"
+    !!   "kmesh_tol"
+    !!   "shell_list"
+    !!   "num_shells"
+    !!   "skip_B1_tests"
     use w90_error, only: w90_error_type, set_error_input, set_error_alloc
     implicit none
     type(kmesh_input_type), intent(inout) :: kmesh_input
@@ -860,7 +881,6 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
         call set_error_input(error, 'Error: number of shell in shell_list must be between zero and kmesh_input%max_shells_h', comm)
         return
       endif
-      !if (allocated(kmesh_input%shell_list)) deallocate (kmesh_input%shell_list)
       allocate (kmesh_input%shell_list(kmesh_input%num_shells), stat=ierr)
       if (ierr /= 0) then
         call set_error_alloc(error, 'Error allocating shell_list in w90_wannier90_readwrite_read', comm)
@@ -874,7 +894,6 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
         return
       endif
     else
-      !if (allocated(kmesh_input%shell_list)) deallocate (kmesh_input%shell_list)
       ! this is the default allocation of the shell_list--used by kmesh_shell_automatic()
       allocate (kmesh_input%shell_list(kmesh_input%max_shells_h), stat=ierr)
       if (ierr /= 0) then
@@ -886,7 +905,8 @@ subroutine w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
     call w90_readwrite_get_keyword(settings, 'num_shells', found, error, comm, i_value=itmp)
     if (allocated(error)) return
     if (found .and. (itmp /= kmesh_input%num_shells)) then
-      call set_error_input(error, 'Error: Found obsolete keyword num_shells. Its value does not agree with shell_list', comm)
+      call set_error_input(error, &
+                           'Error: Found obsolete keyword num_shells. Its value does not agree with shell_list', comm)
       return
     endif
 
@@ -1030,24 +1050,98 @@ end subroutine w90_readwrite_read_lattice
 
   subroutine w90_readwrite_read_atoms(settings, atom_data, real_lattice, bohr, error, comm)
     use w90_error, only: w90_error_type, set_error_input
+    use w90_utility, only: utility_cart_to_frac
+
     implicit none
-    type(atom_data_type), intent(inout) :: atom_data
-    real(kind=dp), intent(in) :: real_lattice(3, 3)
+
+    ! arguments
     real(kind=dp), intent(in) :: bohr
-    type(w90_error_type), allocatable, intent(out) :: error
-    type(w90_comm_type), intent(in) :: comm
+    real(kind=dp), intent(in) :: real_lattice(3, 3)
+    type(atom_data_type), intent(inout) :: atom_data
     type(settings_type), intent(inout) :: settings
+    type(w90_comm_type), intent(in) :: comm
+    type(w90_error_type), allocatable, intent(out) :: error
 
-    integer :: i_temp, i_temp2
-    logical :: found, found2, lunits
+    ! local variables
+    character(len=maxlen), allocatable :: atoms_label_tmp(:)
+    integer :: i_temp, i_temp2, loop, nsymb
+    logical :: found, found2, found3, lunits
+    real(kind=dp), allocatable :: atoms_pos_cart_tmp(:, :)
+    real(kind=dp), allocatable :: atoms_pos_frac_tmp(:, :)
 
-    if (allocated(settings%entries)) return ! don't attempt this read in library mode
+    found = .false.
+    found2 = .false.
+    found3 = .false.
+
+    if (allocated(settings%entries)) then
+      call w90_readwrite_get_vector_length(settings, 'symbols', found, nsymb, error, comm)
+      if (allocated(error)) return
+      call w90_readwrite_get_vector_length(settings, 'atoms_cart', found2, i_temp, error, comm)
+      if (allocated(error)) return
+      call w90_readwrite_get_vector_length(settings, 'atoms_frac', found3, i_temp, error, comm)
+      if (allocated(error)) return
+
+      if (.not. (found .or. found2 .or. found3)) then
+        return ! neither specified, not necessarily an error (only needed if projectors wanted)
+      endif
+
+      ! if supplied, need both entries: labels and positions
+      if (.not. (found .and. (found2 .or. found3))) then
+        call set_error_input(error, 'Error: Must specify both symbols and atoms_frac (or atoms_cart)', comm)
+        return
+      endif
+
+      if (found) atom_data%num_atoms = nsymb ! shape of symbols is n, i_temp returns n
+
+      allocate (atoms_label_tmp(atom_data%num_atoms))
+      allocate (atoms_pos_cart_tmp(3, atom_data%num_atoms))
+
+      ! get symbols list
+      if (found) then
+        call w90_readwrite_get_keyword_vector(settings, 'symbols', found, i_temp, error, comm, &
+                                              c2_value=atoms_label_tmp)
+        if (allocated(error)) return
+      endif
+
+      if (found2) then
+        call w90_readwrite_get_keyword_vector(settings, 'atoms_cart', found, i_temp, error, comm, &
+                                              r2_value=atoms_pos_cart_tmp)
+        if (allocated(error)) return
+      endif
+
+      if (found3) then
+        allocate (atoms_pos_frac_tmp(3, atom_data%num_atoms))
+        call w90_readwrite_get_keyword_vector(settings, 'atoms_frac', found, i_temp, error, comm, &
+                                              r2_value=atoms_pos_frac_tmp)
+        if (allocated(error)) return
+
+        do loop = 1, atom_data%num_atoms
+          call utility_cart_to_frac(atoms_pos_frac_tmp(:, loop), &
+                                    atoms_pos_cart_tmp(:, loop), transpose(real_lattice))
+        enddo
+        deallocate (atoms_pos_frac_tmp)
+      endif
+
+      call w90_readwrite_set_atoms(atom_data, atoms_label_tmp, atoms_pos_cart_tmp, error, comm)
+      if (allocated(error)) return
+
+      deallocate (atoms_label_tmp)
+      deallocate (atoms_pos_cart_tmp)
+
+      return ! no futher action in library mode
+    endif
+
+    i_temp = 0
+    i_temp2 = 0
+    found = .false.
+    found2 = .false.
 
     ! Atoms
     call w90_readwrite_get_block_length(settings, 'atoms_frac', found, i_temp, error, comm)
     if (allocated(error)) return
     call w90_readwrite_get_block_length(settings, 'atoms_cart', found2, i_temp2, error, comm, lunits)
     if (allocated(error)) return
+
     if (found .and. found2) then
       call set_error_input(error, 'Error: Cannot specify both atoms_frac and atoms_cart', comm)
       return
@@ -1549,9 +1643,7 @@ function w90_readwrite_get_convention_type(sc_phase_conv)
     !! associated to a sc_phase_conv integer value.
     integer, intent(in) :: sc_phase_conv
     !! The integer index for which we want to get the string
-    character(len=80)   :: w90_readwrite_get_convention_type
-
-    !character(len=4)   :: orderstr
+    character(len=80) :: w90_readwrite_get_convention_type
 
     if (sc_phase_conv .eq. 1) then
       w90_readwrite_get_convention_type = "Tight-binding convention"
@@ -2545,7 +2637,8 @@ end subroutine w90_readwrite_get_keyword
 
   !================================================!
   subroutine w90_readwrite_get_keyword_vector(settings, keyword, found, length, error, comm, &
-                                              c_value, l_value, i_value, r_value)
+                                              c_value, l_value, i_value, r_value, r2_value, &
+                                              c2_value)
     !================================================!
     !
     !! Finds the values of the required keyword vector
@@ -2572,6 +2665,10 @@ subroutine w90_readwrite_get_keyword_vector(settings, keyword, found, length, er
     !! Keyword data
     real(kind=dp), optional, intent(inout) :: r_value(length)
     !! Keyword data
+    real(kind=dp), allocatable, optional, intent(inout) :: r2_value(:, :)
+    !! Keyword data
+    character(len=*), allocatable, optional, intent(inout) :: c2_value(:)
+    !! Keyword data
     type(settings_type), intent(inout) :: settings
 
     integer           :: kl, in, loop, i
@@ -2594,6 +2691,10 @@ subroutine w90_readwrite_get_keyword_vector(settings, keyword, found, length, er
             i_value = settings%entries(loop)%i1d
           else if (present(r_value)) then
             r_value = settings%entries(loop)%r1d
+          else if (present(r2_value)) then
+            r2_value = settings%entries(loop)%r2d
+          else if (present(c2_value)) then
+            c2_value = settings%entries(loop)%c2d
           else
             call set_error_fatal(error, 'Error: vector sought, but no variable provided to assign to. (readwrite.F90)', comm)
             return
@@ -2667,16 +2768,6 @@ subroutine w90_readwrite_get_vector_length(settings, keyword, found, length, err
 
     found = .false.
 
-    ! get_vector_length only is meaningful for human text in input file
-    ! not suitable for data passed via library interface (data in settings%entries)
-    !if (.not.allocated(settings%in_data)) then
-    !  call set_error_fatal(error, 'w90_readwrite_get_vector_length called with no input file (seeking '//trim(keyword)//')', comm)
-    !  return
-    !elseif (allocated(settings%entries)) then
-    !  call set_error_fatal(error, 'w90_readwrite_get_vector_length called with unspent option arrays', comm)
-    !  return
-    !endif
-
     if (allocated(settings%entries) .and. allocated(settings%in_data)) then
       call set_error_fatal(error, 'Error: (library use) options interface and .win parsing clash.'// &
                            '  See library documentation "setting options." (readwrite.F90)', comm)
@@ -2690,8 +2781,12 @@ subroutine w90_readwrite_get_vector_length(settings, keyword, found, length, err
             length = size(settings%entries(loop)%i1d)
           else if (allocated(settings%entries(loop)%r1d)) then
             length = size(settings%entries(loop)%r1d)
+          else if (allocated(settings%entries(loop)%r2d)) then
+            length = size(settings%entries(loop)%r2d, 1)
+          else if (allocated(settings%entries(loop)%c2d)) then
+            length = size(settings%entries(loop)%c2d, 1)
           else
-            call set_error_input(error, 'lib array not i or r', comm)
+            call set_error_input(error, 'lib array not i or r, r2d or c2d', comm)
           endif
           found = .true.
         end if
@@ -2940,11 +3035,16 @@ subroutine w90_readwrite_get_block_length(settings, keyword, found, rows, error,
     !! Have we found a unit specification
     type(settings_type), intent(inout) :: settings
 
-    integer           :: i, in, ins, ine, loop, line_e, line_s
-    logical           :: found_e, found_s
+    integer :: i, in, ins, ine, loop, line_e, line_s
+    logical :: found_e, found_s
     character(len=maxlen) :: end_st, start_st, dummy
-    character(len=2)  :: atsym
-    real(kind=dp)     :: atpos(3)
+    character(len=2) :: atsym
+    real(kind=dp) :: atpos(3)
+
+    found = .false.
+    rows = 0
+    found_s = .false.
+    found_e = .false.
 
     ! get_block_length only is meaningful for human text in input file
     ! not suitable for data passed via library interface (data in settings%entries)
@@ -2956,13 +3056,8 @@ subroutine w90_readwrite_get_block_length(settings, keyword, found, rows, error,
     !  return
     !endif
 
-    found = .false.
     if (allocated(settings%entries)) return ! don't try to do this in library mode (i.e. when not reading win file)
 
-    rows = 0
-    found_s = .false.
-    found_e = .false.
-
     start_st = 'begin '//trim(keyword)
     end_st = 'end '//trim(keyword)
 
@@ -3230,14 +3325,12 @@ subroutine readwrite_get_atoms(settings, atom_data, lunits, real_lattice, bohr,
   end subroutine readwrite_get_atoms
 
   !================================================!
-  subroutine w90_readwrite_lib_set_atoms(atom_data, atoms_label_tmp, atoms_pos_cart_tmp, &
-                                         real_lattice, error, comm)
+  subroutine w90_readwrite_set_atoms(atom_data, atoms_label_tmp, atoms_pos_cart_tmp, error, comm)
     !================================================!
     !
     !!   Fills the atom data block during a library call
     !
     !================================================!
-
     use w90_utility, only: utility_cart_to_frac, utility_inverse_mat, utility_lowercase
     use w90_error, only: w90_error_type, set_error_alloc
 
@@ -3246,24 +3339,15 @@ subroutine w90_readwrite_lib_set_atoms(atom_data, atoms_label_tmp, atoms_pos_car
     type(atom_data_type), intent(inout) :: atom_data
     type(w90_error_type), allocatable, intent(out) :: error
     type(w90_comm_type), intent(in) :: comm
-    character(len=*), intent(in) :: atoms_label_tmp(atom_data%num_atoms)
+    character(len=*), intent(in) :: atoms_label_tmp(:)
     !! Atom labels
     real(kind=dp), intent(in) :: atoms_pos_cart_tmp(3, atom_data%num_atoms)
-    !! Atom positions
-    real(kind=dp), intent(in) :: real_lattice(3, 3)
+    !! Atom positions, Cartesian, Angstrom
 
-    real(kind=dp) :: inv_lattice(3, 3)
-    real(kind=dp) :: atoms_pos_frac_tmp(3, atom_data%num_atoms)
     integer :: loop2, max_sites, ierr, ic, loop, counter
     character(len=maxlen) :: ctemp(atom_data%num_atoms)
     character(len=maxlen) :: tmp_string
 
-    call utility_inverse_mat(real_lattice, inv_lattice)
-    do loop = 1, atom_data%num_atoms
-      call utility_cart_to_frac(atoms_pos_cart_tmp(:, loop), &
-                                atoms_pos_frac_tmp(:, loop), inv_lattice)
-    enddo
-
     ! Now we sort the data into the proper structures
     atom_data%num_species = 1
     ctemp(1) = atoms_label_tmp(1)
@@ -3292,6 +3376,7 @@ subroutine w90_readwrite_lib_set_atoms(atom_data, atoms_label_tmp, atoms_pos_car
       call set_error_alloc(error, 'Error allocating atoms_symbol in w90_readwrite_lib_set_atoms', comm)
       return
     endif
+
     atom_data%species_num(:) = 0
 
     do loop = 1, atom_data%num_species
@@ -3315,7 +3400,6 @@ subroutine w90_readwrite_lib_set_atoms(atom_data, atoms_label_tmp, atoms_pos_car
       do loop2 = 1, atom_data%num_atoms
         if (trim(atom_data%label(loop)) == trim(atoms_label_tmp(loop2))) then
           counter = counter + 1
-          !atom_data%pos_frac(:, counter, loop) = atoms_pos_frac_tmp(:, loop2)
           atom_data%pos_cart(:, counter, loop) = atoms_pos_cart_tmp(:, loop2)
         end if
       end do
@@ -3337,7 +3421,7 @@ subroutine w90_readwrite_lib_set_atoms(atom_data, atoms_label_tmp, atoms_pos_car
       if ((ic .ge. ichar('a')) .and. (ic .le. ichar('z'))) &
         atom_data%label(loop) (1:1) = char(ic + ichar('Z') - ichar('z'))
     end do
-  end subroutine w90_readwrite_lib_set_atoms
+  end subroutine w90_readwrite_set_atoms
 
   !================================================!
   subroutine w90_readwrite_get_range_vector(settings, keyword, found, length, lcount, error, comm, i_value)
@@ -3377,38 +3461,47 @@ subroutine w90_readwrite_get_range_vector(settings, keyword, found, length, lcou
       return
     endif
 
-    if (allocated(settings%entries)) then ! shortcut for library case
-      if (lcount) then
-        call w90_readwrite_get_vector_length(settings, keyword, found, length, error, comm)
-        return
-      else
-        call w90_readwrite_get_keyword_vector(settings, keyword, found, length, error, comm, &
-                                              i_value=i_value)
-        return
-      endif
-    endif ! end library branch
-
     kl = len_trim(keyword)
 
     found = .false.
-
-    do loop = 1, settings%num_lines
-      in = index(settings%in_data(loop), trim(keyword))
-      if (in == 0 .or. in > 1) cycle
-      if (found) then
-        call set_error_input(error, 'Error: Found keyword '//trim(keyword) &
-                             //' more than once in input file', comm)
-        return
-      endif
-      found = .true.
-      dummy = settings%in_data(loop) (kl + 1:)
-      dummy = adjustl(dummy)
-      if (.not. lcount) settings%in_data(loop) (1:maxlen) = ' '
-      if (dummy(1:1) == '=' .or. dummy(1:1) == ':') then
-        dummy = dummy(2:)
+    if (allocated(settings%entries)) then !  library case
+      do loop = 1, settings%num_entries  ! the first occurance of the variable in settings is used
+        if (settings%entries(loop)%keyword == trim(keyword)) then
+          found = .true.
+          if (allocated(settings%entries(loop)%i1d)) then
+            if (lcount) then
+              call w90_readwrite_get_vector_length(settings, keyword, found, length, error, comm)
+              return
+            else
+              call w90_readwrite_get_keyword_vector(settings, keyword, found, length, error, comm, &
+                                                    i_value=i_value)
+              return
+            endif
+          else
+            dummy = settings%entries(loop)%txtdata
+            dummy = adjustl(dummy)
+          endif
+        endif
+      enddo
+    else ! usual input (.win) file read
+      do loop = 1, settings%num_lines
+        in = index(settings%in_data(loop), trim(keyword))
+        if (in == 0 .or. in > 1) cycle
+        if (found) then
+          call set_error_input(error, 'Error: Found keyword '//trim(keyword) &
+                               //' more than once in input file', comm)
+          return
+        endif
+        found = .true.
+        dummy = settings%in_data(loop) (kl + 1:)
         dummy = adjustl(dummy)
-      end if
-    end do
+        if (.not. lcount) settings%in_data(loop) (1:maxlen) = ' '
+        if (dummy(1:1) == '=' .or. dummy(1:1) == ':') then
+          dummy = dummy(2:)
+          dummy = adjustl(dummy)
+        end if
+      end do
+    endif
 
     if (.not. found) return
 
@@ -3664,6 +3757,10 @@ subroutine w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann
     start_st = 'begin '//trim(keyword)
     end_st = 'end '//trim(keyword)
 
+    ierr = 0
+
+    if (allocated(input_proj)) return ! projectors have already been read, return
+
     if (.not. lcount) then
       allocate (input_proj(num_proj), stat=ierr)
       if (ierr /= 0) then
@@ -3746,13 +3843,12 @@ subroutine w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann
     elseif (allocated(settings%entries)) then ! reading from setopt
       do loop = 1, settings%num_entries
         if (settings%entries(loop)%keyword == 'projections') then
-          counter = counter + 1
           if (settings%entries(loop)%txtdata == 'bohr') lconvert = .true.
           if (settings%entries(loop)%txtdata == 'random') lrandom = .true.
         endif
       enddo
       line_s = 0
-      line_e = settings%num_entries + 1
+      line_e = settings%num_entries
     endif ! reading from input file or entries
 
     counter = 0
@@ -3779,8 +3875,10 @@ subroutine w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann
         else
           dummy = utility_strip(settings%in_data(line))
         endif
+        if (len(trim(dummy)) == 0) cycle
         dummy = adjustl(dummy)
         pos1 = index(dummy, ':')
+
         if (pos1 == 0) then
           call set_error_input(error, &
                                'w90_wannier90_readwrite_read_projection: malformed projection definition: '//trim(dummy), comm)
diff --git a/src/types.F90 b/src/types.F90
index 02c91d65b..d26733080 100644
--- a/src/types.F90
+++ b/src/types.F90
@@ -249,6 +249,7 @@ module w90_types
     ! as different types; otherwise reshape, etc.
     character(len=:), allocatable :: keyword ! token
     character(len=:), allocatable :: txtdata ! text data item
+    character(len=:), allocatable :: c2d(:)
     ! integer data
     integer, allocatable :: i1d(:)
     integer, allocatable :: i2d(:, :)
diff --git a/src/wannier90_readwrite.F90 b/src/wannier90_readwrite.F90
index b4433a875..5876b6500 100644
--- a/src/wannier90_readwrite.F90
+++ b/src/wannier90_readwrite.F90
@@ -122,22 +122,21 @@ subroutine w90_wannier90_readwrite_read_special(settings, atom_data, kmesh_input
     total_bands = 0
     call w90_readwrite_read_total_bands(settings, total_bands, error, comm)
     if (allocated(error)) return
-
     if (total_bands > 0) then
       num_bands = total_bands - num_exclude_bands
     else
-      ! fixme, flag up that it is a mistake to set total bands and num_bands
-      call w90_readwrite_read_num_bands(settings, .false., num_bands, num_wann, stdout, error, comm)
+      call w90_readwrite_read_num_bands(settings, .false., num_bands, num_wann, error, comm)
       if (allocated(error)) return
     endif
 
     disentanglement = (num_bands > num_wann)
 
-    num_proj = num_wann !default, no projections specified
-
     call w90_readwrite_read_mp_grid(settings, .false., mp_grid, num_kpts, error, comm)
     if (allocated(error)) return
 
+    call w90_readwrite_read_distk(settings, distk, num_kpts, error, comm)
+    if (allocated(error)) return
+
     call w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
     if (allocated(error)) return
 
@@ -151,21 +150,22 @@ subroutine w90_wannier90_readwrite_read_special(settings, atom_data, kmesh_input
     call w90_readwrite_read_system(settings, w90_system, error, comm)
     if (allocated(error)) return
 
+    num_proj = num_wann !default, no projections specified
+
     call utility_inverse_mat(real_lattice, inv_lattice)
     call w90_wannier90_readwrite_read_projections(settings, proj, proj_input, use_bloch_phases, &
                                                   lhasproj, wann_control%guiding_centres%enable, &
                                                   select_proj, num_proj, atom_data, inv_lattice, &
                                                   num_wann, gamma_only, w90_system%spinors, bohr, &
                                                   stdout, error, comm)
+    if (allocated(error)) return
+
     if (allocated(proj)) then
       allocate (wann_control%guiding_centres%centres(3, num_proj))
       do ip = 1, num_proj
         wann_control%guiding_centres%centres(:, ip) = proj(ip)%site(:)
       enddo
     endif
-    if (allocated(error)) return
-
-    call w90_readwrite_read_distk(settings, distk, num_kpts, error, comm)
 
   end subroutine w90_wannier90_readwrite_read_special
 
@@ -271,26 +271,10 @@ subroutine w90_wannier90_readwrite_read(settings, band_plot, dis_control, dis_sp
                                     energy_unit, bohr, error, comm)
       if (allocated(error)) return
 
-      !call w90_readwrite_read_num_wann(settings, num_wann, error, comm)
-      !if (allocated(error)) return
-
-      !call w90_readwrite_read_exclude_bands(settings, exclude_bands, num_exclude_bands, error, comm)
-      !if (allocated(error)) return
-
-      !call w90_readwrite_read_num_bands(settings, .false., num_bands, num_wann, stdout, error, comm)
-      !if (allocated(error)) return
-      !disentanglement = (num_bands > num_wann)
-
-      !call w90_readwrite_read_lattice(settings, real_lattice, bohr, error, comm)
-      !if (allocated(error)) return
-
       call w90_wannier90_readwrite_read_wannierise(settings, wann_control, num_wann, &
                                                    w90_extra_io%ccentres_frac, stdout, error, comm)
       if (allocated(error)) return
 
-      !call w90_readwrite_read_mp_grid(settings, .false., mp_grid, num_kpts, error, comm)
-      !if (allocated(error)) return
-
       call w90_readwrite_read_gamma_only(settings, gamma_only, num_kpts, error, comm)
       if (allocated(error)) return
 
@@ -301,9 +285,6 @@ subroutine w90_wannier90_readwrite_read(settings, band_plot, dis_control, dis_sp
       call w90_wannier90_readwrite_read_restart(settings, w90_calculation, seedname, error, comm)
       if (allocated(error)) return
 
-      !call w90_readwrite_read_system(settings, w90_system, error, comm)
-      !if (allocated(error)) return
-
       call w90_readwrite_read_kpath(settings, kpoint_path, has_kpath, w90_calculation%bands_plot, &
                                     error, comm)
       if (allocated(error)) return
@@ -349,7 +330,7 @@ subroutine w90_wannier90_readwrite_read(settings, band_plot, dis_control, dis_sp
     if (allocated(error)) return
 
     if (.not. (w90_calculation%transport .and. tran%read_ht)) then
-      call w90_readwrite_read_dis_manifold(settings, .true., dis_manifold, error, comm)
+      call w90_readwrite_read_dis_manifold(settings, dis_manifold, error, comm)
       if (allocated(error)) return
 
       call w90_wannier90_readwrite_read_disentangle(settings, dis_control, dis_spheres, num_bands, &
@@ -363,35 +344,6 @@ subroutine w90_wannier90_readwrite_read(settings, band_plot, dis_control, dis_sp
                                                     error, comm)
       if (allocated(error)) return
 
-      !call w90_readwrite_read_kmesh_data(settings, kmesh_input, error, comm)
-      !if (allocated(error)) return
-
-      !call utility_recip_lattice(real_lattice, recip_lattice, volume, error, comm)
-      !if (allocated(error)) return
-
-      !call utility_inverse_mat(real_lattice, inv_lattice)
-
-      !call w90_readwrite_read_kpoints(settings, .false., kpt_latt, num_kpts, bohr, error, comm)
-      !if (allocated(error)) return
-
-      !call w90_wannier90_readwrite_read_explicit_kpts(settings, w90_calculation, kmesh_info, &
-      !                                                num_kpts, bohr, error, comm)
-      !if (allocated(error)) return
-
-      !call w90_wannier90_readwrite_read_global_kmesh(global_kmesh_set, kmesh_spacing, kmesh, recip_lattice, &
-      !                             stdout, seedname)
-      !call w90_readwrite_read_atoms(settings, atom_data, real_lattice, bohr, error, comm)
-      !if (allocated(error)) return
-
-      !call w90_wannier90_readwrite_read_projections(settings, proj, proj_input, use_bloch_phases, lhasproj, &
-      !                                              wann_control%guiding_centres%enable, &
-      !                                              select_proj, num_proj, atom_data, inv_lattice, &
-      !                                              num_bands, num_wann, gamma_only, &
-      !                                              w90_system%spinors, bohr, stdout, error, comm)
-      !if (allocated(error)) return
-
-      ! projections needs to be allocated before reading constrained centres
-
       if (wann_control%constrain%constrain) then
         call w90_wannier90_readwrite_read_constrained_centres(settings, w90_extra_io%ccentres_frac, &
                                                               wann_control, real_lattice, &
@@ -1529,12 +1481,15 @@ subroutine w90_wannier90_readwrite_read_projections(settings, proj, proj_input,
                                                       gamma_only, spinors, bohr, stdout, error, &
                                                       comm)
     !================================================!
+    ! Obtain projector definitions
     use w90_error, only: w90_error_type
 
     implicit none
 
+    ! arguments
     integer, intent(in) :: num_wann
     integer, intent(inout) :: num_proj
+    integer, intent(in) :: stdout
     logical, intent(in) :: gamma_only
     logical, intent(in) :: spinors
     logical, intent(in) :: use_bloch_phases, guiding_centres
@@ -1548,7 +1503,7 @@ subroutine w90_wannier90_readwrite_read_projections(settings, proj, proj_input,
     type(w90_comm_type), intent(in) :: comm
     type(w90_error_type), allocatable, intent(out) :: error
 
-    integer, intent(in) :: stdout
+    ! local variables
     integer :: i, j, i_temp, loop, ierr
     logical :: found
     ! projections selection
@@ -1556,6 +1511,8 @@ subroutine w90_wannier90_readwrite_read_projections(settings, proj, proj_input,
     integer, allocatable :: select_projections(:)
     integer :: imap
 
+    lhasproj = .false.
+
     ! Projections
     call w90_readwrite_get_keyword(settings, 'auto_projections', found, error, comm, &
                                    l_value=select_proj%auto_projections)
@@ -1590,15 +1547,13 @@ subroutine w90_wannier90_readwrite_read_projections(settings, proj, proj_input,
         return
       endif
 
-      call w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann, proj_input, &
-                                         recip_lattice, .true., spinors, bohr, stdout, error, comm)
-      ! count number of projections first
-      !write(*,*)"Projections found: ",num_proj
-      call w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann, proj_input, &
-                                         recip_lattice, .false., spinors, bohr, stdout, error, comm)
-      !do ip = 1, num_proj
-      !  write(*,*)"site,l,m: ",proj_input(ip)%site,proj_input(ip)%l,proj_input(ip)%m
-      !enddo
+      if (atom_data%num_species > 0) then
+        ! only read projections if atom positions are known
+        call w90_readwrite_get_projections(settings, num_proj, atom_data, num_wann, proj_input, &
+                                           recip_lattice, .false., spinors, bohr, stdout, error, comm)
+        if (allocated(error)) return
+        lhasproj = .true.
+      endif
     endif
 
     num_select_projections = num_proj !num proj is the size of proj_input
@@ -1683,7 +1638,7 @@ subroutine w90_wannier90_readwrite_read_projections(settings, proj, proj_input,
         endif
         proj(imap) = proj_input(loop)
       enddo
-    endif
+    endif !lhasproj
 
   end subroutine w90_wannier90_readwrite_read_projections