Coverage for cosmolayer / cosmosac / model.py: 97%
33 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-11 14:25 +0000
1"""
2.. module:: cosmolayer.cosmosac.model
3 :synopsis: Model configuration for COSMO-SAC calculations.
5.. moduleauthor:: Charlles Abreu <craabreu@gmail.com>
6"""
8from collections.abc import Callable
9from dataclasses import dataclass, field
11import numpy as np
12from numpy.typing import NDArray
14from .component import Component
15from .constants import (
16 COSMO_SAC_2002_AREA_PER_SEGMENT,
17 COSMO_SAC_2002_AVERAGING_RADIUS,
18 COSMO_SAC_2002_EXPONENTS,
19 COSMO_SAC_2002_F_DECAY,
20 COSMO_SAC_2002_SIGMA_0,
21 COSMO_SAC_2010_AREA_PER_SEGMENT,
22 COSMO_SAC_2010_AVERAGING_RADIUS,
23 COSMO_SAC_2010_EXPONENTS,
24 COSMO_SAC_2010_F_DECAY,
25 COSMO_SAC_2010_SIGMA_0,
26)
27from .interaction_matrices import (
28 create_cosmo_sac_2002_matrix,
29 create_cosmo_sac_2010_matrices,
30)
31from .mixture import Mixture
34@dataclass(frozen=True)
35class Model:
36 r"""Immutable configuration for a COSMO-SAC model variant.
38 Bundles all model-specific parameters into a single object, ensuring that
39 components and mixtures are always created with consistent settings. Use the
40 pre-built :data:`CosmoSac2002Model` and :data:`CosmoSac2010Model` singletons
41 for standard calculations, or construct a custom instance for research purposes.
43 Parameters
44 ----------
45 min_sigma : float
46 Minimum screening charge density in e/Ų.
47 max_sigma : float
48 Maximum screening charge density in e/Ų.
49 num_points : int
50 Number of discrete points in the sigma grid.
51 area_per_segment : float
52 Reference surface area of a single segment in Ų.
53 averaging_radius : float
54 Effective radius for distance-weighted sigma averaging in Å.
55 f_decay : float
56 Decay factor for exponential distance weighting in the sigma averaging
57 procedure.
58 sigma_0 : float or None
59 Standard deviation of the Gaussian probability of a segment to form a
60 hydrogen bond in e/Ų. Set to ``None`` to disable hydrogen-bond
61 splitting (all surface area is assigned to the NHB class).
62 merge_profiles : bool
63 Whether segment-group profiles (NHB, OH, OT) should be merged into a
64 single distribution by default when computing probabilities.
65 temperature_exponents : tuple[int, ...]
66 Exponents applied to the temperature in the interaction energy expression.
67 Must have the same length as the tuple returned by
68 ``interaction_matrix_generator``.
69 interaction_matrix_generator : Callable
70 Function ``f(temperature) -> tuple[NDArray, ...]`` that produces the
71 dimensionless interaction energy matrices ΔW/(RT) at a given temperature.
73 Examples
74 --------
75 Using a pre-built model to create a component:
77 >>> from importlib.resources import files
78 >>> from cosmolayer.cosmosac.model import CosmoSac2002Model
79 >>> path = files("cosmolayer.data") / "C=C(N)O.cosmo"
80 >>> component = CosmoSac2002Model.create_component(path.read_text())
81 >>> component.area
82 97.34554...
84 Inspecting model parameters:
86 >>> CosmoSac2002Model.area_per_segment
87 7.5
88 >>> CosmoSac2002Model.merge_profiles
89 True
90 >>> CosmoSac2002Model.temperature_exponents
91 (1,)
92 """
94 min_sigma: float
95 max_sigma: float
96 num_points: int
97 area_per_segment: float
98 averaging_radius: float
99 f_decay: float
100 sigma_0: float | None
101 merge_profiles: bool
102 temperature_exponents: tuple[int, ...]
103 interaction_matrix_generator: Callable[[float], tuple[NDArray[np.float64], ...]] = (
104 field(repr=False)
105 )
107 def create_interaction_matrices(
108 self, temperature: float
109 ) -> tuple[NDArray[np.float64], ...]:
110 """Create the interaction energy matrices at a given temperature.
112 Parameters
113 ----------
114 temperature : float
115 Temperature in Kelvin.
117 Returns
118 -------
119 tuple[NDArray[np.float64], ...]
120 Dimensionless interaction energy matrices ΔW/(RT). The number of
121 matrices matches the length of :attr:`temperature_exponents`.
123 Examples
124 --------
125 >>> from cosmolayer.cosmosac import CosmoSac2002Model, CosmoSac2010Model
127 COSMO-SAC 2002 produces a single interaction matrix:
129 >>> matrices = CosmoSac2002Model.create_interaction_matrices(298.15)
130 >>> len(matrices)
131 1
132 >>> matrix = matrices[0]
133 >>> matrix.shape
134 (51, 51)
135 >>> print(matrix.min() < 0) # H-bonding can be favorable (negative)
136 True
137 >>> print(matrix.max() > 0) # Misfit interactions are unfavorable
138 True
140 Plotting the COSMO-SAC 2002 interaction matrix:
142 .. plot::
143 :context: close-figs
145 >>> from cosmolayer.cosmosac import CosmoSac2002Model
146 >>> from matplotlib import pyplot as plt
147 >>> matrices = CosmoSac2002Model.create_interaction_matrices(298.15)
148 >>> fig, ax = plt.subplots(figsize=(8, 6))
149 >>> im = ax.imshow(matrices[0], cmap="Spectral")
150 >>> _ = fig.colorbar(im, ax=ax, label="ΔW/(RT)")
151 >>> fig.tight_layout()
153 COSMO-SAC 2010 produces two matrices (one per temperature exponent):
155 >>> matrices = CosmoSac2010Model.create_interaction_matrices(298.15)
156 >>> len(matrices)
157 2
158 >>> all(m.shape == (153, 153) for m in matrices)
159 True
161 Plotting the COSMO-SAC 2010 interaction matrix:
163 .. plot::
164 :context: close-figs
166 >>> from cosmolayer.cosmosac import CosmoSac2010Model
167 >>> from matplotlib import pyplot as plt
168 >>> matrices = CosmoSac2010Model.create_interaction_matrices(298.15)
169 >>> fig, ax = plt.subplots(figsize=(8, 6))
170 >>> im = ax.imshow(sum(matrices), cmap="Spectral")
171 >>> _ = fig.colorbar(im, ax=ax, label="ΔW/(RT)")
172 >>> fig.tight_layout()
174 """
175 return self.interaction_matrix_generator(temperature)
177 def create_component(
178 self,
179 cosmo_string: str,
180 ) -> Component:
181 """Create a :class:`~cosmolayer.cosmosac.component.Component` consistent
182 with this model.
184 Parameters
185 ----------
186 cosmo_string : str
187 Contents of a COSMO output file.
189 Returns
190 -------
191 Component
192 Molecular component configured with the model's parameters.
194 Examples
195 --------
196 >>> from importlib.resources import files
197 >>> from cosmolayer.cosmosac.model import CosmoSac2010Model
198 >>> path = files("cosmolayer.data") / "C=C(N)O.cosmo"
199 >>> component = CosmoSac2010Model.create_component(path.read_text())
200 >>> component.area
201 97.34554...
202 >>> component.probabilities.shape
203 (153,)
204 """
205 return Component(
206 cosmo_string,
207 min_sigma=self.min_sigma,
208 max_sigma=self.max_sigma,
209 num_points=self.num_points,
210 averaging_radius=self.averaging_radius,
211 f_decay=self.f_decay,
212 sigma_0=self.sigma_0,
213 merge_profiles=self.merge_profiles,
214 )
216 def create_mixture(
217 self,
218 components: dict[str, str],
219 ) -> "Mixture":
220 """Create a :class:`~cosmolayer.cosmosac.mixture.Mixture` consistent with
221 this model.
223 Parameters
224 ----------
225 components : dict[str, str]
226 Dictionary mapping component names to COSMO strings (contents of
227 COSMO output files).
229 Returns
230 -------
231 Mixture
232 Mixture configured with the model's parameters.
234 Examples
235 --------
236 >>> from importlib.resources import files
237 >>> from cosmolayer.cosmosac.model import CosmoSac2002Model, CosmoSac2010Model
239 Creating a mixture with the COSMO-SAC 2002 model:
241 >>> source = files("cosmolayer.data")
242 >>> components = {
243 ... "1-aminoethenol": (source / "C=C(N)O.cosmo").read_text(),
244 ... "2-aminoethanol": (source / "NCCO.cosmo").read_text(),
245 ... }
246 >>> mixture = CosmoSac2002Model.create_mixture(components)
247 >>> len(mixture)
248 2
249 >>> mixture.interaction_matrices(298.15)[0].shape
250 (51, 51)
252 Creating a mixture with the COSMO-SAC 2010 model:
254 >>> mixture = CosmoSac2010Model.create_mixture(components)
255 >>> len(mixture)
256 2
257 >>> matrices = mixture.interaction_matrices(298.15)
258 >>> len(matrices)
259 2
260 >>> all(m.shape == (153, 153) for m in matrices)
261 True
262 """
264 return Mixture(
265 components,
266 min_sigma=self.min_sigma,
267 max_sigma=self.max_sigma,
268 num_points=self.num_points,
269 area_per_segment=self.area_per_segment,
270 averaging_radius=self.averaging_radius,
271 f_decay=self.f_decay,
272 sigma_0=self.sigma_0,
273 merge_profiles=self.merge_profiles,
274 interaction_matrix_generator=self.interaction_matrix_generator,
275 temperature_exponents=self.temperature_exponents,
276 )
278 @property
279 def num_segment_types(self) -> int:
280 """Number of segment types."""
281 if self.merge_profiles:
282 return self.num_points
283 return 3 * self.num_points
286CosmoSac2002Model = Model(
287 min_sigma=-0.025,
288 max_sigma=0.025,
289 num_points=51,
290 area_per_segment=COSMO_SAC_2002_AREA_PER_SEGMENT,
291 averaging_radius=COSMO_SAC_2002_AVERAGING_RADIUS,
292 f_decay=COSMO_SAC_2002_F_DECAY,
293 sigma_0=COSMO_SAC_2002_SIGMA_0,
294 merge_profiles=True,
295 temperature_exponents=COSMO_SAC_2002_EXPONENTS,
296 interaction_matrix_generator=lambda T: (create_cosmo_sac_2002_matrix(T),),
297)
299CosmoSac2010Model = Model(
300 min_sigma=-0.025,
301 max_sigma=0.025,
302 num_points=51,
303 area_per_segment=COSMO_SAC_2010_AREA_PER_SEGMENT,
304 averaging_radius=COSMO_SAC_2010_AVERAGING_RADIUS,
305 f_decay=COSMO_SAC_2010_F_DECAY,
306 sigma_0=COSMO_SAC_2010_SIGMA_0,
307 merge_profiles=False,
308 temperature_exponents=COSMO_SAC_2010_EXPONENTS,
309 interaction_matrix_generator=create_cosmo_sac_2010_matrices,
310)