import numpy as np

temps = [0.000000, 100.000000, 200.000000, 300.000000, 400.000000, 500.000000,
         600.000000, 700.000000, 800.000000, 900.000000]
td_ref = [
    [0.00571624, 0.00571624, 0.00571624, 0.00403776, 0.00403776, 0.00403776],
    [0.00877353, 0.00877353, 0.00877353, 0.00654962, 0.00654962, 0.00654962],
    [0.01513305, 0.01513305, 0.01513305, 0.01151749, 0.01151749, 0.01151749],
    [0.02198817, 0.02198817, 0.02198817, 0.01681392, 0.01681392, 0.01681392],
    [0.02898107, 0.02898107, 0.02898107, 0.02220032, 0.02220032, 0.02220032],
    [0.03603064, 0.03603064, 0.03603064, 0.02762357, 0.02762357, 0.02762357],
    [0.04310888, 0.04310888, 0.04310888, 0.03306543, 0.03306543, 0.03306543],
    [0.05020359, 0.05020359, 0.05020359, 0.03851798, 0.03851798, 0.03851798],
    [0.05730864, 0.05730864, 0.05730864, 0.04397723, 0.04397723, 0.04397723],
    [0.06442059, 0.06442059, 0.06442059, 0.04944096, 0.04944096, 0.04944096]]


def test_ThermalDisplacements(ph_nacl):
    ph_nacl.init_mesh([5, 5, 5],
                      with_eigenvectors=True,
                      is_mesh_symmetry=False,
                      use_iter_mesh=True)
    temperatures = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900]
    ph_nacl.run_thermal_displacements(temperatures=temperatures,
                                      freq_min=1e-2)
    td = ph_nacl.thermal_displacements
    # for td_t in td.thermal_displacements:
    #     print(", ".join(["%10.8f" % x for x in td_t]))
    np.testing.assert_allclose(td_ref, td.thermal_displacements, atol=1e-5)

Add ThermalDisplacements test