Proper protocol for Normal Mode Analysis

GROMACS version: 2020.3 (double precision)
GROMACS modification: No
What is the suggested procedure for normal mode analysis in recent version of gromacs (since gromacs 2020)? I understand there are three steps: Energy minimization, Hessian calculation and its diagonalization. The only step I have no trouble with is the Hessian calculation.

In the energy minimization using following parameters I am not able to get the force as low as I would like.

define = -DFLEXIBLE
integrator  = l-bfgs        
emtol       = 0.000001      
emstep      = 0.01         
nsteps      = 100000        

constraints = None
nstlist         = 5         
cutoff-scheme   = verlet   
coulombtype = Reaction-Field
coulomb-modifier = None
epsilon-rf      = 0
rcoulomb        = 1.2
rlist           = 1.4
rvdw            = 1.2      
vdwtype         = Cut-off
vdw-modifier    = Potential-switch
rvdw-switch     = 1.0
pbc             =  xyz      

This way I am able to reach maximum force in range 1e-4 at best. I remember that in previous versions of gromacs with no periodic boundaries and switched potentials for short-range forces it works rather well (different system). However this setup is no longer possible in recent gromacs versions.

Diagonalization of Hessian (gmx nmeig). The system consists of around 11k atoms. When I demand first 10k eigenmodes the diagonalization ends up with the error.

Diagonalizing to find eigenvectors 1 through 10000...
Calculation Ritz values and Lanczos vectors, max 10000000 iterations...
Iteration 20001: 4325 out of 10000 Ritz values converged.

Program:     gmx nmeig, version 2020.3
Source file: src/gromacs/linearalgebra/eigensolver.cpp (line 346)

Fatal error:
Unspecified error from Arnoldi diagonalization:3

However demanding only first 4k it finished properly after approximately 200k iterations providing all eigenvectors and eigenvalues.

Thank you in advance for any reply and suggestions!