Added tests, made DOS splitting possible from QE out files as well

Author: RandomDefaultUser

mala/targets/dos.py

@@ -206,7 +206,15 @@ def feature_size(self):
         if isinstance(self.parameters.ldos_gridsize, int):
             return self.parameters.ldos_gridsize
         elif isinstance(self.parameters.ldos_gridsize, list):
-            return np.sum(self.parameters.ldos_gridsize)
+            # For splits, we sum up the individual grid sizes, BUT we
+            # have to subtract one for each split, as the last energy
+            # of each section gets discarded. So for three sections,
+            # we have to subtract 2.
+            return (
+                np.sum(self.parameters.ldos_gridsize)
+                - len(self.parameters.ldos_gridsize)
+                + 1
+            )
 
     @property
     def data_name(self):
@@ -514,17 +522,40 @@ def read_from_qe_out(self, path=None, smearing_factor=2):
                 "Rerun calculation with verbosity set to 'high'."
             )
 
+        if isinstance(self.parameters.ldos_gridspacing_ev, list):
+            grid_spacings = self.parameters.ldos_gridspacing_ev
+            grid_sizes = self.parameters.ldos_gridsize
+        else:
+            grid_spacings = [self.parameters.ldos_gridspacing_ev]
+            grid_sizes = [self.parameters.ldos_gridsize]
+
         # Get the gaussians for all energy values and calculate the DOS per
         # band.
-        dos_per_band = gaussians(
-            self.energy_grid,
-            atoms_object.get_calculator().band_structure().energies[0, :, :],
-            smearing_factor * self.parameters.ldos_gridspacing_ev,
-        )
-        dos_per_band = kweights[:, np.newaxis, np.newaxis] * dos_per_band
+        dos_data = None
+        previous_beginning = 0
+        for spacing_idx, grid_spacing in enumerate(grid_spacings):
+            size_for_spacing = grid_sizes[spacing_idx] + previous_beginning
+            if spacing_idx != len(grid_spacings) - 1:
+                size_for_spacing -= 1
+
+            dos_per_band = gaussians(
+                self.energy_grid[previous_beginning:size_for_spacing],
+                atoms_object.get_calculator()
+                .band_structure()
+                .energies[0, :, :],
+                smearing_factor * grid_spacing,
+            )
+            dos_per_band = kweights[:, np.newaxis, np.newaxis] * dos_per_band
 
-        # QE gives the band energies in eV, so no conversion necessary here.
-        dos_data = np.sum(dos_per_band, axis=(0, 1))
+            # QE gives the band energies in eV, so no conversion necessary
+            # here.
+            if spacing_idx == 0:
+                dos_data = np.sum(dos_per_band, axis=(0, 1))
+            else:
+                dos_data = np.concatenate(
+                    (dos_data, np.sum(dos_per_band, axis=(0, 1)))
+                )
+            previous_beginning = size_for_spacing
         self.density_of_states = dos_data
         return dos_data
 

mala/targets/ldos.py

@@ -198,7 +198,15 @@ def feature_size(self):
         if isinstance(self.parameters.ldos_gridsize, int):
             return self.parameters.ldos_gridsize
         elif isinstance(self.parameters.ldos_gridsize, list):
-            return np.sum(self.parameters.ldos_gridsize)
+            # For splits, we sum up the individual grid sizes, BUT we
+            # have to subtract one for each split, as the last energy
+            # of each section gets discarded. So for three sections,
+            # we have to subtract 2.
+            return (
+                np.sum(self.parameters.ldos_gridsize)
+                - len(self.parameters.ldos_gridsize)
+                + 1
+            )
 
     @property
     def data_name(self):

test/splitting_test.py

@@ -0,0 +1,84 @@
+import os
+
+import numpy as np
+
+import mala
+from mala.datahandling.data_repo import data_path
+
+
+class TestSplitting:
+    """
+    Class that tests that splitted (L)DOS data can be processed
+
+    Splitted means along the energy domain.
+    """
+
+    def test_ldos_splitting(self):
+        parameters3 = mala.Parameters()
+        parameters3.targets.ldos_gridsize = [
+            3,
+            3,
+            7,
+        ]
+        parameters3.targets.ldos_gridoffset_ev = [-5, 0, 5]
+        parameters3.targets.ldos_gridspacing_ev = [2.5, 2.5, 2.5]
+        cubes_path = os.path.join(data_path, "cubes")
+        ldos_calculator3 = mala.LDOS(parameters3)
+        ldos_calculator3.read_from_cube(
+            [
+                os.path.join(cubes_path, "tmp.pp0*Be_part1_ldos.cube"),
+                os.path.join(cubes_path, "tmp.pp0*Be_part2_ldos.cube"),
+                os.path.join(cubes_path, "tmp.pp0*Be_part3_ldos.cube"),
+            ]
+        )
+
+        parameters = mala.Parameters()
+        parameters.targets.ldos_gridsize = 11
+        parameters.targets.ldos_gridspacing_ev = 2.5
+        parameters.targets.ldos_gridoffset_ev = -5
+        ldos_calculator = mala.LDOS(parameters)
+        ldos_calculator.read_from_cube(
+            os.path.join(cubes_path, "tmp.pp*Be_ldos.cube")
+        )
+
+        print(
+            np.mean(
+                np.abs(
+                    ldos_calculator3.local_density_of_states
+                    - ldos_calculator.local_density_of_states
+                )
+            )
+        )
+
+    def test_dos_splitting(self):
+        parameters3 = mala.Parameters()
+        parameters3.targets.ldos_gridsize = [
+            3,
+            3,
+            7,
+        ]
+        parameters3.targets.ldos_gridoffset_ev = [-5, 0, 5]
+        parameters3.targets.ldos_gridspacing_ev = [2.5, 2.5, 2.5]
+        cubes_path = os.path.join(data_path, "cubes")
+        dos_calculator3 = mala.DOS(parameters3)
+        dos_calculator3.read_from_qe_out(
+            os.path.join(data_path, "Be_snapshot1.out")
+        )
+
+        parameters = mala.Parameters()
+        parameters.targets.ldos_gridsize = 11
+        parameters.targets.ldos_gridspacing_ev = 2.5
+        parameters.targets.ldos_gridoffset_ev = -5
+        dos_calculator = mala.DOS(parameters)
+        dos_calculator.read_from_qe_out(
+            os.path.join(data_path, "Be_snapshot1.out")
+        )
+
+        print(
+            np.mean(
+                np.abs(
+                    dos_calculator3.density_of_states
+                    - dos_calculator.density_of_states
+                )
+            )
+        )