Coverage for cvpack/attraction_strength.py: 98%
48 statements
« prev ^ index » next coverage.py v7.5.1, created at 2024-05-09 16:14 +0000
1"""
2.. class:: AttractionStrength
3 :platform: Linux, MacOS, Windows
4 :synopsis: The strength of the attraction between two atom groups
6.. classauthor:: Charlles Abreu <craabreu@gmail.com>
8"""
10import typing as t
12import numpy as np
13import openmm
14from openmm import unit as mmunit
16from .collective_variable import CollectiveVariable
17from .units import Quantity, ScalarQuantity
18from .utils import evaluate_in_context
20ONE_4PI_EPS0 = 138.93545764438198
23class AttractionStrength(CollectiveVariable, openmm.CustomNonbondedForce):
24 r"""
25 The strength of the attraction between two atom groups:
27 .. math::
29 S_{\rm attr}({\bf r}) = S_{1,2}({\bf r}) = &-\frac{1}{E_{\rm ref}} \Bigg[
30 \sum_{i \in {\bf g}_1}
31 \sum_{\substack{j \in {\bf g}_2 \\ j \neq i}}
32 \epsilon_{ij} u_{\rm disp}\left(
33 \frac{\|{\bf r}_i - {\bf r}_j\|}{\sigma_{ij}}
34 \right) \\
35 &+\sum_{i \in {\bf g}_1}
36 \sum_{\substack{j \in {\bf g}_2 \\ q_jq_i < 0}}
37 \frac{q_i q_j}{4 \pi \varepsilon_0 r_{\rm c}} u_{\rm elec}\left(
38 \frac{\|{\bf r}_i - {\bf r}_j\|}{r_{\rm c}}
39 \right) \Bigg]
41 where :math:`{\bf g}_1` and :math:`{\bf g}_2` are the two atom groups,
42 :math:`r_{\rm c}` is the cutoff distance, :math:`\varepsilon_0` is the
43 permittivity of empty space, :math:`q_i` is the charge of atom :math:`i`, and
44 :math:`E_{\rm ref}` is a reference value (in energy units per mole).
45 The Lennard-Jones parameters are given by the Lorentz-Berthelot mixing rule, i.e.
46 :math:`\epsilon_{ij} = \sqrt{\epsilon_i \epsilon_j}`, and
47 :math:`\sigma_{ij} = (\sigma_i + \sigma_j)/2`.
49 Optionally, one can provide a third atom group :math:`{\bf g}_3` to contrast
50 the attraction strength between :math:`{\bf g}_1` and :math:`{\bf g}_2` with
51 that between :math:`{\bf g}_1` and :math:`{\bf g}_3`. One can also provide
52 a scaling factor :math:`\alpha` to balance the contributions of the two
53 interactions. In this case, the collective variable becomes
55 .. math::
57 S_{\rm attr}({\bf r}) = S_{1,2}({\bf r}) - \alpha S_{1,3}({\bf r})
59 .. note::
61 Groups :math:`{\bf g}_1` and :math:`{\bf g}_2` can overlap or even be the
62 same, in which case the collective variable will measure the strength of
63 the self-attraction of :math:`{\bf g}_1`. On the other hand, the contrast
64 group :math:`{\bf g}_3` cannot overlap with neither :math:`{\bf g}_1` nor
65 :math:`{\bf g}_2`.
67 The function :math:`u_{\rm disp}(x)` is a Lennard-Jones-like reduced potential with
68 a highly softened repulsion part, defined as
70 .. math::
72 u_{\rm disp}(x) = 4\left(\frac{1}{y^2} - \frac{1}{y}\right),
73 \quad \text{where} \quad
74 y = |x^6 - 2| + 2
76 The function :math:`u_{\rm elec}(x)` provides a Coulomb-like decay with
77 reaction-field screening, defined as
79 .. math::
81 u_{\rm elec}(x) = \frac{1}{x} + \frac{x^2 - 3}{2}
83 The screening considers a perfect conductor as the surrounding medium
84 :cite:`Correa_2022`.
86 .. note::
88 Only attractive electrostatic interactions are considered (:math:`q_i q_i < 0`),
89 which gives :math:`S_{\rm attr}({\bf r})` a lower bound of zero. The upper
90 bound will depends on the system details, the chosen groups of atoms, and the
91 adopted reference value.
93 The Lennard-Jones parameters, atomic charges, cutoff distance, boundary conditions,
94 as well as whether to use a switching function and its corresponding switching
95 distance, are taken from :openmm:`NonbondedForce` object.
97 .. note::
98 Any non-exclusion exceptions involving an atom in :math:`{\bf g}_1` and an atom
99 in either :math:`{\bf g}_2` or :math:`{\bf g}_3` are turned into exclusions in
100 this collective variable.
102 Parameters
103 ----------
104 group1
105 The first atom group.
106 group2
107 The second atom group.
108 nonbondedForce
109 The :class:`openmm.NonbondedForce` object from which to collect the necessary
110 parameters.
111 contrastGroup
112 An optional third atom group to contrast the attraction strength between
113 :math:`{\bf g}_1` and :math:`{\bf g}_2` with that between :math:`{\bf g}_1`
114 and :math:`{\bf g}_3`.
115 reference
116 A reference value (in energy units per mole) to which the collective variable
117 should be normalized. One can also provide an :OpenMM:`Context` object from
118 which to extract a reference attraction strength. The extracted value will be
119 :math:`S_{1,2}({\bf r})` for :math:`E_{\rm ref} = 1`, regardless of whether
120 `contrastGroup` is provided or not.
121 contrastScaling
122 A scaling factor :math:`\alpha` to balance the contributions of the two
123 interactions. The default is 1.0.
124 name
125 The name of the collective variable.
127 Examples
128 --------
129 >>> import cvpack
130 >>> from openmm import unit
131 >>> from openmmtools import testsystems
132 >>> model = testsystems.HostGuestExplicit()
133 >>> guest = [a.index for a in model.topology.atoms() if a.residue.name == "B2"]
134 >>> host = [a.index for a in model.topology.atoms() if a.residue.name == "CUC"]
135 >>> forces = {f.getName(): f for f in model.system.getForces()}
136 >>> cv1 = cvpack.AttractionStrength(guest, host, forces["NonbondedForce"])
137 >>> cv1.addToSystem(model.system)
138 >>> platform = openmm.Platform.getPlatformByName("Reference")
139 >>> integrator = openmm.VerletIntegrator(1.0 * mmunit.femtoseconds)
140 >>> context = openmm.Context(model.system, integrator, platform)
141 >>> context.setPositions(model.positions)
142 >>> cv1.getValue(context)
143 4912.5... dimensionless
145 >>> water = [a.index for a in model.topology.atoms() if a.residue.name == "HOH"]
146 >>> cv2 = cvpack.AttractionStrength(guest, water, forces["NonbondedForce"])
147 >>> cv2.addToSystem(model.system)
148 >>> context.reinitialize(preserveState=True)
149 >>> cv2.getValue(context)
150 2063.3... dimensionless
152 >>> cv3 = cvpack.AttractionStrength(guest, host, forces["NonbondedForce"], water)
153 >>> cv3.addToSystem(model.system)
154 >>> context.reinitialize(preserveState=True)
155 >>> cv3.getValue(context)
156 2849.17... dimensionless
157 >>> print(cv1.getValue(context) - cv2.getValue(context))
158 2849.17... dimensionless
160 >>> cv4 = cvpack.AttractionStrength(
161 ... guest, host, forces["NonbondedForce"], water, contrastScaling=0.5
162 ... )
163 >>> cv4.addToSystem(model.system)
164 >>> context.reinitialize(preserveState=True)
165 >>> cv4.getValue(context)
166 3880.8... dimensionless
167 >>> 1 * cv1.getValue(context) - 0.5 * cv2.getValue(context)
168 3880.8...
169 """
171 def __init__( # pylint: disable=too-many-locals
172 self,
173 group1: t.Iterable[int],
174 group2: t.Iterable[int],
175 nonbondedForce: openmm.NonbondedForce,
176 contrastGroup: t.Optional[t.Iterable[int]] = None,
177 reference: t.Union[ScalarQuantity, openmm.Context] = Quantity(
178 1.0 * mmunit.kilojoule_per_mole
179 ),
180 contrastScaling: float = 1.0,
181 name: str = "attraction_strength",
182 ) -> None:
183 group1 = list(group1)
184 group2 = list(group2)
185 contrasting = contrastGroup is not None
186 if contrasting:
187 contrastGroup = list(contrastGroup)
188 cutoff = nonbondedForce.getCutoffDistance().value_in_unit(mmunit.nanometers)
189 expression = (
190 f"-(lj + coul){'*sign1*sign2' if contrasting else ''}/ref"
191 "; lj = 4*epsilon*(1/y^2 - 1/y)"
192 f"; coul = {ONE_4PI_EPS0}*q1q2*(1/x + (x^2 - 3)/2)"
193 f"; x = r/{cutoff}"
194 "; y = abs((r/sigma)^6 - 2) + 2"
195 "; q1q2 = min(0, charge1*charge2)"
196 "; epsilon = sqrt(epsilon1*epsilon2)"
197 "; sigma = (sigma1 + sigma2)/2"
198 "; ref = 1"
199 )
200 super().__init__(expression)
201 self.setNonbondedMethod(
202 self.CutoffPeriodic
203 if nonbondedForce.usesPeriodicBoundaryConditions()
204 else self.CutoffNonPeriodic
205 )
206 self.setCutoffDistance(cutoff)
207 for parameter in ("charge", "sigma", "epsilon") + ("sign",) * contrasting:
208 self.addPerParticleParameter(parameter)
209 contrast_atoms = np.zeros(nonbondedForce.getNumParticles(), dtype=bool)
210 if contrasting:
211 contrast_atoms[contrastGroup] = True
212 for atom, in_group in enumerate(contrast_atoms):
213 charge, sigma, epsilon = nonbondedForce.getParticleParameters(atom)
214 if contrasting:
215 sign = -1 if in_group else 1
216 scale = contrastScaling if in_group else 1
217 self.addParticle([charge * scale, sigma, epsilon * scale**2, sign])
218 else:
219 self.addParticle([charge, sigma, epsilon])
220 for exception in range(nonbondedForce.getNumExceptions()):
221 i, j, *_ = nonbondedForce.getExceptionParameters(exception)
222 self.addExclusion(i, j)
223 self.setUseSwitchingFunction(nonbondedForce.getUseSwitchingFunction())
224 self.setSwitchingDistance(nonbondedForce.getSwitchingDistance())
225 self.setUseLongRangeCorrection(False)
226 self.addInteractionGroup(group1, group2)
227 if isinstance(reference, openmm.Context):
228 reference = evaluate_in_context(self, reference)
229 elif mmunit.is_quantity(reference):
230 reference = reference.value_in_unit(mmunit.kilojoule_per_mole)
231 self.setEnergyFunction(expression.replace("ref = 1", f"ref = {reference}"))
232 if contrasting:
233 self.addInteractionGroup(group1, contrastGroup)
235 self._registerCV(
236 name,
237 mmunit.dimensionless,
238 group1=group1,
239 group2=group2,
240 nonbondedForce=nonbondedForce,
241 contrastGroup=contrastGroup,
242 reference=reference,
243 contrastScaling=contrastScaling,
244 )
247AttractionStrength.registerTag("!cvpack.AttractionStrength")