Coverage for cosmolayer / cosmosac / mixture.py: 96%
73 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 14:25 +0000
1from collections import OrderedDict
2from collections.abc import Callable
4import numpy as np
5from numpy.typing import NDArray
7from .component import Component
8from .constants import (
9 COSMO_SAC_2010_AREA_PER_SEGMENT,
10 COSMO_SAC_2010_AVERAGING_RADIUS,
11 COSMO_SAC_2010_EXPONENTS,
12 COSMO_SAC_2010_F_DECAY,
13 COSMO_SAC_2010_SIGMA_0,
14)
15from .interaction_matrices import create_cosmo_sac_2010_matrices
18class Mixture:
19 """Mixture of molecular components for COSMO-SAC calculations.
21 This class manages a collection of molecular components, each defined by
22 a COSMO output file from quantum mechanical calculations.
24 .. note::
25 The default parameters correspond to the COSMO-SAC 2010 model :cite:`Bell2020`.
27 Parameters
28 ----------
29 components : dict[str, str]
30 Dictionary mapping component names to COSMO strings (i.e., contents of COSMO
31 output files from quantum mechanical calculations).
32 min_sigma : float, optional
33 Minimum screening charge density in e/Ų. Default is -0.025 e/Ų.
34 max_sigma : float, optional
35 Maximum screening charge density in e/Ų. Default is 0.025 e/Ų.
36 num_points : int, optional
37 Number of discrete points in the sigma profile. Default is 51.
38 area_per_segment : float, optional
39 Reference area in Ų. Default is 7.25 Ų.
40 averaging_radius : float, optional
41 Effective radius for distance-weighted sigma averaging in Å.
42 Default is √(7.25 / π) Å.
43 f_decay : float, optional
44 Decay factor for exponential distance weighting. Default is 3.57.
45 sigma_0 : float or None, optional
46 Standard deviation of the Gaussian probability of a segment to form
47 a hydrogen bond in e/Ų. Set to ``None`` to disable hydrogen-bond
48 splitting (all surface area is assigned to the NHB class).
49 Default is 0.007 e/Ų.
50 merge_profiles : bool, optional
51 Whether to merge segment groups (NHB, OH, OT) into a single profile
52 when accessing :attr:`probabilities` and :attr:`sigma_profiles`.
53 Default is False.
54 regularize : float, optional
55 Minimum value for clipping probabilities. Default is 1e-10.
56 interaction_matrix_generator : Callable, optional
57 Function to generate the interaction matrix for the mixture at a given
58 temperature. Default is :func:`create_cosmo_sac_2010_matrices` with default
59 parameters.
60 temperature_exponents : tuple[float, ...], optional
61 Temperature exponents for the interaction matrices. Must be the same length as
62 the tuple returned by ``interaction_matrix_generator``.
63 Default is (1, 3).
65 Raises
66 ------
67 ValueError
68 If no components are provided.
69 FileNotFoundError
70 If any of the specified files do not exist.
72 Examples
73 --------
74 >>> from importlib.resources import files
75 >>> from cosmolayer.cosmosac import Mixture
76 >>> source = files("cosmolayer.data")
77 >>> components = {
78 ... "1-aminoethenol": (source / "C=C(N)O.cosmo").read_text(),
79 ... "2-aminoethanol": (source / "NCCO.cosmo").read_text(),
80 ... }
81 >>> mixture = Mixture(components)
82 >>> len(mixture)
83 2
84 >>> mixture["1-aminoethenol"].area
85 97.34554...
86 >>> mixture["2-aminoethanol"].area
87 103.51765...
88 >>> mixture.component_names
89 ('1-aminoethenol', '2-aminoethanol')
90 >>> areas = mixture.areas
91 >>> areas.shape
92 (2,)
93 >>> float(areas[0])
94 97.34554...
95 """
97 def __init__( # noqa: PLR0913
98 self,
99 components: dict[str, str],
100 min_sigma: float = -0.025, # e/Ų
101 max_sigma: float = 0.025, # e/Ų
102 num_points: int = 51,
103 area_per_segment: float = COSMO_SAC_2010_AREA_PER_SEGMENT, # Ų
104 averaging_radius: float = COSMO_SAC_2010_AVERAGING_RADIUS, # Å
105 f_decay: float = COSMO_SAC_2010_F_DECAY,
106 sigma_0: float | None = COSMO_SAC_2010_SIGMA_0, # e/Ų
107 merge_profiles: bool = False,
108 interaction_matrix_generator: Callable[
109 [float], tuple[NDArray[np.float64], ...]
110 ] = create_cosmo_sac_2010_matrices,
111 temperature_exponents: tuple[float, ...] = COSMO_SAC_2010_EXPONENTS,
112 ) -> None:
113 if not components:
114 raise ValueError("At least one component must be provided.")
116 self._components: OrderedDict[str, Component] = OrderedDict()
117 self._min_sigma = min_sigma
118 self._max_sigma = max_sigma
119 self._num_points = num_points
120 self._area_per_segment = area_per_segment
121 self._averaging_radius = averaging_radius
122 self._f_decay = f_decay
123 self._sigma_0 = sigma_0
124 self._merge_profiles = merge_profiles
125 self._interaction_matrix_generator = interaction_matrix_generator
126 self._temperature_exponents = temperature_exponents
128 for name, cosmo_string in components.items():
129 self.add_component(name, cosmo_string)
131 def __len__(self) -> int:
132 """Return the number of components in the mixture."""
133 return len(self._components)
135 def __getitem__(self, name: str) -> Component:
136 """Get a component by name.
138 Parameters
139 ----------
140 name : str
141 Component name.
143 Returns
144 -------
145 Component
146 The requested component.
148 Raises
149 ------
150 KeyError
151 If component name not found.
152 """
153 return self._components[name]
155 def _create_component(self, cosmo_string: str) -> Component:
156 return Component(
157 cosmo_string,
158 min_sigma=self._min_sigma,
159 max_sigma=self._max_sigma,
160 num_points=self._num_points,
161 averaging_radius=self._averaging_radius,
162 f_decay=self._f_decay,
163 sigma_0=self._sigma_0,
164 merge_profiles=self._merge_profiles,
165 )
167 def add_component(self, name: str, cosmo_string: str) -> None:
168 """Add a component to the mixture.
170 Parameters
171 ----------
172 name : str
173 Component name.
174 cosmo_string : str
175 Contents of a COSMO output file.
176 """
177 self._components[name] = self._create_component(cosmo_string)
179 def remove_component(self, name: str) -> None:
180 """Remove a component from the mixture.
182 Parameters
183 ----------
184 name : str
185 Component name.
186 """
187 del self._components[name]
189 def replace_component(
190 self, old_name: str, new_name: str, cosmo_string: str
191 ) -> None:
192 """Replace a component in the mixture.
194 The new name must not already exist in the mixture, unless it is the same as
195 the old name. In this case, the component data is updated using the new COSMO
196 string.
198 Parameters
199 ----------
200 old_name : str
201 Name of the component to replace.
202 new_name : str
203 Name of the new component.
204 cosmo_string : str
205 Contents of the new component's COSMO output file.
207 Raises
208 ------
209 ValueError
210 If the old name is not found in the mixture.
211 If the new name already exists in the mixture and is not the same as the
212 old name.
213 """
214 if old_name not in self._components:
215 raise ValueError(f"Component {old_name} not found in mixture.")
216 if new_name != old_name and new_name in self._components:
217 raise ValueError(f"Component {new_name} already exists in mixture.")
218 if new_name == old_name:
219 self._components[new_name] = self._create_component(cosmo_string)
220 else:
221 old_components = self._components
222 self._components = OrderedDict()
223 for name, component in old_components.items():
224 if name == old_name:
225 self._components[new_name] = self._create_component(cosmo_string)
226 else:
227 self._components[name] = component
229 @property
230 def area_per_segment(self) -> float:
231 """Reference area per segment used by the COSMO-SAC model, in Ų."""
232 return self._area_per_segment
234 @property
235 def merge_profiles(self) -> bool:
236 """Whether segment groups (NHB, OH, OT) are merged for :attr:`sigma_profiles`
237 and :attr:`probabilities`."""
238 return self._merge_profiles
240 @property
241 def component_names(self) -> tuple[str, ...]:
242 """Names of all components in the order they were provided."""
243 return tuple(self._components.keys())
245 @property
246 def areas(self) -> NDArray[np.float64]:
247 """Cavity surface areas for all components in Ų. Shape: (n_components,)."""
248 return np.array([component.area for component in self._components.values()])
250 @property
251 def volumes(self) -> NDArray[np.float64]:
252 """Cavity volumes for all components in ų. Shape: (n_components,)."""
253 return np.array([component.volume for component in self._components.values()])
255 @property
256 def probabilities(self) -> NDArray[np.float64]:
257 """Normalized segment-type probabilities for each component.
259 Stack of each component's :attr:`Component.probabilities`. Shape is
260 ``(n_components, num_points)`` if :attr:`merge_profiles` is True, else
261 ``(n_components, 3*num_points)``.
263 Returns
264 -------
265 np.ndarray
266 Probabilities; each row sums to 1.0.
267 """
268 return np.stack(
269 [component.probabilities for component in self._components.values()],
270 axis=0,
271 )
273 @property
274 def sigma_profiles(self) -> NDArray[np.float64]:
275 """Per-component sigma profiles (surface area vs. charge density), in Ų.
277 Stack of each component's :attr:`Component.sigma_profile`. Shape is
278 ``(n_components, num_points)`` when :attr:`merge_profiles` is True,
279 ``(n_components, 3, num_points)`` when False (NHB, OH, OT).
281 Returns
282 -------
283 np.ndarray
284 Sigma profile array; last dimension is the sigma grid.
285 """
286 return np.stack(
287 [component.sigma_profile for component in self._components.values()],
288 axis=0,
289 )
291 def interaction_matrices(
292 self, temperature: float
293 ) -> tuple[NDArray[np.float64], ...]:
294 """COSMO-SAC interaction matrices for the mixture at the given temperature.
296 Parameters
297 ----------
298 temperature : float
299 Temperature in K; used to scale the matrices.
301 Returns
302 -------
303 tuple of np.ndarray
304 Matrices used in the COSMO-SAC activity coefficient calculation.
305 Length and shapes match the generator (e.g. sigma–sigma and
306 sigma'–sigma' for the 2010 model).
307 """
308 return self._interaction_matrix_generator(temperature)
310 @property
311 def temperature_exponents(self) -> tuple[float, ...]:
312 """Exponents used to scale each interaction matrix with temperature.
314 Each entry scales the corresponding matrix from
315 :meth:`interaction_matrices` as T^exponent (e.g. 1 and 3 for the
316 COSMO-SAC 2010 model).
318 Returns
319 -------
320 tuple of float
321 One exponent per interaction matrix.
322 """
323 return self._temperature_exponents