Skip to content

hessian #

Functionality for computing the Hessian, closely based on https://github.com/qubekit/QUBEKit/blob/bdccc28d1ad8c0fb4eeb01467bfe8365396c411c/qubekit/engines/openmm.py#L174.

Functions:

  • calculate_hessian –

    Using finite displacement calculate the hessian matrix of the molecule

_get_forces #

_get_forces(
    simulation: Simulation, positions: ndarray
) -> ndarray

Compute the forces on the system at the given positions

Source code in presto/hessian.py 
15
16
17
18
19
20
21
22
23
24
25
def _get_forces(simulation: Simulation, positions: np.ndarray) -> np.ndarray:
    """
    Compute the forces on the system at the given positions
    """
    simulation.context.setPositions(positions * _LENGTH_UNIT)
    forces: np.ndarray = (
        simulation.context.getState(getForces=True)
        .getForces(asNumpy=True)
        .value_in_unit(_FORCE_UNIT)
    )
    return forces

calculate_hessian #

calculate_hessian(
    simulation: Simulation,
    input_coords: Quantity,
    finite_step: Quantity = 0.0005291772 * nanometers,
) -> Quantity

Using finite displacement calculate the hessian matrix of the molecule using a seminumerical scheme. See https://mattermodeling.stackexchange.com/questions/10358/numerical-evaluation-of-hessian for a brief description of the method.

Args: simulation (Simulation): The OpenMM Simulation object. input_coords (unit.Quantity): The coordinates of the molecule. finite_step (unit.Quantity): The finite step size for displacement.

Returns: unit.Quantity: The Hessian matrix in appropriate units.

Source code in presto/hessian.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def calculate_hessian(
    simulation: Simulation,
    input_coords: unit.Quantity,
    finite_step: unit.Quantity = 0.0005291772 * unit.nanometers,
) -> unit.Quantity:
    """
    Using finite displacement calculate the hessian matrix of the molecule
    using a seminumerical scheme. See https://mattermodeling.stackexchange.com/questions/10358/numerical-evaluation-of-hessian
    for a brief description of the method.

    Args:
        simulation (Simulation): The OpenMM Simulation object.
        input_coords (unit.Quantity): The coordinates of the molecule.
        finite_step (unit.Quantity): The finite step size for displacement.

    Returns:
        unit.Quantity: The Hessian matrix in appropriate units.
    """

    # Perform all operations in _LENGTH_UNIT
    input_coords = input_coords.value_in_unit(_LENGTH_UNIT)
    finite_step = finite_step.value_in_unit(_LENGTH_UNIT)
    n_atoms = simulation.topology.getNumAtoms()

    # Check we have coordinates corresponding to a single conformation
    if input_coords.shape != (n_atoms, 3):
        raise ValueError(
            f"Input coordinates must have shape ({n_atoms}, 3), got {input_coords.shape}"
        )

    # Calculate the full Hessian elements with finite difference
    hessian = np.zeros((3 * n_atoms, 3 * n_atoms))

    for j in range(0, n_atoms):
        for _dimension, offset_j in _DIMENSION_TO_OFFSET.items():
            coords = input_coords.copy()
            # Get gi + hej
            coords[j, offset_j] += finite_step
            gi_plus = -_get_forces(simulation, coords)
            # Get gi - hej
            coords[j, offset_j] -= 2 * finite_step
            gi_minus = -_get_forces(simulation, coords)

            for i in range(0, n_atoms):
                for _dimension, offset_i in _DIMENSION_TO_OFFSET.items():
                    # Compute hessian element
                    hessian[i * 3 + offset_i, j * 3 + offset_j] = (
                        gi_plus[i, offset_i] - gi_minus[i, offset_i]
                    ) / (2 * finite_step)

    # Symmetrise the Hessian by averaging with its transpose
    hessian = 0.5 * (hessian + hessian.T)

    return hessian * _HESSIAN_UNIT