1
"""
2
This class contains a simulation analysis wrapper for use with the OpenFF Evaluator.
3
"""
4

5 38
import logging
6 38
from typing import Any, Dict, List
7

8 38
import numpy as np
9

10 38
from paprika.analysis import fe_calc
11 38
from paprika.restraints import DAT_restraint
12

13 38
logger = logging.getLogger(__name__)
14

15

16 38
class Analyze:
17
    """
18
    The ``Analyze`` class provides a wrapper function around the analysis of simulations.
19
    """
20

21 38
    @classmethod
22 38
    def compute_phase_free_energy(
23
        cls,
24
        phase: str,
25
        restraints: List[DAT_restraint],
26
        windows_directory: str,
27
        topology_name: str,
28
        trajectory_mask: str = "*.dcd",
29
        bootstrap_cycles: int = 1000,
30
        analysis_method: str = "ti-block",
31
    ) -> Dict[str, Any]:
32
        """Computes the free energy of a specified phase of an APR calculation.
33

34
        Parameters
35
        ----------
36
        phase
37
            The phase to analyze. This should be one of ``'attach'``, ``'pull'`` or
38
            ``'release'``.
39
        restraints
40
            A list of the restraints which were used as part of the phase being
41
            analysed.
42
        windows_directory
43
            The directory which contains the final trajectories and topologies
44
            from the simulated phase.
45
        topology_name
46
            The expected file name (not path) of the coordinate file which contains
47
            topological information about the system.
48
        trajectory_mask
49
            The pattern to use when searching for the simulated trajectories.
50
        bootstrap_cycles
51
            The number of bootstrap iterations to perform.
52
        analysis_method
53
            The analysis method to use.
54

55
        Returns
56
        -------
57
            The computed free energies (and their uncertainties) in units of kcal / mol
58
        """
59 0
        analysis = fe_calc()
60

61 0
        analysis.prmtop = topology_name
62 0
        analysis.trajectory = trajectory_mask
63 0
        analysis.path = windows_directory
64

65 0
        analysis.restraint_list = restraints
66

67 0
        analysis.methods = [analysis_method]
68 0
        analysis.bootcycles = bootstrap_cycles
69

70 0
        analysis.collect_data(single_prmtop=False)
71 0
        analysis.compute_free_energy(phases=[phase])
72

73 0
        return analysis.results
74

75 38
    @classmethod
76 38
    def compute_ref_state_work(
77
        cls,
78
        temperature: float,
79
        guest_restraints: List[DAT_restraint],
80
    ) -> float:
81
        """Computes the reference state work of the attach phase.
82

83
        Parameters
84
        ----------
85
        temperature
86
            The temperature at which the calculation was performed in units of kelvin.
87
        guest_restraints
88
            The guest restraints which were applied during the pull phase.
89

90
        Returns
91
        -------
92
            The reference state work in units of kcal / mol.
93
        """
94

95 38
        analysis = fe_calc()
96 38
        analysis.temperature = temperature
97

98 38
        distance_restraint = next(
99
            (
100
                restraint
101
                for restraint in guest_restraints
102
                if "DM" in restraint.mask1
103
                and restraint.mask2 is not None
104
                and restraint.mask3 is None
105
                and restraint.mask4 is None
106
            ),
107
            None,
108
        )
109

110 38
        theta_restraint = next(
111
            (
112
                restraint
113
                for restraint in guest_restraints
114
                if "DM" in restraint.mask1
115
                and "DM" in restraint.mask2
116
                and restraint.mask3 is not None
117
                and restraint.mask4 is None
118
            ),
119
            None,
120
        )
121 38
        beta_restraint = next(
122
            (
123
                restraint
124
                for restraint in guest_restraints
125
                if "DM" in restraint.mask1
126
                and "DM" not in restraint.mask2
127
                and restraint.mask3 is not None
128
                and restraint.mask4 is None
129
            ),
130
            None,
131
        )
132

133 38
        if not distance_restraint or not theta_restraint or not beta_restraint:
134

135 0
            raise RuntimeError(
136
                "Could not determine the r, θ, or β restraint for computing the "
137
                "analytic release step."
138
            )
139

140 38
        analysis.compute_ref_state_work(
141
            [distance_restraint, theta_restraint, None, None, beta_restraint, None]
142
        )
143

144 38
        return analysis.results["ref_state_work"]
145

146 38
    @classmethod
147 38
    def symmetry_correction(
148
        cls,
149
        n_microstates: int,
150
        temperature: float,
151
    ) -> float:
152
        """Computes the free energy corrections to apply to symmetrical guest
153
        when the guest is restrained to just one of several possible symmetrical
154
        configurations (e.g. butane restrained in a cyclic host).
155

156
        Parameters
157
        ----------
158
        temperature
159
            The temperature at which the calculation was performed in units of kelvin.
160
        n_microstates
161
            The number of different symmetrical microstates that the guest can exist
162
            in (e.g. for butane this is two).
163

164
        Returns
165
        -------
166
            The symmetry correction to apply in units of kcal / mol.
167
        """
168

169 38
        assert n_microstates > 0
170

171 38
        if n_microstates == 1:
172 38
            return 0.0
173

174 38
        return -temperature * 0.001987204258640832 * np.log(n_microstates)

Read our documentation on viewing source code .

Loading