sphinx eigenvalue parser problem fix (again)

[ahmedabdelkawy]

In addition to the explanation in #1159 (PR was closed as the branch was damaged), the problem in

def _parse_band(self, term):
       arr = np.loadtxt(re.findall(term, self.log_main, re.MULTILINE))
       shape = (-1, len(self.k_points), arr.shape[-1])
       if self.spin_enabled:
           shape = (-1, 2, len(self.k_points), shape[-1])
       return arr.reshape(shape)

apparently became more severe while doing minimization. This function tried to parse all eigenvalues at the end of each ionic steps (something still not implemented in the vasp parser as far as I know). This then can not be fixed with the fix I introduced in #1159 (arr = np.vstack((arr[::2], arr[1::2]))).

The reason again is that the eigen values in the Sphinx log has the following order:

1. Ionic step = 1
    1.1 KPOINTS = 1
        1.1.1 SPIN channel = 0
        1.1.2 SPIN channel = 1
    1.2 KPOINTS = 2
        1.2.1 SPIN channel = 0
        1.2.2 SPIN channel = 1
    1.3 KPOINTS = 3
        1.1.1 SPIN channel = 0
        1.1.2 SPIN channel = 1
    1.4 KPOINTS = 4
        1.2.1 SPIN channel = 0
        1.2.2 SPIN channel = 1
2. Ionic step = 2
    2.1 KPOINTS = 1
        2.1 SPIN channel = 0
        2.1.2 SPIN channel = 1
    2.2 KPOINTS = 2
        2.2.1 SPIN channel = 0
        2.2.2 SPIN channel = 1
    2.3 KPOINTS = 3
        2.1.1 SPIN channel = 0
        2.1.2 SPIN channel = 1
    2.4 KPOINTS = 4
        2.2.1 SPIN channel = 0
        2.2.2 SPIN channel = 1

with the shape = (n_ionic_steps, n_kpoints, n_spin_channnels, n_bands)

while the parser expects the following order:

1. Ionic step = 1
    1.1 SPIN channel = 0
        1.1.1 KPOINTS = 1
        1.1.2 KPOINTS = 2
        1.1.3 KPOINTS = 3
        1.1.4 KPOINTS = 4
    1.2 SPIN channel = 1
        1.2.1 KPOINTS = 1
        1.2.2 KPOINTS = 2
        1.2.3 KPOINTS = 3
        1.2.4 KPOINTS = 4
2. Ionic step = 2
    2.1 SPIN channel = 0
        2.1.1 KPOINTS = 1
        2.1.2 KPOINTS = 2
        2.1.3 KPOINTS = 3
        2.1.4 KPOINTS = 4
    2.2 SPIN channel = 1
        2.2.1 KPOINTS = 1
        2.2.2 KPOINTS = 2
        2.2.3 KPOINTS = 3
        2.2.4 KPOINTS = 4

with the shape = (n_ionic_steps, n_spin_channnels, n_kpoints, n_bands)

[ahmedabdelkawy]

To reproduce (and notice) the problem, the system must contain two different spin channels and minimization calculation should be run. Here is an example where the reference DOS is VASP:

structure_fe = pr_ferro.create.structure.bulk('Fe', 'bcc', orthorhombic=True)
structure_fe.positions[0,0] += 0.1
structure_fe.set_initial_magnetic_moments([3,3])

sphinx_job = pr_ferro.create.job.Sphinx(job_name='Fe_min_spx', delete_existing_job=True)
sphinx_job.structure = structure_fe.copy()
sphinx_job.set_encut(450)
sphinx_job.set_kpoints([13,13,13])
sphinx_job.set_occupancy_smearing(smearing='FermiDirac', width= 0.05, order=0)
sphinx_job.set_empty_states(7)
sphinx_job.server.cores = 5
sphinx_job.calc_minimize()
sphinx_job.executable = 'latest'
sphinx_job.run(run_mode='queue')

vasp_job = pr_ferro.create.job.Vasp(job_name='Fe_min_Vasp', delete_existing_job=True)
vasp_job.structure = structure_fe.copy()
vasp_job.set_encut(450)
vasp_job.set_kpoints([13,13,13])
vasp_job.input.incar["ISMEAR"] = -1      #fermi smearing 
vasp_job.input.incar["SIGMA"] = 0.05 
vasp_job.input.incar["NBANDS"] = 20   
vasp_job.input.incar['LORBIT']  = 11
vasp_job.input.incar["ALGO"]= 'Normal'
vasp_job.executable = '5.4.4'
vasp_job.server.cores = 5
vasp_job.calc_minimize()
vasp_job.run(delete_existing_job=True, run_mode='queue')

This is also a nice example for whoever in the future to compare VASP and SPHINX. Modifications in NBANDS, ISMEAR, SIGMA in VASP and set_empty_states and set_occupancy_smearing is for consistency between both codes.

If you now plot the DOS from both calculations:

plt.plot(fe_vasp_job.get_density_of_states()['grid'], fe_vasp_job.get_density_of_states()['dos'][0], color='blue', label='vasp')
plt.plot(fe_vasp_job.get_density_of_states()['grid'], -fe_vasp_job.get_density_of_states()['dos'][1], c='blue')
plt.plot(fe_sphinx_job.get_density_of_states()['grid'], fe_sphinx_job.get_density_of_states()['dos'][0], c='r', label='sphinx')
plt.plot(fe_sphinx_job.get_density_of_states()['grid'], -fe_sphinx_job.get_density_of_states()['dos'][1], c='r')
plt.axvline(0, color="black", linestyle="dashed")
plt.axhline(0, color="black", linestyle="dashed")
plt.title('Fe DOS')
plt.legend()

you will end up with this:

[ahmedabdelkawy]

After #1209 I get this: