run_tests.py

"""
This is a module for unit testing the CP2K parser. The unit tests are run with
a custom backend that outputs the results directly into native python object for
easier and faster analysis.

Each property that has an enumerable list of different possible options is
assigned a new test class, that should ideally test through all the options.

The properties that can have any value imaginable will be tested only for one
specific case inside a test class that is designed for a certain type of run
(MD, optimization, QM/MM, etc.)
"""
import os
import unittest
import logging
import numpy as np
from cp2kparser import CP2KParser
from nomadcore.unit_conversion.unit_conversion import convert_unit

# Setup the logger so that it doesn't spam messages during tests
logger = logging.getLogger("nomad")
logger.setLevel(logging.CRITICAL)


#===============================================================================
def get_results(folder, metainfo_to_keep=None):
    """Get the given result from the calculation in the given folder by using
    the Analyzer in the nomadtoolkit package. Tries to optimize the parsing by
    giving the metainfo_to_keep argument.

    Args:
        folder: The folder relative to the directory of this script where the
            parsed calculation resides.
        metaname: The quantity to extract.
    """
    dirname = os.path.dirname(__file__)
    filename = os.path.join(dirname, folder, "unittest.out")
    parser = CP2KParser(filename, metainfo_to_keep)
    results = parser.parse()
    return results


#===============================================================================
def get_result(folder, metaname):
    results = get_results(folder, metaname)
    result = results[metaname]
    return result


#===============================================================================
class TestErrors(unittest.TestCase):
    """Test misc. error stuations which may occur during the parsing.
    """
    def test_no_file(self):
        self.assertRaises(IOError, get_result, "errors/no_file", "XC_functional")

    def test_invalid_file(self):
        self.assertRaises(RuntimeError, get_result, "errors/invalid_file", "XC_functional")

    def test_invalid_run_type(self):
        self.assertRaises(KeyError, get_result, "errors/invalid_run_type", "XC_functional")

    def test_unknown_version(self):
        get_result("errors/unknown_version", "XC_functional")


#===============================================================================
class TestXCFunctional(unittest.TestCase):
    """Tests that the XC functionals can be properly parsed.
    """

    def test_pade(self):
        xc = get_result("XC_functional/pade", "XC_functional")
        self.assertEqual(xc, "1*LDA_XC_TETER93")

    def test_lda(self):
        xc = get_result("XC_functional/lda", "XC_functional")
        self.assertEqual(xc, "1*LDA_XC_TETER93")

    def test_blyp(self):
        xc = get_result("XC_functional/blyp", "XC_functional")
        self.assertEqual(xc, "1*GGA_C_LYP+1*GGA_X_B88")

    def test_b3lyp(self):
        xc = get_result("XC_functional/b3lyp", "XC_functional")
        self.assertEqual(xc, "1*HYB_GGA_XC_B3LYP")

    def test_olyp(self):
        xc = get_result("XC_functional/olyp", "XC_functional")
        self.assertEqual(xc, "1*GGA_C_LYP+1*GGA_X_OPTX")

    def test_hcth120(self):
        xc = get_result("XC_functional/hcth120", "XC_functional")
        self.assertEqual(xc, "1*GGA_XC_HCTH_120")

    def test_pbe0(self):
        xc = get_result("XC_functional/pbe0", "XC_functional")
        self.assertEqual(xc, "1*HYB_GGA_XC_PBEH")

    def test_pbe(self):
        xc = get_result("XC_functional/pbe", "XC_functional")
        self.assertEqual(xc, "1*GGA_C_PBE+1*GGA_X_PBE")


#===============================================================================
class TestSCFConvergence(unittest.TestCase):
    """Tests whether the convergence status and number of SCF step can be
    parsed correctly.
    """

    def test_converged(self):
        result = get_result("convergence/converged", "single_configuration_calculation_converged")
        self.assertTrue(result)

    def test_non_converged(self):
        result = get_result("convergence/non_converged", "single_configuration_calculation_converged")
        self.assertFalse(result)


#===============================================================================
class TestForceFiles(unittest.TestCase):
    """Tests that different force files that can be output, can actually be
    found and parsed.
    """

    def test_single_point(self):

        # results = get_results("force_file/single_point")
        # print results._results
        result = get_result("force_file/single_point", "atom_forces")
        expected_result = convert_unit(
            np.array([
                [0.00000000, 0.00000000, 0.00000000],
                [0.00000000, 0.00000001, 0.00000001],
                [0.00000001, 0.00000001, 0.00000000],
                [0.00000001, 0.00000000, 0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
            ]),
            "forceAu"
        )
        self.assertTrue(np.array_equal(result, expected_result))


#===============================================================================
class TestSelfInteractionCorrectionMethod(unittest.TestCase):
    """Tests that the self-interaction correction can be properly parsed.
    """

    def test_no(self):
        sic = get_result("sic/no", "self_interaction_correction_method")
        self.assertEqual(sic, "")

    def test_ad(self):
        sic = get_result("sic/ad", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_AD")

    def test_explicit_orbitals(self):
        sic = get_result("sic/explicit_orbitals", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_EXPLICIT_ORBITALS")

    def test_mauri_spz(self):
        sic = get_result("sic/mauri_spz", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_MAURI_SPZ")

    def test_mauri_us(self):
        sic = get_result("sic/mauri_us", "self_interaction_correction_method")
        self.assertEqual(sic, "SIC_MAURI_US")


#===============================================================================
class TestStressTensorMethods(unittest.TestCase):
    """Tests that the stress tensor can be properly parsed for different
    calculation methods.
    """
    def test_none(self):
        get_results("stress_tensor/none", "section_stress_tensor")

    def test_analytical(self):
        results = get_results("stress_tensor/analytical", ["stress_tensor_method", "stress_tensor"])
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Analytical")

    def test_numerical(self):
        results = get_results("stress_tensor/numerical", ["stress_tensor_method", "stress_tensor"])
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Numerical")

    def test_diagonal_analytical(self):
        results = get_results("stress_tensor/diagonal_analytical", ["stress_tensor_method", "stress_tensor"])
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Diagonal analytical")

    def test_diagonal_numerical(self):
        results = get_results("stress_tensor/diagonal_numerical", ["stress_tensor_method", "stress_tensor"])
        method = results["stress_tensor_method"]
        results["stress_tensor"]
        self.assertEqual(method, "Diagonal numerical")


#===============================================================================
class TestConfigurationPeriodicDimensions(unittest.TestCase):
    """Tests that the self-interaction correction can be properly parsed.
    """

    def test_default(self):
        result = get_result("configuration_periodic_dimensions/default", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((True, True, True))))

    def test_none(self):
        result = get_result("configuration_periodic_dimensions/none", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((False, False, False))))

    def test_x(self):
        result = get_result("configuration_periodic_dimensions/x", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((True, False, False))))

    def test_y(self):
        result = get_result("configuration_periodic_dimensions/y", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((False, True, False))))

    def test_z(self):
        result = get_result("configuration_periodic_dimensions/z", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((False, False, True))))

    def test_xy(self):
        result = get_result("configuration_periodic_dimensions/xy", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((True, True, False))))

    def test_xyz(self):
        result = get_result("configuration_periodic_dimensions/xyz", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((True, True, True))))

    def test_xz(self):
        result = get_result("configuration_periodic_dimensions/xz", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((True, False, True))))

    def test_yz(self):
        result = get_result("configuration_periodic_dimensions/yz", "configuration_periodic_dimensions")
        self.assertTrue(np.array_equal(result, np.array((False, True, True))))


#===============================================================================
class TestEnergyForce(unittest.TestCase):
    """Tests for a CP2K calculation with RUN_TYPE ENERGY_FORCE.
    """

    @classmethod
    def setUpClass(cls):
        cls.results = get_results("energy_force", "section_run")

    def test_energy_total_scf_iteration(self):
        energy_total = self.results["energy_total_scf_iteration"]
        expected_result = convert_unit(np.array(-32.2320848878), "hartree")
        self.assertTrue(np.array_equal(energy_total[0], expected_result))

    def test_energy_change_scf_iteration(self):
        energy_change = self.results["energy_change_scf_iteration"]
        expected_result = convert_unit(np.array(-3.22E+01), "hartree")
        self.assertTrue(np.array_equal(energy_change[0], expected_result))

    def test_energy_XC_scf_iteration(self):
        result = self.results["energy_XC_scf_iteration"]
        expected_result = convert_unit(np.array(-9.4555961214), "hartree")
        self.assertTrue(np.array_equal(result[0], expected_result))

    def test_energy_total(self):
        result = self.results["energy_total"]
        expected_result = convert_unit(np.array(-31.297885372811063), "hartree")
        self.assertTrue(np.array_equal(result, expected_result))

    def test_electronic_kinetic_energy(self):
        result = self.results["electronic_kinetic_energy"]
        expected_result = convert_unit(np.array(13.31525592466418), "hartree")
        self.assertTrue(np.array_equal(result, expected_result))

    def test_atom_forces(self):
        atomic_forces = self.results["atom_forces"]
        expected_result = convert_unit(
            np.array([
                [0.00000000, 0.00000000, 0.00000000],
                [0.00000000, 0.00000001, 0.00000001],
                [0.00000001, 0.00000001, 0.00000000],
                [0.00000001, 0.00000000, 0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
                [-0.00000001, -0.00000001, -0.00000001],
            ]),
            "forceAu"
        )
        self.assertTrue(np.array_equal(atomic_forces, expected_result))

    def test_atom_label(self):
        atom_labels = self.results["atom_labels"]
        expected_labels = np.array(8*["Si"])
        self.assertTrue(np.array_equal(atom_labels, expected_labels))

    def test_simulation_cell(self):
        cell = self.results["simulation_cell"]
        n_vectors = cell.shape[0]
        n_dim = cell.shape[1]
        self.assertEqual(n_vectors, 3)
        self.assertEqual(n_dim, 3)
        expected_cell = convert_unit(np.array([[5.431, 0, 0], [0, 5.431, 0], [0, 0, 5.431]]), "angstrom")
        self.assertTrue(np.array_equal(cell, expected_cell))

    def test_number_of_atoms(self):
        n_atoms = self.results["number_of_atoms"]
        self.assertEqual(n_atoms, 8)

    def test_atom_position(self):
        atom_position = self.results["atom_positions"]
        expected_position = convert_unit(np.array([4.073023, 4.073023, 1.357674]), "angstrom")
        self.assertTrue(np.array_equal(atom_position[-1, :], expected_position))

    def test_cp2k_filenames(self):
        input_filename = self.results["x_cp2k_input_filename"]
        expected_input = "si_bulk8.inp"
        self.assertTrue(input_filename, expected_input)

        bs_filename = self.results["x_cp2k_basis_set_filename"]
        expected_bs = "../BASIS_SET"
        self.assertEqual(bs_filename, expected_bs)

        geminal_filename = self.results["x_cp2k_geminal_filename"]
        expected_geminal = "BASIS_GEMINAL"
        self.assertEqual(geminal_filename, expected_geminal)

        potential_filename = self.results["x_cp2k_potential_filename"]
        expected_potential = "../GTH_POTENTIALS"
        self.assertEqual(potential_filename, expected_potential)

        mm_potential_filename = self.results["x_cp2k_mm_potential_filename"]
        expected_mm_potential = "MM_POTENTIAL"
        self.assertEqual(mm_potential_filename, expected_mm_potential)

        coordinate_filename = self.results["x_cp2k_coordinate_filename"]
        expected_coordinate = "__STD_INPUT__"
        self.assertEqual(coordinate_filename, expected_coordinate)

    def test_target_multiplicity(self):
        multiplicity = self.results["spin_target_multiplicity"]
        self.assertEqual(multiplicity, 1)

    def test_total_charge(self):
        charge = self.results["total_charge"]
        self.assertEqual(charge, 0)

    def test_single_configuration_calculation_converged(self):
        result = self.results["single_configuration_calculation_converged"]
        self.assertTrue(result)

    def test_scf_dft_number_of_iterations(self):
        result = self.results["number_of_scf_iterations"]
        self.assertEqual(result, 10)

    def test_single_configuration_to_calculation_method_ref(self):
        result = self.results["single_configuration_to_calculation_method_ref"]
        self.assertEqual(result, 0)

    def test_single_configuration_calculation_to_system_description_ref(self):
        result = self.results["single_configuration_calculation_to_system_ref"]
        self.assertEqual(result, 0)

    def test_stress_tensor(self):
        result = self.results["stress_tensor"]
        expected_result = convert_unit(
            np.array([
                [7.77641684, -0.00000106, -0.00000106],
                [-0.00000106, 7.77641703, -0.00000106],
                [-0.00000106, -0.00000106, 7.77641703],
            ]),
            "GPa"
        )
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_eigenvalues(self):
        result = self.results["x_cp2k_stress_tensor_eigenvalues"]
        expected_result = convert_unit(np.array([7.77641809, 7.77641797, 7.77641485]), "GPa")
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_eigenvectors(self):
        result = self.results["x_cp2k_stress_tensor_eigenvectors"]
        expected_result = np.array([
            [0.00094549, -0.79967815, 0.60042815],
            [-0.70749682, 0.42379757, 0.56554741],
            [0.70671590, 0.42533573, 0.56536905],
        ])
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_determinant(self):
        result = self.results["x_cp2k_stress_tensor_determinant"]
        expected_result = convert_unit(4.70260626E+02, "GPa^3")
        self.assertTrue(np.array_equal(result, expected_result))

    def test_stress_tensor_one_third_of_trace(self):
        result = self.results["x_cp2k_stress_tensor_one_third_of_trace"]
        expected_result = convert_unit(7.77641697E+00, "GPa")
        self.assertTrue(np.array_equal(result, expected_result))

#===============================================================================
if __name__ == '__main__':
    pass
    logger = logging.getLogger("cp2kparser")
    logger.setLevel(logging.ERROR)

    suites = []
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestErrors))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestXCFunctional))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestEnergyForce))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestStressTensorMethods))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSelfInteractionCorrectionMethod))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestConfigurationPeriodicDimensions))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestSCFConvergence))
    suites.append(unittest.TestLoader().loadTestsFromTestCase(TestForceFiles))
    alltests = unittest.TestSuite(suites)
    unittest.TextTestRunner(verbosity=0).run(alltests)