anchor.diff:          case default
   
               ! Find the number of determinants with the required spin.
  -            if (abs(mod(Ms,2)) /= mod(sys%nel,2)) call stop_all('set_spin_polarisation','Required Ms not possible.')
  +            if (abs(mod(Ms,2)) /= mod(sys%nel,2) .or. abs(Ms) > sys%nel) then
  +                write (err, err_fmt) Ms
  +                call stop_all(proc_name, err)
  +            end if
   
               sys%nbeta = (sys%nel - Ms)/2
               sys%nalpha = (sys%nel + Ms)/2
  
positive.code:             ! Spin polarization is different (see comments in system) as the
              ! Heisenberg model is a collection of spins rather than electrons.
              ! See comments in init_system and at module-level.
              sys%nel = (sys%lattice%nsites + Ms)/2
              sys%nvirt = (sys%lattice%nsites - Ms)/2
              ! The Heisenberg model doesn't use these values, but they need to be
  
negative.code:         integer, intent(in) :: nbasis, Ms
          type(sys_t), intent(inout) :: sys
  
          character(len=*), parameter :: proc_name = 'set_spin_polarisation'
          character(len=*), parameter :: err_fmt = '("Required Ms not possible: ",i11, ".")'
          character(len=40) :: err
  
          select case(sys%system)
  
          case(heisenberg)
  