mdevaluate/test/test_free_energy_landscape.py

import os

import pytest
import numpy as np

import mdevaluate
import mdevaluate.extra.free_energy_landscape as fel


@pytest.fixture
def trajectory(request):
    return mdevaluate.open(os.path.join(os.path.dirname(__file__), "data/pore"))


def test_get_fel(trajectory):
    test_array = np.array(
        [
            184.4909136,
            233.79320471,
            223.12003988,
            228.49746397,
            200.5626769,
            212.82484221,
            165.10818396,
            170.74123681,
            175.86672931,
        ]
    )

    OW = trajectory.subset(atom_name="OW")
    box = trajectory[0].box
    box_voxels = (np.diag(box) // [0.05, 0.05, 0.05] + [1, 1, 1]) * [0.05, 0.05, 0.05]
    occupation_matrix = fel.occupation_matrix(OW, skip=0, segments=10)
    radius_maxima = 0.05 * 3 ** (1 / 2) + 0.05 / 100
    maxima_matrix = fel.find_maxima(
        occupation_matrix,
        box=box_voxels,
        radius=radius_maxima,
        pore_geometry="cylindrical",
    )
    maxima_matrix = fel.add_distances(maxima_matrix, "cylindrical", np.diag(box) / 2)
    r_bins = np.arange(0, 0.5, 0.02)
    distance_bins = np.arange(1.8, 1.9, 0.01)
    energy_df = fel.distance_resolved_energies(
        maxima_matrix, distance_bins, r_bins, box, "cylindrical", 225
    )
    result = fel.find_energy_maxima(energy_df, r_min=0.05, r_max=0.15)
    assert (np.round(np.array(result["energy"])) == np.round(test_array)).all()