From 0ea657322bce81c7516fbbc789e4ac97a18ff3a9 Mon Sep 17 00:00:00 2001 From: Marvin Friede <51965259+marvinfriede@users.noreply.github.com> Date: Sun, 3 Aug 2025 19:13:27 +0200 Subject: [PATCH 1/2] Batched functorch Hessian (#89) --- .github/workflows/release.yaml | 4 +- .pre-commit-config.yaml | 10 +-- docs/source/modules/data/index.rst | 1 - docs/source/modules/data/radii.rst | 2 - docs/source/modules/typing/builtin.rst | 2 - docs/source/modules/typing/index.rst | 1 - examples/batch.py | 24 ++---- examples/hessian.py | 105 +++++++++++++++++++++++++ examples/single.py | 4 +- setup.cfg | 2 +- src/tad_dftd3/damping/atm.py | 7 +- src/tad_dftd3/damping/rational.py | 6 +- src/tad_dftd3/data/__init__.py | 3 +- src/tad_dftd3/data/r4r2.py | 89 +++++++++++++-------- src/tad_dftd3/data/radii.py | 63 --------------- src/tad_dftd3/data/vdw-d3.pt | Bin 87275 -> 0 bytes src/tad_dftd3/disp.py | 43 ++++++---- src/tad_dftd3/model/c6.py | 102 ++++++++++++++++++++---- src/tad_dftd3/model/weights.py | 11 ++- src/tad_dftd3/reference.py | 5 +- src/tad_dftd3/typing/__init__.py | 1 - src/tad_dftd3/typing/builtin.py | 24 ------ test/conftest.py | 7 +- test/test_disp/test_dftd3.py | 9 ++- test/test_disp/test_disp.py | 19 +++-- test/test_disp/test_special.py | 9 ++- test/test_grad/test_hessian.py | 70 ++++++++++++++++- test/test_grad/test_param.py | 6 +- test/test_grad/test_pos.py | 6 +- test/test_model/samples.py | 3 +- test/test_model/test_c6.py | 52 +++++++++++- test/test_model/test_load.py | 7 +- test/test_model/test_reference.py | 5 +- test/test_model/test_weights.py | 2 +- 34 files changed, 471 insertions(+), 233 deletions(-) delete mode 100644 docs/source/modules/data/radii.rst delete mode 100644 docs/source/modules/typing/builtin.rst create mode 100644 examples/hessian.py delete mode 100644 src/tad_dftd3/data/radii.py delete mode 100644 src/tad_dftd3/data/vdw-d3.pt delete mode 100644 src/tad_dftd3/typing/builtin.py diff --git a/.github/workflows/release.yaml b/.github/workflows/release.yaml index 9028cd2..813bb58 100644 --- a/.github/workflows/release.yaml +++ b/.github/workflows/release.yaml @@ -140,7 +140,7 @@ jobs: merge-multiple: true - name: Publish to Test PyPI - uses: pypa/gh-action-pypi-publish@release/v1 + uses: pypa/gh-action-pypi-publish@76f52bc884231f62b9a034ebfe128415bbaabdfc # v1.12 with: repository-url: https://test.pypi.org/legacy/ @@ -164,4 +164,4 @@ jobs: merge-multiple: true - name: Publish to PyPI - uses: pypa/gh-action-pypi-publish@release/v1 + uses: pypa/gh-action-pypi-publish@76f52bc884231f62b9a034ebfe128415bbaabdfc # v1.12 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 81965b1..667bdbb 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -42,7 +42,7 @@ repos: ] - repo: https://github.com/asottile/pyupgrade - rev: v3.19.1 + rev: v3.20.0 hooks: - id: pyupgrade args: [--py37-plus, --keep-runtime-typing] @@ -52,21 +52,21 @@ repos: hooks: - id: isort name: isort (python) - args: ["--profile", "black", "--filter-files"] + args: ["--profile", "black", "--line-length", "80", "--filter-files"] - repo: https://github.com/psf/black rev: 25.1.0 hooks: - id: black - stages: [pre-commit] + args: ["--line-length", "80"] - repo: https://github.com/woodruffw/zizmor-pre-commit - rev: v1.5.2 + rev: v1.11.0 hooks: - id: zizmor - repo: https://github.com/pre-commit/mirrors-mypy - rev: v1.15.0 + rev: v1.17.1 hooks: - id: mypy pass_filenames: false diff --git a/docs/source/modules/data/index.rst b/docs/source/modules/data/index.rst index 0593c57..d25a5b6 100644 --- a/docs/source/modules/data/index.rst +++ b/docs/source/modules/data/index.rst @@ -5,4 +5,3 @@ .. toctree:: r4r2 - radii diff --git a/docs/source/modules/data/radii.rst b/docs/source/modules/data/radii.rst deleted file mode 100644 index e0e9ded..0000000 --- a/docs/source/modules/data/radii.rst +++ /dev/null @@ -1,2 +0,0 @@ -.. automodule:: tad_dftd3.data.radii - :members: diff --git a/docs/source/modules/typing/builtin.rst b/docs/source/modules/typing/builtin.rst deleted file mode 100644 index 8f79631..0000000 --- a/docs/source/modules/typing/builtin.rst +++ /dev/null @@ -1,2 +0,0 @@ -.. automodule:: tad_dftd3.typing.builtin - :members: diff --git a/docs/source/modules/typing/index.rst b/docs/source/modules/typing/index.rst index 6c66f2e..b9a2dff 100644 --- a/docs/source/modules/typing/index.rst +++ b/docs/source/modules/typing/index.rst @@ -4,6 +4,5 @@ .. toctree:: - builtin d3 pytorch diff --git a/examples/batch.py b/examples/batch.py index 935ba27..e87d8a2 100644 --- a/examples/batch.py +++ b/examples/batch.py @@ -60,30 +60,18 @@ ) ) ref = d3.reference.Reference() -rcov = d3.data.COV_D3[numbers] -rvdw = d3.data.VDW_D3[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] -r4r2 = d3.data.R4R2[numbers] +rcov = mctc.data.COV_D3()[numbers] +rvdw = mctc.data.VDW_PAIRWISE()[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] +r4r2 = d3.data.R4R2()[numbers] param = { "a1": torch.tensor(0.49484001), "s8": torch.tensor(0.78981345), "a2": torch.tensor(5.73083694), } -cn = mctc.ncoord.cn_d3( - numbers, positions, counting_function=mctc.ncoord.exp_count, rcov=rcov -) -weights = d3.model.weight_references(numbers, cn, ref, d3.model.gaussian_weight) -c6 = d3.model.atomic_c6(numbers, weights, ref) -energy = d3.disp.dispersion( - numbers, - positions, - param, - c6, - rvdw, - r4r2, - d3.disp.rational_damping, -) +energy = d3.dftd3(numbers, positions, param) torch.set_printoptions(precision=10) -print(torch.sum(energy, dim=-1)) +print("Expected:", torch.tensor([-0.0014092578, -0.0057840119])) +print("Actual :", torch.sum(energy, dim=-1)) # tensor([-0.0014092578, -0.0057840119]) diff --git a/examples/hessian.py b/examples/hessian.py new file mode 100644 index 0000000..f77b68f --- /dev/null +++ b/examples/hessian.py @@ -0,0 +1,105 @@ +# SPDX-Identifier: CC0-1.0 +from __future__ import annotations + +from typing import Callable + +import tad_mctc as mctc +import torch + +import tad_dftd3 as d3 + +sample1 = dict( + numbers=mctc.convert.symbol_to_number("Pb H H H H Bi H H H".split()), + positions=torch.tensor( + [ + [-0.00000020988889, -4.98043478877778, +0.00000000000000], + [+3.06964045311111, -6.06324400177778, +0.00000000000000], + [-1.53482054188889, -6.06324400177778, -2.65838526500000], + [-1.53482054188889, -6.06324400177778, +2.65838526500000], + [-0.00000020988889, -1.72196703577778, +0.00000000000000], + [-0.00000020988889, +4.77334244722222, +0.00000000000000], + [+1.35700257511111, +6.70626379422222, -2.35039772300000], + [-2.71400388988889, +6.70626379422222, +0.00000000000000], + [+1.35700257511111, +6.70626379422222, +2.35039772300000], + ] + ), +) +sample2 = dict( + numbers=mctc.convert.symbol_to_number( + "C C C C C C I H H H H H S H C H H H".split(" ") + ), + positions=torch.tensor( + [ + [-1.42754169820131, -1.50508961850828, -1.93430551124333], + [+1.19860572924150, -1.66299114873979, -2.03189643761298], + [+2.65876001301880, +0.37736955363609, -1.23426391650599], + [+1.50963368042358, +2.57230374419743, -0.34128058818180], + [-1.12092277855371, +2.71045691257517, -0.25246348639234], + [-2.60071517756218, +0.67879949508239, -1.04550707592673], + [-2.86169588073340, +5.99660765711210, +1.08394899986031], + [+2.09930989272956, -3.36144811062374, -2.72237695164263], + [+2.64405246349916, +4.15317840474646, +0.27856972788526], + [+4.69864865613751, +0.26922271535391, -1.30274048619151], + [-4.63786461351839, +0.79856258572808, -0.96906659938432], + [-2.57447518692275, -3.08132039046931, -2.54875517521577], + [-5.88211879210329, 11.88491819358157, +2.31866455902233], + [-8.18022701418703, 10.95619984550779, +1.83940856333092], + [-5.08172874482867, 12.66714386256482, -0.92419491629867], + [-3.18311711399702, 13.44626574330220, -0.86977613647871], + [-5.07177399637298, 10.99164969235585, -2.10739192258756], + [-6.35955320518616, 14.08073002965080, -1.68204314084441], + ] + ), +) +numbers = mctc.batch.pack( + ( + sample1["numbers"], + sample2["numbers"], + ) +) +positions = mctc.batch.pack( + ( + sample1["positions"], + sample2["positions"], + ) +) + +param = { + "a1": torch.tensor(0.49484001), + "s8": torch.tensor(0.78981345), + "a2": torch.tensor(5.73083694), +} + + +def _energy(numbers: torch.Tensor, positions: torch.Tensor) -> torch.Tensor: + """ + Closure over non-tensor argument `param` for `dftd3` function. + + Returns the energy as a scalar, which is required for Hessian computation + to obtain the correct shape of ``(..., nat, 3, nat, 3)``. + """ + return d3.dftd3(numbers, positions, param).sum(-1) + + +def hessian( + fn: Callable[..., torch.Tensor], argnums: tuple[int] | int = 0 +) -> Callable: + """ + Compute the Hessian using reverse-mode autodiff twice. + (Functorch's `hessian` uses forward and backward mode, but forward is + not implemented for the custom autograd functions in DFT-D3.) + """ + return torch.func.jacrev( + torch.func.jacrev(fn, argnums=argnums), argnums=argnums + ) + + +hess_fn_single = hessian(_energy, argnums=1) +hess_fn_batch = torch.func.vmap(hess_fn_single, in_dims=(0, 0)) + +pos = positions.clone().requires_grad_(True) +hess = hess_fn_batch(numbers, pos) + +print(f"Shape of numbers : {numbers.shape}") +print(f"Shape of positions: {positions.shape}") +print(f"Shape of Hessian : {hess.shape}") diff --git a/examples/single.py b/examples/single.py index 24c2c3f..f0f14c3 100644 --- a/examples/single.py +++ b/examples/single.py @@ -4,7 +4,9 @@ import tad_dftd3 as d3 -numbers = mctc.convert.symbol_to_number(symbols="C C C C N C S H H H H H".split()) +numbers = mctc.convert.symbol_to_number( + symbols="C C C C N C S H H H H H".split() +) positions = torch.tensor( [ [-2.56745685564671, -0.02509985979910, 0.00000000000000], diff --git a/setup.cfg b/setup.cfg index 952b8da..1cdd4a3 100644 --- a/setup.cfg +++ b/setup.cfg @@ -34,7 +34,7 @@ project_urls = packages = find: install_requires = numpy<2 - tad-mctc + tad-mctc==0.5.3 torch>=1.11,<3 python_requires = >=3.8, <3.13 include_package_data = True diff --git a/src/tad_dftd3/damping/atm.py b/src/tad_dftd3/damping/atm.py index 52647a6..cd4d38a 100644 --- a/src/tad_dftd3/damping/atm.py +++ b/src/tad_dftd3/damping/atm.py @@ -33,9 +33,9 @@ import torch from tad_mctc import storch from tad_mctc.batch import real_pairs, real_triples +from tad_mctc.typing import DD, Tensor from .. import defaults -from ..typing import DD, Tensor __all__ = ["dispersion_atm"] @@ -143,7 +143,10 @@ def dispersion_atm( ) ang = torch.where( - mask_triples * (r2ij <= cutoff2) * (r2jk <= cutoff2) * (r2jk <= cutoff2), + mask_triples + * (r2ij <= cutoff2) + * (r2jk <= cutoff2) + * (r2jk <= cutoff2), 0.375 * s / r5 + 1.0 / r3, torch.tensor(0.0, **dd), ) diff --git a/src/tad_dftd3/damping/rational.py b/src/tad_dftd3/damping/rational.py index 5d08877..739e7ef 100644 --- a/src/tad_dftd3/damping/rational.py +++ b/src/tad_dftd3/damping/rational.py @@ -25,12 +25,12 @@ \dfrac{R^n_{\text{AB}}}{R^n_{\text{AB}} + \left( a_1 R_0^{\text{AB}} + a_2 \right)^n} """ -from typing import Dict +from __future__ import annotations import torch +from tad_mctc.typing import DD, Tensor from .. import defaults -from ..typing import DD, Tensor __all__ = ["rational_damping"] @@ -39,7 +39,7 @@ def rational_damping( order: int, distances: Tensor, qq: Tensor, - param: Dict[str, Tensor], + param: dict[str, Tensor], ) -> Tensor: """ Rational damped dispersion interaction between pairs. diff --git a/src/tad_dftd3/data/__init__.py b/src/tad_dftd3/data/__init__.py index f5645e9..f9ce090 100644 --- a/src/tad_dftd3/data/__init__.py +++ b/src/tad_dftd3/data/__init__.py @@ -16,7 +16,6 @@ Atomic data =========== -Data arrays for atomic constants like covalent radii or van-der-Waals radii. +Data arrays for atomic constants. """ from .r4r2 import * -from .radii import * diff --git a/src/tad_dftd3/data/r4r2.py b/src/tad_dftd3/data/r4r2.py index 559567c..ec6f86a 100644 --- a/src/tad_dftd3/data/r4r2.py +++ b/src/tad_dftd3/data/r4r2.py @@ -21,40 +21,65 @@ Also new super heavies Cn, Nh, Fl, Lv, Og and Am-Rg calculated at 4c-PBE/Dyall-AE4Z (Dirac 2022). """ +from __future__ import annotations + import torch __all__ = ["R4R2"] -# fmt: off -r4_over_r2 = torch.tensor([ - 0.0000, # None - 8.0589, 3.4698, # H,He - 29.0974,14.8517,11.8799, 7.8715, 5.5588, 4.7566, 3.8025, 3.1036, # Li-Ne - 26.1552,17.2304,17.7210,12.7442, 9.5361, 8.1652, 6.7463, 5.6004, # Na-Ar - 29.2012,22.3934, # K,Ca - 19.0598,16.8590,15.4023,12.5589,13.4788, # Sc- - 12.2309,11.2809,10.5569,10.1428, 9.4907, # -Zn - 13.4606,10.8544, 8.9386, 8.1350, 7.1251, 6.1971, # Ga-Kr - 30.0162,24.4103, # Rb,Sr - 20.3537,17.4780,13.5528,11.8451,11.0355, # Y- - 10.1997, 9.5414, 9.0061, 8.6417, 8.9975, # -Cd - 14.0834,11.8333,10.0179, 9.3844, 8.4110, 7.5152, # In-Xe - 32.7622,27.5708, # Cs,Ba - 23.1671,21.6003,20.9615,20.4562,20.1010,19.7475,19.4828, # La-Eu - 15.6013,19.2362,17.4717,17.8321,17.4237,17.1954,17.1631, # Gd-Yb - 14.5716,15.8758,13.8989,12.4834,11.4421, # Lu- - 10.2671, 8.3549, 7.8496, 7.3278, 7.4820, # -Hg - 13.5124,11.6554,10.0959, 9.7340, 8.8584, 8.0125, # Tl-Rn - 29.8135,26.3157, # Fr,Ra - 19.1885,15.8542,16.1305,15.6161,15.1226,16.1576,14.6510, # Ac-Am - 14.7178,13.9108,13.5623,13.2326,12.9189,12.6133,12.3142, # Cm-No - 14.8326,12.3771,10.6378, 9.3638, 8.2297, # Lr- - 7.5667, 6.9456, 6.3946, 5.9159, 5.4929, # -Cn - 6.7286, 6.5144,10.9169,10.3600, 9.4723, 8.6641, # Nh-Og -]) # fmt: on -"""Actually calculated r⁴ over r² expectation values.""" -# fmt: on - -R4R2 = torch.sqrt(0.5 * (r4_over_r2 * torch.sqrt(torch.arange(r4_over_r2.shape[0])))) -"""r⁴ over r² expectation values.""" +def R4R2( + dtype: torch.dtype | None = None, device: torch.device | None = None +) -> torch.Tensor: + """ + Returns the r⁴ over r² expectation values as a tensor. + + Parameters + ---------- + dtype : torch.dtype | None, optional + The desired data type of the returned tensor. Defaults to None. + device : torch.device | None, optional + The desired device of the returned tensor. Defaults to None. + + Returns + ------- + Tensor + A tensor containing the r⁴ over r² expectation values. + """ + + # Actually calculated r⁴ over r² expectation values + # fmt: off + _r4_over_r2 = [ + 0.0000, # None + 8.0589, 3.4698, # H,He + 29.0974,14.8517,11.8799, 7.8715, 5.5588, 4.7566, 3.8025, 3.1036, # Li-Ne + 26.1552,17.2304,17.7210,12.7442, 9.5361, 8.1652, 6.7463, 5.6004, # Na-Ar + 29.2012,22.3934, # K,Ca + 19.0598,16.8590,15.4023,12.5589,13.4788, # Sc- + 12.2309,11.2809,10.5569,10.1428, 9.4907, # -Zn + 13.4606,10.8544, 8.9386, 8.1350, 7.1251, 6.1971, # Ga-Kr + 30.0162,24.4103, # Rb,Sr + 20.3537,17.4780,13.5528,11.8451,11.0355, # Y- + 10.1997, 9.5414, 9.0061, 8.6417, 8.9975, # -Cd + 14.0834,11.8333,10.0179, 9.3844, 8.4110, 7.5152, # In-Xe + 32.7622,27.5708, # Cs,Ba + 23.1671,21.6003,20.9615,20.4562,20.1010,19.7475,19.4828, # La-Eu + 15.6013,19.2362,17.4717,17.8321,17.4237,17.1954,17.1631, # Gd-Yb + 14.5716,15.8758,13.8989,12.4834,11.4421, # Lu- + 10.2671, 8.3549, 7.8496, 7.3278, 7.4820, # -Hg + 13.5124,11.6554,10.0959, 9.7340, 8.8584, 8.0125, # Tl-Rn + 29.8135,26.3157, # Fr,Ra + 19.1885,15.8542,16.1305,15.6161,15.1226,16.1576,14.6510, # Ac-Am + 14.7178,13.9108,13.5623,13.2326,12.9189,12.6133,12.3142, # Cm-No + 14.8326,12.3771,10.6378, 9.3638, 8.2297, # Lr- + 7.5667, 6.9456, 6.3946, 5.9159, 5.4929, # -Cn + 6.7286, 6.5144,10.9169,10.3600, 9.4723, 8.6641, # Nh-Og + ] + # fmt: on + + sqrtz = torch.sqrt( + torch.arange(len(_r4_over_r2), device=device, dtype=dtype) + ) + return torch.sqrt( + 0.5 * (torch.tensor(_r4_over_r2, device=device, dtype=dtype) * sqrtz) + ) diff --git a/src/tad_dftd3/data/radii.py b/src/tad_dftd3/data/radii.py deleted file mode 100644 index ac729c3..0000000 --- a/src/tad_dftd3/data/radii.py +++ /dev/null @@ -1,63 +0,0 @@ -# This file is part of tad-dftd3. -# SPDX-Identifier: Apache-2.0 -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -""" -Data: Radii -=========== - -Data arrays for atomic constants like covalent radii or van-der-Waals radii. - -The `vdw_rad_d3` were previously stored explicitly in one list and then -reshaped to the required `(MAX_ELEMENT, MAX_ELEMENT)` tensor. For the old -version, see older commits (e.g. https://github.com/dftd3/tad-dftd3/blob/ecc50f19adb8aa8baa38a188d04228c4f26975d6/src/tad_dftd3/data/radii.py) -""" -import os.path as op -from typing import Optional - -import torch -from tad_mctc._version import __tversion__ -from tad_mctc.data.radii import COV_D3 - -from ..typing import Any, Tensor - -__all__ = ["COV_D3", "VDW_D3"] - - -def _load_vdw_rad_d3( - dtype: torch.dtype = torch.double, device: Optional[torch.device] = None -) -> Tensor: - """ - Load reference VDW radii from file. - - Parameters - ---------- - dtype : torch.dtype, optional - Floating point precision for tensor. Defaults to `torch.double`. - device : Optional[torch.device], optional - Device of tensor. Defaults to None. - - Returns - ------- - Tensor - VDW radii. - """ - kwargs: dict[str, Any] = {"map_location": device} - if __tversion__ > (1, 12, 1): # pragma: no cover - kwargs["weights_only"] = True - - path = op.join(op.dirname(__file__), "vdw-d3.pt") - return torch.load(path, **kwargs).type(dtype=dtype) - - -VDW_D3 = _load_vdw_rad_d3() diff --git a/src/tad_dftd3/data/vdw-d3.pt b/src/tad_dftd3/data/vdw-d3.pt deleted file mode 100644 index 32626cd35548ba7fbb3cd6c1c49c16cb7a6c151f..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 87275 zcmeFadGrs}{|A1HvXsaYp-iYyBuOdi3PnPn=ibqm+x#uhgY?{oM8IAnOALL)nO?$x75 z*UtTJ>eaJv(aUb_()HG^U7Fm~xnI#fg{xoQr%1m(MKAAjM6FgW3Kc%`yHlq=dZ4<0 z{NkIg#*?8-{1fB!2k8Y|N9K{zI4k`H$hVUhhyJr1xUD~Ez@qdo{VI_|OCX@fFazi|-Uo~-_G5W*Czh{$fbAJjmhJrbaN@EjSbxL0s4tAzh4a&+ zK1X}9)=P0+eA4&8Y{@u2$F&{G{?pFruCGx~Xjua2N>1VLbJ_v@-!+Imjs&_TXYoFR z2IBna%*LoEDcYX(wcUsF!m|Z|QOg!+UsQVDvAEw!C-c0K#B6zEKF7}TY{&H9a9vz! zEimc+66>oxiO17M?ANA6*`KY7qCS6r6a4OPu1#E51(;rZ3^18J0GO5ek=SJ-&<{Qz z7>0k(dY@|xOujC|^A>%8a`<2jkI$-%`zDjeqP>3pXxu+LZwcFRb`L&B%l^1N?fnX_ ziz-yd^Sefl=N)l3>%V>|>mRk5{d)GF#BF7G-F2q}-N*-UT$p7`Q5L$-d$C-17SB7P zC9ca{K=VZ3s-)eT{nJ+orMJ z&DR2b*Qz`|r~=OOx1NdXUGM9FarN@N->KX1ds1y2FuI{0Fx#*U_l@eFgmDyK^E~7` zs`~)t{L7BSWoHB3?()FUVIRj^w{5_<|9hy%FY8Nuy#p{b=JWXi&H21z>+t-`ssY2d zukpN-djj2!zw`MwJjmy0U7tAVbF??zI}*>Em21K0u2_WSHm?F*ol-bIzUMyHQ}HHX zGUsxZe=ZG-+P=ens5H0yKjX|b9m3}`IdV4?LVL2Yop8Q4x$+D7z3`o1W<8^i`VVy5 zo1!d?zSxQTxcdCwMa5V1ch!GX%HRKK!-2nizA?E;8}~r_(jSJipWdAfj4rOi_wv)5 zs4scyS3JLO%kOmbdozqT*StXmf1uvz!`s>4H!sI^Nr`JJ{&PI7I1}&PuY4B2r#BAe zb2oq8+S&BDf6nuFwZwVhlu_)bY$h;jaWvb}<71AWLvBGmuI(AXWKcz5RH+2(f3*yM ze|;bBlT5px*N>fs_PE=c@p-#){H62n!};N^O?>Wd_Z|OFd!izx@m&6yayUOO@Eq^| z=O`R^GfNShbmQ;InLu|^Euf!MgXP|*1GAY^*^V!+KsmI%4)=-gseyKdQG;+i+V&vs zlb!!1&iA{2tn^QN+;?Mee)wV}?wg#}0+?4(j(9igV_kaIn91qh!;P0w``eXS0 zz^(D8b(-cI$xYI@JH}DcxjDv>pSh0XspmKRu1u~9OuoH}cupgJPol{je>Z-I{z>i~ z0CXLotK{*Y9|rKO&S}*S&1Vu!1~q?@E*tzn6PTVtIY>{rme+MJ&+{*PkNr??B#&>r6qtUvf$wGF zNiL5A^T!Ls{s4KLO>q!qvnb`}w=--?0kqko;(W zqjg%a`9^XRpLH3YGYq^O&!1%bAusWr1&J-zLT=KL+c^Kc-HXp-^Gh_U8lEeB(hc{G zt~dtu#r00Xc=PQN)Eh14dOsS|4&%T-^asulN3SPtzM9wfeH+(jhaHRSv!NyU9y&eA z@w>zJBep*XM;(ddNsTqQk1O>%`XfoEpuXhsM{<=x3Z4uIUGK%Sr?Nwzt`ya|`o1 znq5F#6Ht%a_Y1FIb}=x1U=pwY@*H6Fl&#Z;t;PM*3JZau*dqKMntjgqe*7D#FZ=on zw9Bo(l%~p&Z?KDxbH{Vz#r*L#%h+5R{YpyUV3fM*8;$B z(-Z86vF8EP(R>eS^#tW?&d;bP+wd3aiO1ChMjfB$^{2N2Mz=el=9#$CXShz|-1Uv| zd)$37-$U7TK>yYbVy|MXXYvgwM~2Cg$59qWKej+w>%)A%rTV|!jdn|Zq!-ZsM(Z@q zH@VzA{WrcF+1LFbH*P`?uHQ~y1i8spp2__}g$LQ*Mzw)y+10>kR&BK3RePNGU-K}= zk#D^U^(MPc2S#g$;=JT4kA6&V-ivZpa}uxba6a+cnmn%**TLzi3cSzL>rj6L!7!0^T#UiUEfFYe_T zxGpQRg!en=Vtl9KVU2*%A#>Qi=Z65prkL&Q)D)Nvsm^kZ62NqIC)5*HeF66m3pNpJ zPhvZ4o#!s;&hi)65!dni9^wTkM_sQ6#=9@Z_3^v+qTNvyj5r*J;cu0I*=OYS|I?YV^GDXcpl^=Lhnyz>;!(|nPSbK4h59{x+7Y+QOn z$&t`MT$=4QIgB_ls?x;?FLGa#U*<q`8+i@;&?Xi zE?{_~95B9S3+@vY?+#4r9?Np2Pf@R%Ig|IPdM7Y`wF0lV?_%=7M3#5F$?IlQFX!ss zkMq;-GW_lz>cIDWO$)x)!@fhivL2kj!p*1f{eE*3&d;9uiS@Oe1kB#P66nT6z;yK4 zKyUl`=-el8y<5lm)!lsr>-(?=>dWomDS{nraPp8tLun4S10&JSPh z;5xAN3LFmu%AsGv11r$4>7_UGeK~%QCvOB&9?_3>Hv z0Miqi5<3=Vy%RXUX0L2vI~q{$>K|E*-@`}E_#6YSW_|XZiw5rD@3wvptM1`>eJ4U* z!rS91KYQ!2ytWj}wqJBDwxXR;(Um~gbP+Hq$@M^Jzu;=HTDvMIM^_F9y0tBVY4>tG-f=U|Pp)|y<*>yj zKj(Kx_f5OaM!oLId3?T4rUL!x4*-+p_wsp*Pv`HCay+Js>Z2Te{R6H~3b3B&{ySM; zyHURG5weJ`&wUV{_y`p zZX^e?8wn+ELU)i{9XwA(J6}=*0AC(hF#RldsckzLDIxqSYWj@dL{- zFC{lw>Fu0{{07c1YFBt>8qeEf?K}A>l#?%hU_FgZ z|M?6RaKfiP7IkQ;*;!%9AW_JO@lUx_ajjrQ$>2&^n*#uyG?1en($#Tv82~%|3#-?{l1&>z7Th7QNF1^Pu+q zl9PPg+dAxE`O4+Y<|WNfk|WJqS`Vsz$yY8v<}V<-r|dT_n#}hieVXY7w7-dLou>IF zmz$>xVSdd)^mJ2VESYOmP>Fw?jJe{<>XoJC(=i+!S(L!GSD|gYgXbu z=@(USJdLX1{OFfE(XKF~74%$D&y&#pbXR@sTeDU3iH&I2@_oKTIeP6v?rYyE!{2`? z3XIPv08E$A4&us8#r=}P5spW9ZsT)5vK04`-ZJ@^`qpsWc%1KBcY@s`{df-bZ*5cd zx3wc3UY-5?+kVuWlwHsIi%>ro{W%`x@W_1B7oW2Y7!Fy>>rCGe=GMk}@u>BDF z+|w7Jzg)MS=s(x<5B8(!P2I-js3-j5foY$6iFbU<`@VD$&{ceq_x<=5*5AD~pQl+R z)_dOwo_FhP+%LRNJCM8Zrpo@n;LT(Qx@rRYBm4FkV0_vL)RR8_Brv&h3CpJUPEI+K z#|{0E9w-ardyhpqUeEWQ4cvvY>JI~#<2^}!71W1K10!6_}y=;isRb9roFF0zok1D@;y|07v=C)YxZaT zI>2cCR(_{S+|T}fKLX}E^d_o*y%u{nPkJzV|uDp&v5SYo__6lszLNIjwB@&Kuk4*OB0hiX2Eg2&YIarW^MtwG)$B|0zP7wh^_yK)@{{|g znV(INzku02^&ONRTKcrq^a9%7B(_e|d?UF@*S5m_vx-OZ{WPNf(^ne6@jRp>Fnfu5 zwxkR7l7931=&$tIS?o8H8`pX!`ZZhgGv7<|!%>b;oQ8V*5Xy<$O}l~rmh)9|JNF6M z@`~)w#tr$N-#CW%Dcyke2I{Snzd27O>$h=S9916gC9dg#`p){gJ8)gPu1$(YN5S=;JB7ab0a-4CeW-8c#8U*UY^Znt&A{+@hZ*P7>* zWIrX>R%7{nUKegX8~62+7;eU1f2tp>9{-llyonIH2z zaOrxsa~9VrS%Zl8JNjkxLsIX0_D9*@AUB#<5|fM6{0Fo@id!zi_3<9gt6KL(A5!lx zOiUjklsz^Y*p1I2e+8>wb|=YCFuzXiJLG30e}UZYY4#iGp`}lgUO@YZe4S?VO)fVL zpGJH9oVU@BelX=HJU$U}klRJ-{(ye;Ndvy`hv*;TKBv7r`KC9oZ#4wrA8uev=6$Yj)hhm(hCmhdl z`peP4Xi-_*KRW3l_Ve2{QIA`63F=9=Q16v>GrRd9j+?Y>chno-d;`!8U(M^C$+gXQ z?(wHlUv%g`)_1`@z-ZbT7;mBAG|JtQi-10x&30OQ-DjOp&dlDLJ$f?AahDQ6zk31d zkJnxX%s!Ze{?3-oVE8&GH$&fMLQaoS&S#2JMc|e}L`(Z8z(wc@r>vx*qKf zzt#uFjcEr@mi1%*4!w}r{d-&=3f_cvMeFvWoV>vAZKnGE{2ujWb8dtDXdNZF$>qW9 zLP_6%c4)lE@Be|a|DEd?VbpRYjtiw1kEU~7>2Kpao!4)E43Z!DLnr3fDgQKQem3$K zklj=E8|k5A)2B%<5ZnGHU#DrlF}aD4zMbP>64y%yJj*CQ@)w#zyLnXkXgpW))@wN5 z&29rs3eldbdfYnNoqemXQBV5KG+Pa6zl>Pf`N#X_H;(aDN zX@3Z>4C49q3gCVEaokV%{q%2%@4Xw>Y2Op4_p+Y)+Wa+jQZKUP=rJ6!n++8h1y0`FJhH}+$b7#N;? z5Avh^iq>oYHNV*Wp?M{rmrRawd6NH+q2x_+qj@Z^-{ej5qd0G>x0cuX={}zcgj}?sdSh z>jTzXg!YU0j2lsJ^v7h3(|FNEz;xyYw8z&xAHOF>=%i~!6`UX4 z#=UIcr<>W{Im=mY_AIg56kvA73ba3K-4*E4%Xt5!9)9<0BVahZILn(40VWOJ#QAZz z1-x!lN1!jw?{s?W$9x|Zmht&Zv|+#Ael?yydF2N@fAo7RzMu2CuJg0G@5vrL3D?W- zG*n8_F8MWuvK_eZ(SF3%R~qk{H#DEfPeJ1*eylv@(opk|Fh0oN?O=J5UtKOg7Jn02 z9D)3yo%wajKTUo%@)wZZ)0zE7dT1won&|~n+uvxNmd`gPH`$Ea(4X19+o`AAF^c{A zSZO?`TlPHq&rP}*C_k&{_xZTam+gc1Kx_BwaS^tmNUuWF`OlOw_ z`s%M?Tt+t^g7c!c>u{WQ<##=4_dMUzg8snp=6NiSr$1;E{$P1kC0yseSrwKCP(PktJcRG@9{Pc~i+@9Veyy=_>)?^1kaYdn9}t~gNkQGHLc`it3bb()}_tWkSluK(F_29Mi*)Ris*`O&&f z`}~9DL*rj^p?N2_yPDs%`9^Y+^(=(`aaXm(dLh1OEWV@eofiC_4yNC2 zw3v2}bUfo;;?cBE#7BIG_wIVruH*O7FVVH1%;TRo<9M5273-eZ>?EN+$5mW<7T%Yi zeh4r<_bs$198m%L^k}{5L(i$ebyIRHFkD;-=;krr#0~F`_NT)cj}orDg3q({Qm%u3 zxe(72)vC+$CRV|HvWu=nJ?YV#@q2tB^h}eo;(a!p$LoHs!t-Yy#q<8$kM@NtepX-%wwgHKTue^vl=iU)PC#QOV6W zq8z?8J7TY+P)B~1MeZO+m`R~Na|VKjjK>jnw*03lk9Mu=c;{%`}<*aIZg&d zXqS8L2_Ap*E68QA_?Pf3_fgra10jFW_Pc<|9@?9d?+cD z`~v5*$(%pK>Ua1&)&2&?g?9k6y*1e{TaQ5f>6>?<9`{#k_RESIkejUiF_ry+5iC1u z7`g%HXJ!}5WLG_PI@@{Zb41%eX&)tho-j1Jl=uDdCZ50HHKOWYQWa$g&w>BE@e3$U zTJfHV#jPp+MsWo4hn8PwWd3RLvys0*FuSMhH?irVrB91YFQENRuyvZ|8_A6?dL#NR zT=5zExfk`|Njt`AWM>xPdpP8AjCaLj2D5Xgn`_{BHtT1+zvyL~zij=Vt>S*t-FG$a z6IW;4nASVl%))$sn{VOxxXU&s{+Kc6t}}`HA&Mx_=|zi`l2bripyMtuF!n9rv?87Jb9} zwOqoo=7F-b0|na;NgtQ%wNKv4_FH_a{0fq%`|`M9`r#aug$LV9<*%Un6@Mi85idt< z`~uG6q!sTOEN)HlH<85=$RApMozDEzazC4c`~_scksez5H0cGjztK8P^Nr*t9NrK8 zm%MX6`}dZ=xUXn35$_=~eO0pC{AMoY`aiVecQ5Sj2l>gKr(N5hW$`9Yyo}??&-}ip zhh2sG-0k%!Z@Zqs`%9+P`Wx7lmmug9>w zr$6dR7hHsL7}^{2uI3MoZ`o0E{{r(LkY9oJOR^UUNN&A&i#5|Wp^ezV6)esVq= z=18qFkgy9>Vwz!`pG3A3?vOv@`80SsUt|veHlR zeHEo$DA~P*@B1?9bK?aZFL9yAIey-%#rIm`XPlQdroJPYVS2yQsGm(Q>`FQMr3f(V zNPWE@bT!*O{2|IkiN(O^%ra<4RI@+I{g0Ns}tu>3sj|H)+!vLFAZelkkG#P=ig+{X7ZdIp{=Ub&rU{y)C+Ybb}~ z8JC{!`<~DL?Ddo%TW96-iu`o`2k${HMnC_8`fRc%U-lv&oE? zBN4xV#YsDh_srwg6n|rJ1oDRt=GQ6zH2K+h^B0iaQ}!F_p>utj=>@dE(K^lBd?UGW z*B3zl`)Tbt552t?daclf@%Gs(&$GOB2K8O$cOIR+0pmC=xsv13;!33-aZMh^@4mvX zs5g6y>y@Ny2i|YC#cS>FOSv)sg=lUAzTef?P)`5Yhu?jNwZLfSYdrszg*@KhQLZ}D ze&lBl#r5*HjA#9YdJ=)~g*5}~;-nSt=`C(e@i&ph5y&4pF~3gf`{ZXMe*xJ&WxvVw(56q5UO@XBZ|gM8H*>euLvt@r(l>r2O@v-8Q@Le9HOUzah{5=EpcLDp22?Y~{LA z_Lb>#S@Z^;RJ~NclRh8!$ZC1Q>U&OzhA7T=LV+ zd@j?cIrBRU*Koh(OI?F`Bs6RSOn>8i8-ITUu1oJ}kK<8U^DljU59_Np9_UB#Jhy{> z_|cZ*&_Bs-=dvGWUyA36%ZJ&NC-A8BvWzEysD@tRKj-Sh#PceGwR*uOyb)ZCxV z{N8eVs`>fIJ}Nxe-Xs0E>X$xS^G{^FgW^vU50v<967h0~U%=v|oyB`9ZY__$u{eTY z{?PL4lz*E1Z1R3HyJu?l8|k5?PYb3O(Edj2G;i~b5Hvs+01;A+YqrmXuZyb;2f9>YpiswqU zm&EVxG3vuYk&QfFWf?HuF&Wpn`?!vaYWK(Stn!)s-TY&s0b_aI)W+x+*LD@2KkEK7 zzaPJ3=>M$R1U#pIfc7aj@?DnAFE^2%-~5AeJ(;XhZet?=>@dE(K;=% z`6ig$M5A`%z39Ha1NC8H?=HMY=_%sgzoQ%`(Z9ykA`hDX;}2fnoV+L6Z@FXd$GDOn z%zw`HM)XsC)DwSG6qsB*7~>?p*K@p#H~AU$J>_O8<40w`Ox|otxp|NA*U82>GLv1ms&;bzEBdc#__Yu<^FpXil8Q4UXk3>42(Tw)x`e%NjNeJ=eF z;@9Xe7q&Fx_%}UrvXK4_QPagZ9xa^#3_a=Rmkj3oo1S9+4tM+vdGYIC2F5F?A4pE8 ze`7R+-;=PKakg5I_;~3;vLNS)5HU1d^VBsa)@6*anjD>J>{>FTKtXd z$nuAlU#I-jg4%|j;I8Ai0=H1iiVvM`AR3(g~Wc4FBp3NwUda_c~s|E4= zR3aa4a)HJByiWOvFB${vz`^3f{<%<{MxH0>`KjwJK%lp&^hTbQl95uL}a`A3u z$WOHH5tO5f)T2e?Uq?B4{UNruNI#Y>p24j-o5w5A?^yFjGU6~kSNa^~<<$p)>2}IV zF#UQ|?HoRL-3s_UuJ$3y;gT^xe+>66ae?hM zHqc*5Jx^S|JImM8-#g#e{mcA!#1HL_=T?5B;vE!!ns}gcJ{#lZ$oT~|F>b=v+c=KO@4>k6zrDkG=S%XJWS3i< zbjdd`Ucx)nYa}Cwv%Gg4O{O!Q((6;=&?U}&nrV7N%y0YGy^pgwYGnA7#^!JD-kjF0`)sAvh zrzy}KeJsvPxBUckr%VIJH+BGKeJksVzPs_qQ>B>3kkB|4{k$#em3$Kklj=E8|k6N;~>3&_BUFmX}-zj=ExqDlVSAB%8K4kc`8YLTRz{u(Hi44 ze10Z=cgHT``1txKjw8iY{>1NaJe%`Q{5$>LqWj4goSgeU&sF`g^;s~cYi^=DH`1Et>!$!WZ=`2#4vT6*a` zPR`;Gv>(d%|KjZt@1kL@XRtV*gXKkf1Jy6Nl>9`Nk3hWY;)f36xi!8{@eYbVO*~NI zvk@I`6a(`&^>vZOyCO;eb3&`&2%zh(1wDf7x3*`G7Tc;&9 z-*}T7_vv>SS4l<2L%7mk@;&`T{X_E1ZFrw?!&1bv=ns&{p0REh##6j9;kc@q;d<@E z++ON`lTnnztbYxTkB^H1!}t3sw}*25mYi1r^GH^*6WjGQ$FDnr_SLio?OLJLT-+yZ z(Uta(2T$YtWA^^2=7p4(XO0EBb8P>A3;kr`8|dd2E}~sO`-A>x?scxeq8&DmXx=|; z2jn5=NBv_3%28)tmo(Xn`$T2VV!O@$nJqjG{gPHraeg{Qq2|S);F@cg*<>%r2ALa}-Z%JUz14X#7hbCzN~% zB~Ow+&2M@A=7%l$36{qrvwQ@Zc-4&`T0FPj_&RewR^v|-50v<9#LFRmfkbiAiubg* zHO1d3jvz6A=pesN^G}nXz5E4a_ss1#riYe3O?rXY_BUFm<@1fnO;YA`yf5j+!#Q_Q zj;?&2Mt{xvI&nP}x37hI-L5&P*ZtfT<1KBum~waPjlksMK>hS3kFvcV z-^ubx*8yF%M!@9q8a&=;CEI&@4D|6^Xos7!6qqcd--Vk$mF3e%12bFCx^j#M${M^% zc`@FSuxAqL3q1<6e(R6)u!#L@d@ydv&FtT}GU|G@`p~% zuT%bMnfck|{sLzA%N#c!;h7G{`%|4~xJB%&B_xslz zjE|^9Cyr}-ukkb77y0evm6zW@_6O~dzAX8H#FOjmQNI>i9gX{hWj!fZeV1c=#BaX9 zaoC9Q%dYu)j=wHNi0|~|_?vkL##ear3t%#zya~SkK3t!4_=NRzd!D#v9PxANr_=dY zqwE*mg?fGCV!-T*nLu~;79PKR888}caiHa=qnuV>$lsrwm-(x_q@dU zGWnYGPqLW$dht2BZy3KF*Ivtd=l_QKLKE7ZqU-{e?_)dl-ojhOG5?AeNBbc8-{kYm zj9+klo^L^ZIhvoeo)b@kq2?j2+oT^?{hFUN|HR50W%+MX%j2Q=6!EId&pH^-ZO+SS zyo2IT6AzU5Y*ORp$oT~Il zaua{v2k$k0#p2bBuR%OK&iH4tQ>x+p`X@_qe7`W4`-|goIpZa=Q|X@>?IGV$T9x1VY#;Y^Y3uGd-<`nuDyx18%OBHkJpQW; z%JC@rt7jGaP=2D*h?aLmJS=|xTl{VquXl9BVmwd!d`moE5_SS(lVilE*x7k`@ZK&c^hmp7-s z*&V~MPKqzMitDC3`*Phh=@^#(9F6r+WITq+O)XGPC!9kpbO!OemY65vk=)-Wcl4uP z$>M2()>ns8Kao_Z3VDdzU0B)w^YZ4NLHsPib=mkEfqvB`Z09xgfKihZfa#bcfl0Z& zXs5s7Q{3O(QUvXZ*D-D^`E3e*kM<7&y1#o7PascQsA+KxHM*f+-2Ejup5Ix<`<-q(RGvcRTl1DTO8IY;$3yuDQsY&R#1Cyex8m!J&Hf<%H1R-*&qll);ulby zwBkLT#jWM>Hx@@Ae`xu2%0EqhHu4vc-7_}(jr7pcr%5l6?{6gUHs45Y;&R{P{YJYR z0sYY5F&@HQ^m}oB8@|u!S6siiajW<{t$D(=Db57I+uJ`-wvm`y$o*DKygao>t_mHg!Sq%AM5@_Q;zVX%B_${Xb^ z|Bdo^C?7#&yz1hIj*RD4e4UCb5f6-bpd#b5iHw&+`~rC#i^X%r7PqGO8*gz0@`p~% zuT%bM^0Se@fcRt+v)@P$Ej_990@~kbofd4qk=#U$x8r>#e{sH%e`{iVZP5(ITSlK% zh1~elDgdJgX(x5ZrCisTAGe>&@sf;O2Kh)@y#&nfa~Sz}lLFk=xakuyuCnM}yk9?U zA;)RKaU6f`7{?txzmCUGXB>An^c-H-;#6SrHvQ*ZF~;Bd9p~b?{Kj`UUZ+3B^4d|r zXu<%Z#T)6nlDtxi^>+T6zrSqrj5q#;JAUVPWj)us{wRy1c=t8Pl|R1&Bm3T@Nr>FYG zuPyme9(LuER$kmZzo+FXj4j_lE*x_|wD#B|e+P zcsax`kXoE{9`9*!Yq7=O-#8?4GjUNDr;`p!5Ry{>Ij6ns0Ks zSxLcpUR>v^fJucal%ET(EVofO+KrH~3;vLNS(=6{^WPCQ>csXL@7jPCQohshb>`97;P#l5$p>w}Z^G}nX zjg!BC**#^yksez5H0j||+uvxNrujy46FqS%-fL8AGRAjy&eQnbW`ob=dZi8b>24DJ zDU+)D+FtZw&Qq*;SwLx}VxoK3ado@jsaH2H7cZ z;rvwgp60wh>B4yT_}+5(-NodCi~GFHalK^-P;m-zg)waR)3@?Drgs8H&s<9Ul>Qp& zO6ITd=2w?AS&sUB)#HHi5yNrj58h zMATon=lK-#F=^#gA%w0qt+JPSboNxzV~Y ze9Cd>OwXR~eTU=p>KPoT%Bx`U6Y;Ex7$1uNkF#4TM{WP+xOuoZf8Ws=*Js13QGVLB zg#1LaExzDU@(v|!7>5&WDui*FuB1IQ&zoTW0a;zMQ_0`E(?Fb`=L;~up=1T^PjNTO zf0#c4?U4RA9C|XY%PwdRbh92sd!q*n0^N1JfbpIJKz}{so8qSB@cl^JeU0Pd9f~Kf z!t=Z3*RsCJJ$b#^m9!7hdP?(+#=HFCG|qGWf`j6_57w(mzn1I4EkDD-@+5n#A%20(;-r-qDN)>-#ox%^H8p=|`E|-aEj9a|`~_t93}(O4Je!(6O?mC()jh z=P?e~avYr60_dkxzwTDk-!u7;^K4{(+#zm_`-wL`s{aVCi?_YOdMY*ox>pz<61_4S z&!0>^57(tt$T#LDehRrsZ~qI&-6@nq#Yt)3ulwC1-iqz2%cuUR>q(RGvcRTXU8-O8IZR@=IGj0%yGH;)fQ`ZO+$e zyo2IT6AzU5Y@G3OIO7+HElyhTo{C$GE&fJv1oDTLUuW*0W_~t`Ba+=y_8aM;rB9Py zK;H$e)AISof`$uH{d-PVr#DKf$u`sJau2qQ3M!+Id~; z&v0G%f%`Yb-TT58N73Rlo^?+&@izHs&uNyQl0oiRq!GPm^B2+x|xDw8-Wg$xYVt zS&SF$rxZ^ZP1=lpmH$vywjRcdyZk-ko%1-Z&Y~Z`^bb*|!Q7|+{Rh6&(U-LQD_>mJ zn0_?Tm+Nr7zyAebRDCq-tI-P>?^?s}|1jpO$ObOKePhehlf6v;h;*It3yH7S@{goP za30Hwlm8>k<31{x>v+HK=_i}kqun!YIuHGjoIxJRtWIId&qZ&uJ?1|bx^g|3n%`f% zgza@TsfWq_W%A?yt^xURo3CZw-#In#JZ?l0%s+nETFyg9m*e-+@~|hCPueLjuI2Ys zoC{+%aQX7Se$fb@t%rXQ~Zs$ zI09WSzs}r0&HQX~e*v?5%6^mUp-rDAy+CaH8?Dnc-$-skd-9QFJom`A><9>4u|2ada!nBT~aV*Emw{s`tDzmWR9Z1GzhXO_=i zyid^xkK>v*qxySMUl!B9D0D8(d8iHJ1kx91pNjW!yhm4ay_r<&jO)S=Kj8fI`thhY zOB({UUQKI0j`PwEX8@!0WnO3BkNCwsyuQ|}e15AZdg2<$P1=Wc=x_r4iIR4)F^Fi<4HoXJT<{ioa1D zK|Zd{uQT^gGd~;o3&>vU&3=>Xp-rEb>ji9oqjg$p^Gz_h3FmQt6z}|j<9?H^W8Qxg z<@Be?z}WZ|{rY=Qb{*!@u2O$D*KzL-LfJQ3Lw(#OPq6$Rc|nukFGtx`KMk0)C`DX# z73Hth@zi(i;l5q@=He@_=JD(3XCpt3oF}2oL|#9!J}_MLGseI481Y=LN2Bk1QC_Om z1csjzj3en$qAP#Jb7hzO#^*KOJMD|Isz;-KcLeuI;&)GtpIUnBoKNI`@s27Vs`u{>Y6C8dg#!WVEC%?;0$Ojb$yb8H-wY%cFs5|YZ>EDc32&*Tdo_I3zA&UPk zD`9yfCKkl+>F$A$k7#dwp!nUQ4-=G=CY^x(FZw;i)5$v$z8-?(*-YkZNN&3W^~rA_ zD?xDU_GvmI1@jwD(4jF-|c%n|9RDk8$SFL!Bu<#ygQ&UIne= zWZz81cdfiY^2;`!YuR^ly|(2UmV5|v{kO@R>`8h3#&0V5iOy^cd2qe2<9yPV{)}_@RyGR(ze}9nAUDj0Z}5HnJZD;}=kzwBkJ#x2E_T#SzFKT7I4KPm`Zb zWc~uOd&+*3nI2mDH0cH8SEzMbu=&PIZhoB!xe#xDs(jZUHlZB-@fyeB8|xuA^0W2h z&F^${ZOY*^`kjTI^gD>Q>_WY1<@LZUyg~UJwFmVo4?;S5GJhX+1u)seyvCslBsZ>&7WKH zlN4>w`Lud@jQ{wa8hl@WUkUk9-aJ2U1)kTfdmHlO23$?Kdub5npJdLDxZXdq6xT<} z!>)YNmKRs~J(Z_0v3zUF8>Rd=vZpB@fq2zr?-b9i_&UWqDE>6%ffApMcsUaB3s{_V zYVn?_#jQDuztMV7`2}*nPV-NbpN;$lWcQT)MtW%J(<0Lg#J0bQZJnn1MsgFivvIkN z>ndO4LcWh*W?=leySH(ib$k&RS>BPT%`A+&G-cd&*w>ola!74pw1N2+qA|<^BOdEe zYc%SQN0~jTabL(mI^arRG@=pEUBr0zWLHs?L-EJ>`+4U9lQ$T5?u&6>mekva_R4-6 zmm+_6s6@YvJWg^W=Q;PP#hI^Xz9RS7DAsTMbLkP~(awLH{6uZ{aonG}4f~%kY7pgQ z;W>~WxBgP@+uoQ%xv}`CWZavOAAkEpT*uX3$@yngHT0YI)5^oHeA2Sti%(Q}3YBk7 zd83s7CbK*q&hin6S6%$j;<*)Hr+5eR{kQQziO)v79I^2WC{9}Op3dUd6n~>Q0{KJB zuT%bM^0SG}UqE(G*>9Zbp`}k#yYl^wtnu$69gRbR5L-csKVq*~J?% z-z3L90E`d27Z`m`f9dq6_9zF-7nII8A7%NqM|CcT{JFcE1Cs|VA5nMitFoiFP)>)x zPdt_LmLKyv<=OJ|#1GQHIzE~CBa%yMTG{e~adkA2!_GmFt<)?!)>1 z%LaJf@FDr?QsuKBV17Y|{ek-ab^D)Z3u4?S)y6@7qHPaSUZnT@F{9i(T%PN;4|_mv zqJeuLH_DglH_#57d`x{zeDW83-YMjvR33I``J|N>H?jPl%2TL(Ys#~({5Og6cvwCH z@v4g-Ix?QyoUhY(2V>(;6AzU5Y^3*#ZN{D5z@ zi1+XPCarQZ)W4noCFCUQ$Nf*b;sMIf%QGl9<}VipUJm(j^<&OI zn>s=c!*#S1WW^Svzp~Z4fawL)A1M#J@<}T%uJU{4c?vDxn*7R?|3-N{oaG}Bue$i5 zbDmrAbsFzrF#a_0K#9*LHC~RKU%=v|74NCIHO1d3j=-BgwEQ{~`KQUx#{30j_mus{ zn;u&FwAAzhneA^ftDX}{<~K=PPKlsc~E=^VR{{&|LSqTcq;R7XI1Cox~$?XU|gBP2W2k<07d`-Uqk+8(wGi$S>ZnqhCVW_zC>3e2R+0QQSdp z4^m!k%Zv42ex$Y^I@qsG`yA;Jq|ZCpA6xq(?SG=uv%i8Gvr;~#Xu?$Z%SO}AV7+&LgnIKl?3Pbjd2#dn zp0@6YE#F$6H_GzgD36EyO~tD&9>JXF*7!QbJ1G7%@j!X;*%<$GX8ZzqoV3Mz#)@0B z_#4F$r1FP0zfSq5$3OZJn0SHzqf3Ohv?xMd{JN>|^Hn zOov~FeooJ199{B$d6bjG4s0m!22p6ZlJ)l;1oSVQ2~6j2!+FsuOL_e#*8-DoKF59IM=VbF$U}Ji zB*y(k&s9g+zj++aa~Iu!??Cd|V#u-ov+_y*+&5Aj>JcrO4|Tz6_RC9mGrn%*1H|hl zv!CkDA~tXHpZI$(mW7sAHLG+zj=Lk~^LMM??PUKSXafEsUYN zs(%yv`J9U}?)}g7195%LZ==;F%umU-g^;K8<2g9r_ua?qj=vum3QptklDD8AqKPwb zpUm>GE1$IT;wrzV@)RoHn({`amj6b1JW|U?AYS#@_@Tvf>y58781JC?)5HTMKAY5d zIb!1%$m66f-cxaFioeM$jvzCCX!&(!=ARbK&qn?NvU_G`zmXnV`n1^e0@~kbotDox zlAALh#(RyYcdYFHd6PxQ0pkMwfZ@bzfaw#A_YVGal(V@f0FyP9FwR4RFM#ngPXUvU zDG%vy3;6r-PXm3Y$B7Sf-{mWOk8*b3$2@P@B;qfffzgSl^Y?G2<2{EbpT>L7_KxRw z`Ngpu7vmUz9&KivY%-sEI3Jh8^{#Il$djAEJgon={K3)4lR1xE^E~7r&mY?6Jv^uT z(&p)5oTuaE4=)6|B^RN;k~ycKAJbxtJIr2~ zj{7JNyYfjZFRt=?Do>&Et%)yRc?P}mcvwCH@u~;oht7F!m7i0*gT|jG9w_nIh?m0| zzkuST74PXRZY{C+8^saGA6kB$@=w!#Oa21E?4GjUNDnQ2n)Cws{>Ij6nr|dG>CXc> z&Kg`vEZqR;%TxtshaC-!ckSaiT*rAR{N9u0aif6o%NC#Ak$JV<(uZ+8*))>py?PJO zk7~j5N7Q2dwaNg)8*6bq>-HGWpU1q$@h^AbeY^9?L+76#iFKI#3ghOCANLV?4}6an z7$@;U<~>gv+Wd4AdAXd;C(64mdr%%vr+8H3>B;@I6_5SwduX3{u*7F76hEWZhl=aY z>o>i*v%hS4W31bvy4>f;kJp=BI_DQWWDcG?{@_roTjRZTh+`7SkNdVWx}1Ce4XMQOpQNHJW%4Z5if@`eu3EHr1N-Bi(88<{zmq9^`rbc<)0=$ z8~F>!?wOkXMtbO6pJsXi?QgVB(|qGhZjuuQP);^n3rsdP14d)(aeRMSg~v}i1{i;M z7RHhPp7G%Re9A}EWi{ns$r6sUqAvr}x?I`m!3m z`MtclH{bL5mth>GT_<3C$X*z#**sI|Qyfo=T3mVy>g%&|L0CtTwBl z8~37L6-Vu#r@y_MUL4<_bU(j`iRF`4UR?1a~XR2L&wH* zE56Recn8Iw78wte_-r!cojNcO>Ap5ghGY(A@17wj4-lhRY4LOoHp*Kxo2(8)aC_9f!SQ9PaInSB1xJTr6y>e2j? z=e4(aM1Bmy%;Lc_%X24x@4Wu99Z;W5KQ2uF4g3CV+1Zc-`IW?7S5SV+w1?b;P2@>V zb}qm^N%^_s`1AHCr-%QACB72 zetMg{7tZpqE1$HpytvBmnOdGgjdw6I{Z-f>;J%e&-arb72-H3OM6jVj_a-HuP1q4uL;D(_W_e_ck%oeUZ(s! zSqSo@bzwYl8Go-wIg7u(kk<`yc+c{iOF#M!-(C6T`Zvjk>iUx}#4Q?v_wT3k{@GIU zX~i`Ma9n8rFCMD)_n;l(&5<1ZPdSo&2qh=NIk%zRxtu9ZTx7iOX$>(iL9*E5Jcy3x9%F7() zb53e89^d`A5dG@<(f%AA-2l%K&L58Vl~^8j<&$=n7gzpW%2O!+jMVZ*DgTXlw1V;x z7_YkKkJxx_z43KM#ycqfG-o_evGLi6mqYvlxt-nOJr%d6_#4F$$RApMooc`QY=Zd< z$nKfjZ={DdeVX(F+TUoMmd`gPH`(Wu8`pU|`ZHd?-tc>V*KPk1FKEc)Uvm5>y$7M3 zzSYmAJ2V%I+y!>$zyct^E}VwEXwMj9KKnH_pLk{p{@Ce8^4IwDRGnIalW(f ztgkhk<(H`^j9)ds4C%+j(`NbXHBV^X_+N6QxbGZo9+JGsFIDr9_Cb05mj6!jlXQOx z@+1H1xcOAB%hnFVb8BCc4PhK``0Q@XKT*+@)YtW6Tu`)`_SC4zbC9=Wz%Tgjxae?T z_9*$z;!jJV-n7}XoPYK+K3_be$)Q{qhsHI~e_5fbI~Jd}?hHC}b`L+3oV#@8v{!N~a2!~^Ax&qll);ulbywBkJ#x2E_T#SzFK zT7I4KPm`aG`~@P}JydCG=2Yr?QgVBb2i^dZZfw9{ptI^&3?UQ1^Uz7brH)& zP5~zGb3GT@^+Gv3*%=sr)ee~c)Sd6)u}v&rL%o6Y-Py_;a9s1V`=0h9cl65`7g3`} zF}{5IJ>GvP&b$_B-J@9F)bsJ(4I^!x{M7ZhUwS8bzFf|iShcc0Fo-kbb^EXNnaw8$ z%aP#V75P5DKLqE8pR2Gu_!eGo`Z(!@vW|y9 zZ=c&yjn6li6WgEU@}c#Zt=BaEbF}qX&JSn&apH~ByrlZYyCC_A(?;++i2J_>Jwh^b zI^;q5ThrbF<63df$#%2%l&lRo%06ctST^<%Tpu1`+@j_gK7VY@>?BY1cF!EAfqLzw^K3MskqLg{{Lh-{ffdRdadLJhtIA z$g_CnBtL1>qKw;K))(VH`i*(0v!$Im-YvgvT=fOWOZdGB#%*Z8d;+P(FGT(6{}!LK zjdEpq(fuprRf*S9ZyEV#(VyapNUGn5vfsmen7-gslpo{sQ=CC~{#*33n>!mAU+@y! zaoRB`XLdj3lU82b$ntx}mZwnp)>7q-vivvlTTwoO$avMo51o~zokDz_iSZ7`#-ApC zZ}HiPmm@KL0mVrx-c$Ag#os87K>pD3>kQ_frsJvk3&`&2%zh(1bY}WAC%u5}Z?sN} zZN8D*M1K{=xX7+N9^*B0)N8xT9^myqlDF6OuZ?ol z%jBzW2fW{;+|?Kt@f603yB^yqKResu{U!V9ALJA#b`RwzG(R2p%l&KrW%Ezgig}{_ z>m49J+41Iw)SmH9@rurPzOaBi`B4GJ&t;9OV*IAvX5+f_oLSsQ9kLev5-+lK@9l;7 zJ^5e~>Tz4DLM~m(ydCb9j+Cc|&tv~?I1lejc0@npB=+CdFDWNWZU)B7PhdT>K1Tg+ z?yW%O$C18Ad2u7l?Bkt>sKg}Btl=y7K z%OQS&#Nwo##d|7lP4PE*9D)3y&95^y|1|m8$X_5|XUTqJdT8m>q!-9+f0NldEw%YZ za+40}#_xT@1?>NAkYIX`(NugtwV)}F}@X# z5RcyuxrzVogz@g4oQLbe;!%(z>FL6|j9bk1_vCXgn}L2vw%(0#pC z;WkyIY6x247*aJ0Le2Nh-rvu@^|P$I?z-#V<&Wep-S@oj@Atgtd7eEz``P=ubJY<< zKD7E=k)KAMO{RPS?mfA`Ar4&-Pg7ig^9_C)_NI{B6u)wO?CdG8i_O;VV)mpuXlb`I z7vE$1(JUHm&*>(Q7o;2gd(HcjgMYX@mR-G%3i>I#-TRx=-&7tfIE9B_X!~7`yWII= zfcK9e-kQz#cvG6Oz~yJ)Urm0vmudCS@W1fO!mn$=yGBl+_-Ejnw#Mne%?kNBWMhvX zwBOO=&ojPZ@{;cKzDW5H)pHY9c*gs(RUN!fRQ32&m%AN~HTmJby_&wB^K;t8_&HzK z>$9?v`|Asyu(T_m5>fl4?!nvg`A{Xu;_-Wuk zfoB6Q2lxWiNmK7h-5T{b)De&mO`MDTH1cd3{ZhF20ty6y+Ll8 z%{sc=eDop5=GokQrSUcQ zENSP_&fllsZ1$&-zm@kp&F#sHZ+gkY>^)l}7mbd5jcZz>F z-dOtZQ~O?VpXUpflTNYyDqkFKwCL-3_@!{j>aDi3&r@%oQLnh*D3_np_cB_m4?B1c z^dY3bCw&U(w?^M6`rpvUBc~sMaMc^(LxXb*z5%#{;HTBXfdbDa36}$W0qUfw_bgPm z7WFr(BdC-QT`RAX{IptmHjVNHD&2c>e?uIacp7m5&NujJVQ(ZiRr&UaTiXSR}VBR8|U$y?ZHt%etPNh;3}#{F!$4oOI60-+bHtZoT^(_H*hxv!3Hz-%a&=@O*pS zyY?^qJ>T8XJ`dh_cF#1Imp%s=tu`5GH0d0cs}Fmre$s{d;?m!fK85sKOX?dH{cqIA zgMI|ys)G*=&TRo+r*H?sPXh-EJR5L1z!#uSx=_6*b!*h$gdeGV=%l<(^3!VN*^n;~ z_nx}HAr4JEjkrLi^9_DlEql{QZc=r-$piTQ#amk#e^)nu$7p@A*R2wdD|;?$pKtcN z$M(DG`h@FW)l1S3v!ajKZSv=yb9t&xTFu5ywcl5icX=LN)DO`2M|~KGr-#2qegXHR z#L+k(i6S4+kUO29h+D__^8b*Z^6De~{PXhd#{c=>R+9Ye>+_mM zI?najqon?x^eLp@8hxY4|E7;erhWwAs)G+*3Fj6&nS?tCep z>ei^g$*ChyJ~Vlqx$@JJ@@xw71$6Jp{S9$w;%USM8l7+O)37(lO?mY0wlCGVF)klJ z`OKfsyVv#HjA2I8x4bV&TKjcd_vYojjHWln+CDcU9<%+a&b`$372M41GRIR${X5vt z)3M*Q@mQDW0|$Te^XFLpnPcZI&^ zcvR7!ME!uMSEcS0e=Wur9#+WD0q2|jNh22-&HMe!`TU1}w)ui9bi~sX7vOw@pVr9UAUE~9 z9~ggUD;(qeKVem)^#Sg$mz#V4mZqKKhtwzC;rsl%%X~kcyvD!x*V4y%X_Dl%NI#xs`xSWy%l$onr+U}(M6-6CZQb&&o(Gq2SGb%$>G=$G zmw_ff#cm$&DGqte_v8IOM#~rP@^#(gNE@fR!*R^GFU^mB(|%s>JjrO_zns|7p22CO z&Y~c%-#O`3%il@6Y+?Lgoag7g=BXdJy_tQo%g_DeO>XK9o^^li{SEB7Z14A{ z#kA@6{iJ@4^cAf;Ij&06ca*JHb@n|rZZ+-^=jYiyd|h>(iG5r!!R&fXz32nJFX@GC zZGDS_ANB8TevaFb2OQrhf5QDk`olFp>*f2@)9WLeSsS??Qy(ers|xjJDzc6Iz0`+2 zqo1_;;%4gaS*uSW{nqFkMgNyndN} ze$C5nck1pvT@O$0X6xMaKifW^P29+6z5hRaKexWb<>i*$?C1Fgy?x)*ADFlq{ua0% z^!W)KUhsN_8%ke#X!Hl_<>!sSiJ$hCzlZLp7#}|!|E;=lPnVwy=Nrw(ZDsSS{^B_L z%{%RFzS&y;@;Jfr!;DskE^G3WoqCMhpQnyB`KjN^?eEgBw|e?z8!z~(oA-Cfy5Hw= zx!&K5);CYH@zPw+x6TXq3!+c+qnp?|(m%UedyvcbDIPz_jy%i$PVu1QWfko`UYhp4 z&gFG|kB3yfD}SFuJ2;(jz5ibQw2G5ovgh(J?{h^T_ImjqejZa_T>5)PpF;Ipqi+=b zZ!-FL(2qd4>fl2YcL%Q%+`+(46Al!3Hi63_e1W7oY3e=eUL)LotNw;M0`j5B>m)yo zJR9-_YTbLr{f*+##M6iiaK6D$!`>8$@&hH!i-uAcJ&vEeUdAHg=SI1vz>r7sJ zHhrpnKAo|%(dG`ncPKvl#Gdmz?zjD-FE;0WaCyNE1)nLHZ$h3fbGa-Tq#O0c&D7sB`V^|)8uIle(Ca*IoKaD(_$QR(=Q};K- zp(CEAxB%xH{50%MjocjH(f4VJ@o9PRcd8?b`&vHCF;1DnzjTS$;(){R5eaEU}pLc%#XYSws%;V?j!kuhivS)kxID5|Z z_dH>W(KO*?r$cr%TA%j5&qs9+)zX=Lzr1ZjqxA*eUnie=uaCFX$M(Il*DqY2-+0{i zvr-@SLj9xz4^{m=EA=VN)NhTxQS`r|4^X9k1mLQJ4-L*Oc%9%5f}d7BzOTtMcsAg2 zRKgdaPMY{Qb!*h$aNf(vhgM!E`Dx_YkT1Z!C-*nRp=-s{hzoGO!B4~9B;|AUzryR6@-esB_o<_; zTkwg~b<=&^)&2fE|75nUOWtWUU!P@;w(rwVn!dR7_oPoD{nqFkMgJT6c+igkTy^lF z!MP2*PT>xMp9T&TcsAs9fGb^Up^zTm>r&$&@tu3XCVw(D1}vvqFLc22)y{pZT7dYjxd z!)6=JcDmE`Z_h(~p6hzQmZtr0jsL3g%Ns3T9qIPxox_agC!FJQ)ML8ObED42XZ3%) z@6Qtk87I$El%@w6FxLJx8QYxJJ<+64IHSz zvk@)__yU~MsQ099jrtqv2*9O_yiVn(70Rf4RTXn zf1s^*zWGBo-}JG^gUYQRwCAkwI1~4d<#jjs=e}b7>&q!;nLMOTRyH}vZuj$cx$Z7D zzvf!SgFkhCtgn9B^hW)h_oXXl?rGmI2P%)@k;B~2IaTA#Ec||K^sDaQb{gRJX&cw$ z_C)Hh!)9eCD4QenS1hruJN)p*{rj9KSJp?n7JGrq8)X^ZSSTI{&nn zQS$ldM_Qlc{Wz+kxBZ=BWkK+^e!yb$aSa#PIo z{zc7R|Md9xDvmQ+kVmij$b7|*ZJ(%H;2t>Q2)h3ShnKhleMgAX;qO7g>4iqWd*S?k z;ht4}Pw?OJnP-?h6vv(Hc4XHFY#!)`=1I?+%`bBQr8@p`SJ> zGkO=3mvY8*=U?II*I)B~PWev5O%C#_onP}Vo^Mxf;`TZ9Dg3;y^{fAWzI~p1o8rL* zwvO2`^NptOb#nec?Na-D#adoZRQ`Qqdv0d0VwC)#eDCW`Zn7(Hw0TxD-?j0xHSe@> z%KB~p-Y_4h{?l}OpYnpQ8%=Eq0E_+c>Y?Ab{9W^;QTnaWhn;@XoHvLs(5H}oYxIqx z{|$XS=tlsqI{45fIJd&<1a}bpG;pB6vjLX_e1S}L(z)tAsavD|rc@n4qkL%cIt%5e zk!RD8FQ9u*@+62u6Hg;9!1)G04SN%EbH3+WSJNCnCV%J(TUYLZs2l#O$xpNOT`mVb zMw$FHOP^%?RIPZb*`MO)xyv{A@%uU>m?a<5SjXqDRZ_DO= zY~zAs+g$uDAMc#*K7Rk9ViH_^vK{Rp@}1s^(aZiNq! zz#SBR8aPnk*;L?i2w#9YY3e-#_e=FRNp%F|Ll?^HBtMNjn@srv+m%Q`+AZGcSGUpD%>67*wpusw?;lrqh5?spl*PKEvX?xF5~UAwEvPgOl*d{aJfn(=e?V~^jL$?r%H}9q}~91vuZ}r(tgz$xU_NAY1o*WM9`)%KQDQQ*wv9^Bn)l47!+}5#PV-43!%O7mt2cLw#!}XuryPWOk z`QG(mXFESDU#=Q|l=Jnhi;UJMA8WMfdZ_>WoWVx(S9bMztnK&T)r4WLfA09Ry;rr2 z*Hbj3Jid{pjIr;f^*#Pi{Jpr#`@7SJoqp2v#ihR|eF~{xO6nU${~P*v(2oFIb?~7J z;oO4P3GQGLep(|OsKA#NE(iDm)JaqCN!=RtH`Ec34_zp)ll(ODY#QYYaPP_eO{q9E z@wA8w=zN3xk#~*UG_QJH12}u&j)EH!^$M!np`Y*<@NbBva~w={VvF1*I=&o7e8&%8BEe)3HZHo3`H|BJ~>Gx`Om zpUgM;tMBkSo?@Ek&6Jyay=nd7t1d6fccU+4_K?RVxVLORebo8zqg`C@47tno`pw%r z{lx7@wV7~!9vxwyC$Fw9|LlB!!z80kS1bFwMHlZMRepM{y;u3P`>A!du8l|j5B-aqU%p}Mmik`n>v!K6Ur+U6r=N6E zUtIco(x)(2zcp}z=zl{W5Bd?5!c}jC4-L*Oc%9%5f}aKs6nK*0aumWBh&pN2ds4SX z{S9>lx$>dO>m)yoJR9-_a@~8z{f*+##M2@!@V>`gll(O5Opu%E0>|qCH@+5bV$|md zHwQnCd};WeI0JG8w6xDKI*z$I$I zAp+ls`!}AcpO5jtX%(&*eP-&@pYitrzxAps&TP5$W6ZgP2bGh*eiSQQky^fLcbp774Hil zfxD>sZ{`>C7p_mDz7?Dj@J_X^**?GXeQ$dQ``162JxK#kF#AIN9XPl}ucuAFlDClW zxVh`w?t7cw&UMZPuNL{Eo)&$;c@4S84$=pVzTyqNYsv zF>b^&b&hP+Gw3JALvNsWQug8J=XG7h2$hHTGQmE7&yc1Rq>5^?$6!3&}cSC_ruzM)(t); ze-Ak!pNM={@+X2HWCyWN_@~5I7>D}`#wSmMaj}bGKh@Vfo>dPO<9*`yJj7wp+p!Th^o^d$Nb|AqS&)`NK=XHfJc^Tn?~&tPxpW5T|K{nxzWeZ@yH z|CaeNe<<^f_0amo_>ybnr`d0U=M5aZtes<+FLoyV%y9MA>D2@5*|h*McwUTPgl$9P$InR|4+bv)3hu4!19F4hhhCH3K`(_~Q(hJ~1x# z8a;}i7WPK_g1m=cj$Lc@ADNdZ^M~R$#d_d>BX`0@Kz`cf+4FbK8ew*b^G*IU$2&k? z@bB`;U0qHt{Eo*FR^H7h_d(bl^f~;E{|BEVFZ7WM|4)7h^9X$gPMG+eagYzjNA4Ln zU$l$)o7hv%Z9F4KF`n{tv5%bR@FUQh$Tj>PeDlF$7H^1oD^DJM!QVqZ(Z}d5D8FZX z^m(TKp7?7aAJqH*H*$e~M-F(V&L96B8hnd?`X%+A@xIElL?7W-uw5!q>o1{1|b87W*6VALN|(#oqiM z;qTNhb~FDipXGfR;7{y9uJa9c2tN(E$woVVNFjeN?63F$ehj~A(XU-U&$+_oNA{r* zUPqkoG+)kR><@CpzQC7^kGvp9$SLy1{_)(p&+s>V5a$g1XvOzJPk=|!Qil@rmOVmV z`FqG)=wJL(@~(zB-Z%OKek0z-I^nm3z2MyaWpdQYUwZz69D(1Ycn$gtdjwzeUNIl} zzvv(I73|xGZr2;#K(j_$|CIei;5)=n2gaJpx5<#Co8I@B@jD zEY_)Q|FWO-bNSi-CiCBdKV=Uh?_TkOY|6nNFOb}X-==d=_{HJ__%Zw}>i={;f#31B zhyyWS>^b{`9I-F(CF6&@=sbmAhrF?WJh$#M{EZ*m%sHob-nLG0{uDojp5Q)C`hxM8 zH#FpxdPvDz=wI0(&OO*~^auRb(1)G)5OHGoT=78o3BAm@;0yePo^Fw&;2-!*`V2XR zuX(R{U-{SQpH{hHJ(w4A0!1!kT-|$e-i&yT;)(1F=esYJ8_5Clg2sAq&SSn>57v+I z(I0Ir;r}L{=5-5abTat~{u#TK#oLe<{5JfwH2Fo(H&EOz=NzPSE_Mw5#vZ`u$V=Ec z#la)ap?RR6*$?=WagYzj#}8y&>>Bo|)m|Y-j0c~y@5F&a&vVbMxHl)zjGbXT z{6hSUMtKHtKR|q*_^Zyn*bT-DJEQp{pZp&CjNZi#LHRxY9d@``G}iBPTktpGcl=)b ze)J}O3UYv6=DF1#p^wqS=x^i+znb^reeqxLld(6fJL|%{utQM%I>wE7g8Vr01-SQR zUzj)Z;(fz@;GYp^fnuldf3d^#HED@+W9P9yOXNqm5^FiG)e`wp-(T_968UkvgM7Bi zPntj0{kP}d_Bg^c$GsxofOs4CKEw;Mm;dDFADx4+L!p_5}*RV4v_`;+#Uf93^cJ%auZ`!BS`bKHO1 zVm|Qh*OjAgC9pcn|b2Tz%Rtv;9vM1{^q`x{37-pJJ0zPzG9q^1I62r zOU4aHg?J=! z8r89p=SrNGePQ0M@&X^l{PEAA_-Rn)%X$!pQT!48iT{>54l%ovb$;07hy8;O$=l%` z2zkN|5s$+@vA(P;_5%K7UE!xz`w@Pv^eTJ`pT+q=_p70SJJ;RzU*DPm$*#IYWMy(M1&zBSb!C&YP^c(sz^sn-*kQek4`Xl@l?29P=A>*Qt!{1aL0e(Tm@ySnXk(c0q z^5oGMEpY+l4Y|RdAvfG#pg)nHdikE_zcrmIvqR|ReBoTrv(~+6*csjPnU=#d!F2 z_#3VHtmsGX%dv;>3F9$uDDg&qAM)CA-w1#4d&Y+^$zS2#9sS9BvA@`fFZgHj@5A5I zbI61AJ$#Iv;e14V6+6j$@xI{?aL!Sl3G2>!Ffa5}Q2A|)8~(fOZ`9wYz6|+b-Z3xi zmG}txV_s0?6Uuy9k61t1ca!^}gO}~lW%GWYrCr;l#s3W*-2Q7Rb?AUl{TJ<=2BZ<= z_8&Rw$OAVz^1yLF8#Q*!pl!EF!~b{xHeG)48~pV*2k*XA+n-we-_XJB7L#~~6<*n3 z>cwr(i~sAg_!po1#~U3m;)oF&9zOce)^Gp&+&X>ko3FJ!{rk5!+PL+*|Nd5;Klk0C zyS0n&YAuEiUS^pM^r1D{_l=LY*Z&8drRPW5hFJWs?R(2ErW1$x^Wy&8ynpsT0alKo4gdfE diff --git a/src/tad_dftd3/disp.py b/src/tad_dftd3/disp.py index e686267..1beca60 100644 --- a/src/tad_dftd3/disp.py +++ b/src/tad_dftd3/disp.py @@ -54,17 +54,19 @@ """ from __future__ import annotations +from typing import Any + import torch from tad_mctc import storch +from tad_mctc.autograd import is_functorch_tensor from tad_mctc.batch import real_pairs -from tad_mctc.data import pse +from tad_mctc.data import pse, radii from . import data, defaults, model, ncoord from .damping import dispersion_atm, rational_damping from .reference import Reference from .typing import ( DD, - Any, CountingFunction, DampingFunction, Tensor, @@ -125,22 +127,25 @@ def dftd3( """ dd: DD = {"device": positions.device, "dtype": positions.dtype} - if torch.max(numbers) >= defaults.MAX_ELEMENT: - raise ValueError( - f"No D3 parameters available for Z > {defaults.MAX_ELEMENT-1} " - f"({pse.Z2S[defaults.MAX_ELEMENT]})." - ) + if not is_functorch_tensor(numbers): + if torch.max(numbers) >= defaults.MAX_ELEMENT: + raise ValueError( + f"No D3 parameters available for Z > {defaults.MAX_ELEMENT-1} " + f"({pse.Z2S[defaults.MAX_ELEMENT]})." + ) if cutoff is None: cutoff = torch.tensor(defaults.D3_DISP_CUTOFF, **dd) if ref is None: ref = Reference(**dd) if rcov is None: - rcov = data.COV_D3.to(**dd)[numbers] + rcov = radii.COV_D3(**dd)[numbers] if rvdw is None: - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] if r4r2 is None: - r4r2 = data.R4R2.to(**dd)[numbers] + r4r2 = data.R4R2(**dd)[numbers] cn = ncoord.cn_d3( numbers, positions, counting_function=counting_function, rcov=rcov @@ -202,7 +207,7 @@ def dispersion( if cutoff is None: cutoff = torch.tensor(defaults.D3_DISP_CUTOFF, **dd) if r4r2 is None: - r4r2 = data.R4R2.to(**dd)[numbers] + r4r2 = data.R4R2(**dd)[numbers] if numbers.shape != positions.shape[:-1]: raise ValueError( @@ -212,11 +217,13 @@ def dispersion( raise ValueError( "Shape of expectation values is not consistent with atomic numbers.", ) - if torch.max(numbers) >= defaults.MAX_ELEMENT: - raise ValueError( - f"No D3 parameters available for Z > {defaults.MAX_ELEMENT - 1} " - f"({pse.Z2S[defaults.MAX_ELEMENT]})." - ) + + if not is_functorch_tensor(numbers): + if torch.max(numbers) >= defaults.MAX_ELEMENT: + raise ValueError( + f"No D3 parameters available for Z > {defaults.MAX_ELEMENT-1} " + f"({pse.Z2S[defaults.MAX_ELEMENT]})." + ) # two-body dispersion energy = dispersion2( @@ -226,7 +233,9 @@ def dispersion( # three-body dispersion if "s9" in param and param["s9"] != 0.0: if rvdw is None: - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] energy += dispersion3(numbers, positions, param, c6, rvdw, cutoff) diff --git a/src/tad_dftd3/model/c6.py b/src/tad_dftd3/model/c6.py index a146202..78805a6 100644 --- a/src/tad_dftd3/model/c6.py +++ b/src/tad_dftd3/model/c6.py @@ -23,13 +23,15 @@ """ from __future__ import annotations +from typing import NamedTuple, Protocol + import torch from tad_mctc._version import __tversion__ from tad_mctc.math import einsum from tad_mctc.tools import memory +from tad_mctc.typing import Callable, Tensor from ..reference import Reference -from ..typing import Callable, Protocol, Tensor __all__ = ["atomic_c6"] @@ -260,13 +262,18 @@ def _atomic_c6_chunked( # (..., n1, n2, r1, r2) * (..., n1, r1) * (..., n2, r2) -> (..., n1, n2) contribution = _einsum(rc6_chunk, weights_chunk, weights) - # Add contributions to the correct slice of the output tensor - c6_output[..., start:end, :] += contribution + # Add contributions to the correct slice of the output tensor with + # out-of-place aggregation along the "i" axis. In-place aggregation: + # c6_output[..., start:end, :] += contribution + idx = torch.arange(start, end, device=numbers.device) + c6_output = torch.index_add( + c6_output, dim=-2, index=idx, source=contribution + ) return c6_output -# custom autograd functions +# typing class CTX(Protocol): @@ -276,6 +283,14 @@ class CTX(Protocol): reference: Reference +class VmapInfo(NamedTuple): + batch_size: int + randomness: str + + +# custom autograd functions + + class AtomicC6Base(torch.autograd.Function): """ Base class for the version-specific autograd function for atomic C6. @@ -283,7 +298,9 @@ class AtomicC6Base(torch.autograd.Function): """ @staticmethod - def backward(ctx: CTX, grad_out: Tensor) -> tuple[None, Tensor, None, None]: + def backward( + ctx: CTX, grad_outputs: Tensor + ) -> tuple[None, Tensor, None, None]: numbers, weights = ctx.saved_tensors chunk_size = ctx.chunk_size ref = ctx.reference @@ -304,14 +321,14 @@ def backward(ctx: CTX, grad_out: Tensor) -> tuple[None, Tensor, None, None]: g_jb = einsum("...ijab,...ia->...ijb", rc6, weights) # vjp: (..., n1, n2) * (..., n1, n2, r2) -> (..., n2, r2) - _gj = einsum("...ij,...ijb->...jb", grad_out, g_jb) + _gj = einsum("...ij,...ijb->...jb", grad_outputs, g_jb) # ∂c_ij/∂w_ia = ∑b w_jb * c_ijab # (..., n1, n2, r1, r2) * (..., n2, r2) -> (..., n1, n2, r1) g_ia = einsum("...ijab,...jb->...ija", rc6, weights) # vjp: (..., n1, n2) * (..., n1, n2, r1) -> (..., n1, r1) - _gi = einsum("...ij,...ija->...ia", grad_out, g_ia) + _gi = einsum("...ij,...ija->...ia", grad_outputs, g_ia) weights_bar = _gi + _gj @@ -322,13 +339,14 @@ def backward(ctx: CTX, grad_out: Tensor) -> tuple[None, Tensor, None, None]: ####################### nat = weights.shape[-2] - weights_bar = torch.zeros_like(weights) + gi_accum = torch.zeros_like(weights) + gj_accum = torch.zeros_like(weights) for start in range(0, nat, chunk_size): end = min(start + chunk_size, nat) # Numbers and derivatives for this chunk - grad_chunk = grad_out[..., start:end, :] # (..., chunk_size, nat) + grad_chunk = grad_outputs[..., start:end, :] # (..., c_size, nat) num_chunk = numbers[..., start:end] # (..., chunk_size) # Chunked indexing into reference.c6: (..., chunk_size, nat, 7, 7) @@ -346,11 +364,14 @@ def backward(ctx: CTX, grad_out: Tensor) -> tuple[None, Tensor, None, None]: g_jb = einsum("...ijab,...ia->...ijb", rc6_chunk, weights_chunk) _gj = einsum("...ij,...ijb->...jb", grad_chunk, g_jb) - # Accumulate gradients for the current chunk - weights_bar[..., start:end, :] += _gi - weights_bar += _gj + # Accumulate gradients for current chunk with using out-of-place ops + # to allow vmap. Old version: weights_bar[..., start:end, :] += _gi + idx = torch.arange(start, end, device=weights.device) + gi_accum = torch.index_add(gi_accum, dim=-2, index=idx, source=_gi) - return None, weights_bar, None, None + gj_accum = gj_accum + _gj + + return None, gi_accum + gj_accum, None, None class AtomicC6_V1(AtomicC6Base): @@ -383,16 +404,17 @@ class AtomicC6_V2(AtomicC6Base): This is supposed to reduce memory usage. """ - generate_vmap_rule = True - # https://pytorch.org/docs/master/notes/extending.func.html#automatically-generate-a-vmap-rule - # should work since we only use PyTorch operations + generate_vmap_rule = False + # Auto-generation should work since we only use PyTorch operations, + # however, it does not: PyTorch throws an internal error when indexing + # `reference.c6` with the `numbers` tensor. @staticmethod def forward( numbers: Tensor, weights: Tensor, reference: Reference, - chunk_size: None | int = None, + chunk_size: int | None = None, ) -> Tensor: if chunk_size is None: return _atomic_c6_full(numbers, weights, reference) @@ -410,3 +432,49 @@ def setup_context( ctx.save_for_backward(numbers, weights) ctx.chunk_size = chunk_size ctx.reference = reference + + @staticmethod + def vmap( + info: VmapInfo, + in_dims: tuple[int | None, ...], + numbers: Tensor, + weights: Tensor, + reference: Reference, + chunk_size: int | None, + ) -> tuple[Tensor, int]: + bd_n, bd_w, bd_ref, bd_cs = in_dims + + # Only numbers and weights are batched + if bd_ref is not None or bd_cs is not None: # pragma: no cover + raise ValueError( + "`Reference` and `chunk_size` must be static under vmap." + ) + + # Move batch dimensions to the front if necessary + if info.batch_size != numbers.shape[0]: + raise ValueError( + f"Batch size mismatch: expected {info.batch_size}, got " + f"{numbers.shape[0]} in `numbers`. The first dimension " + "should be the batch dimension." + ) + + if info.batch_size != weights.shape[0]: + raise ValueError( + f"Batch size mismatch: expected {info.batch_size}, got " + f"{weights.shape[0]} in `weights`. The first dimension " + "should be the batch dimension." + ) + + if bd_n not in (0, None): # pragma: no cover + raise ValueError( + f"Batch dimension for `numbers` must be 0 (first dimension) " + f"or `None`, got {bd_n}." + ) + if bd_w not in (0, None): # pragma: no cover + raise ValueError( + f"Batch dimension for `weights` must be 0 (first dimension) " + f"or `None`, got {bd_w}." + ) + + out = AtomicC6_V2.forward(numbers, weights, reference, chunk_size) + return out, 0 diff --git a/src/tad_dftd3/model/weights.py b/src/tad_dftd3/model/weights.py index c6b7b83..d493aa2 100644 --- a/src/tad_dftd3/model/weights.py +++ b/src/tad_dftd3/model/weights.py @@ -42,11 +42,15 @@ """ from __future__ import annotations +from typing import Any + import torch from tad_mctc import storch +from tad_mctc.autograd import is_functorch_tensor +from tad_mctc.typing import Tensor from ..reference import Reference -from ..typing import Any, Tensor, WeightingFunction +from ..typing import WeightingFunction __all__ = ["gaussian_weight", "weight_references"] @@ -141,7 +145,10 @@ def weight_references( # back to real dtype gw_temp = storch.divide(weights, norm, eps=small).type(cn.dtype) - assert torch.isnan(gw_temp).sum() == 0 + + # If the tensor is not a grad tracking tensor, we can check for NaN's + if not is_functorch_tensor(gw_temp): + assert torch.isnan(gw_temp).sum() == 0 # The following section handles cases with large CNs that lead to zeros in # after the exponential in the weighting function. If this happens all diff --git a/src/tad_dftd3/reference.py b/src/tad_dftd3/reference.py index ba277bf..80b556e 100644 --- a/src/tad_dftd3/reference.py +++ b/src/tad_dftd3/reference.py @@ -20,12 +20,11 @@ C6 dispersion coefficients. """ import os.path as op -from typing import Optional +from typing import Any, NoReturn, Optional import torch from tad_mctc._version import __tversion__ - -from .typing import Any, NoReturn, Tensor, get_default_device, get_default_dtype +from tad_mctc.typing import Tensor, get_default_device, get_default_dtype __all__ = ["Reference"] diff --git a/src/tad_dftd3/typing/__init__.py b/src/tad_dftd3/typing/__init__.py index 45c1d5b..51bbf70 100644 --- a/src/tad_dftd3/typing/__init__.py +++ b/src/tad_dftd3/typing/__init__.py @@ -18,6 +18,5 @@ All type annotations for this project. """ -from .builtin import * from .d3 import * from .pytorch import * diff --git a/src/tad_dftd3/typing/builtin.py b/src/tad_dftd3/typing/builtin.py deleted file mode 100644 index aab741b..0000000 --- a/src/tad_dftd3/typing/builtin.py +++ /dev/null @@ -1,24 +0,0 @@ -# This file is part of tad-dftd3. -# SPDX-Identifier: Apache-2.0 -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -""" -Type annotations: Built-ins -=========================== - -Built-in type annotations are imported from the *tad-mctc* library, which -handles some version checking. -""" -from tad_mctc.typing import Any, Callable, NoReturn, Protocol, TypedDict - -__all__ = ["Any", "Callable", "NoReturn", "Protocol", "TypedDict"] diff --git a/test/conftest.py b/test/conftest.py index 6f2c6fa..4fd143f 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -96,7 +96,8 @@ def pytest_addoption(parser: pytest.Parser) -> None: default=6, type=int, help=( - "Number of digits of precision for floating point output " "(default = 4)." + "Number of digits of precision for floating point output " + "(default = 4)." ), ) @@ -176,4 +177,6 @@ def pytest_runtest_setup(item: pytest.Function) -> None: for _ in item.iter_markers(name="cuda"): if not torch.cuda.is_available(): - pytest.skip("Torch not compiled with CUDA or no CUDA device available.") + pytest.skip( + "Torch not compiled with CUDA or no CUDA device available." + ) diff --git a/test/test_disp/test_dftd3.py b/test/test_disp/test_dftd3.py index 778f9c8..3d0d69d 100644 --- a/test/test_disp/test_dftd3.py +++ b/test/test_disp/test_dftd3.py @@ -18,6 +18,7 @@ import pytest import torch from tad_mctc.batch import pack +from tad_mctc.data import radii from tad_dftd3 import damping, data, dftd3, model, reference from tad_dftd3.ncoord import exp_count @@ -59,9 +60,11 @@ def test_single(dtype: torch.dtype, name: str) -> None: positions = sample["positions"].to(**dd) ref = (sample["disp2"] + sample["disp3"]).to(**dd) - rcov = data.COV_D3.to(**dd)[numbers] - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] - r4r2 = data.R4R2.to(**dd)[numbers] + rcov = radii.COV_D3(**dd)[numbers] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] + r4r2 = data.R4R2(**dd)[numbers] cutoff = torch.tensor(50, **dd) param = { diff --git a/test/test_disp/test_disp.py b/test/test_disp/test_disp.py index 5820b2c..4e67dc7 100644 --- a/test/test_disp/test_disp.py +++ b/test/test_disp/test_disp.py @@ -20,6 +20,7 @@ import pytest import torch from tad_mctc.batch import pack +from tad_mctc.data import radii from tad_dftd3 import damping, data, disp from tad_dftd3.typing import DD @@ -40,7 +41,9 @@ } # TPSS0-D3BJ parameters -param_noatm = {k: torch.tensor(0.0) if k == "s9" else v for k, v in param.items()} +param_noatm = { + k: torch.tensor(0.0) if k == "s9" else v for k, v in param.items() +} def test_fail() -> None: @@ -73,8 +76,10 @@ def test_disp2_single(dtype: torch.dtype, name: str) -> None: positions = sample["positions"].to(**dd) ref = sample["disp2"].to(**dd) c6 = sample["c6"].to(**dd) - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] - r4r2 = data.R4R2.to(**dd)[numbers] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] + r4r2 = data.R4R2(**dd)[numbers] cutoff = torch.tensor(50.0, **dd) par = {k: v.to(**dd) for k, v in param_noatm.items()} @@ -147,7 +152,9 @@ def test_atm_single(dtype: torch.dtype, name: str) -> None: c6 = sample["c6"].to(**dd) ref = sample["disp3"].to(**dd) - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] par = {k: v.to(**dd) for k, v in param.items()} @@ -200,7 +207,9 @@ def test_atm_batch(dtype: torch.dtype, name1: str, name2: str) -> None: par = {k: v.to(**dd) for k, v in param.items()} - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] energy = damping.dispersion_atm( numbers, diff --git a/test/test_disp/test_special.py b/test/test_disp/test_special.py index 7fc697e..31d4f69 100644 --- a/test/test_disp/test_special.py +++ b/test/test_disp/test_special.py @@ -19,6 +19,7 @@ import pytest import torch from tad_mctc.batch import pack +from tad_mctc.data import radii from tad_dftd3 import damping, data, dftd3, model, reference from tad_dftd3.ncoord import exp_count @@ -38,9 +39,11 @@ def test_single(dtype: torch.dtype, name: str) -> None: positions = sample["positions"].to(**dd) ref = sample["disp2"].to(**dd) - rcov = data.COV_D3.to(**dd)[numbers] - rvdw = data.VDW_D3.to(**dd)[numbers.unsqueeze(-1), numbers.unsqueeze(-2)] - r4r2 = data.R4R2.to(**dd)[numbers] + rcov = radii.COV_D3(**dd)[numbers] + rvdw = radii.VDW_PAIRWISE(**dd)[ + numbers.unsqueeze(-1), numbers.unsqueeze(-2) + ] + r4r2 = data.R4R2(**dd)[numbers] cutoff = torch.tensor(50, **dd) # GFN1-xTB parameters diff --git a/test/test_grad/test_hessian.py b/test/test_grad/test_hessian.py index abc0c90..829b8dc 100644 --- a/test/test_grad/test_hessian.py +++ b/test/test_grad/test_hessian.py @@ -19,7 +19,8 @@ import pytest import torch -from tad_mctc.autograd import hessian +from tad_mctc._version import __tversion__ +from tad_mctc.autograd import hess_fn_rev, hessian from tad_mctc.batch import pack from tad_mctc.convert import reshape_fortran @@ -88,11 +89,59 @@ def test_single(dtype: torch.dtype, name: str) -> None: positions.detach_() -# TODO: Figure out batched Hessian computation +@pytest.mark.skipif(__tversion__ < (2, 0, 0), reason="Requires PyTorch>=2.0.0") +@pytest.mark.parametrize("dtype", [torch.double]) +@pytest.mark.parametrize("name", sample_list) +def test_single_v2(dtype: torch.dtype, name: str) -> None: + dd: DD = {"device": DEVICE, "dtype": dtype} + + sample = samples[name] + numbers = sample["numbers"].to(DEVICE) + positions = sample["positions"].to(**dd) + + # GFN1-xTB parameters + param = { + "s6": torch.tensor(1.00000000, **dd), + "s8": torch.tensor(2.40000000, **dd), + "s9": torch.tensor(0.00000000, **dd), + "a1": torch.tensor(0.63000000, **dd), + "a2": torch.tensor(5.00000000, **dd), + } + + ref = reshape_fortran( + sample["hessian"].to(**dd), + torch.Size(2 * (numbers.shape[-1], 3)), + ) + + # variable to be differentiated + positions.requires_grad_(True) + + def _energy(numbers: torch.Tensor, positions: torch.Tensor) -> torch.Tensor: + """ + Closure over non-tensor argument `param` for `dftd3` function. + + Returns energy as scalar, which is required for Hessian computation + to obtain the correct shape of ``(..., nat, 3, nat, 3)``. + """ + return dftd3(numbers, positions, param).sum(-1) + + pos = positions.clone().requires_grad_(True) + hess = hess_fn_rev(_energy, argnums=1)(numbers, pos) + assert isinstance(hess, Tensor) + + assert pytest.approx(ref.cpu(), abs=tol, rel=tol) == hess.detach().cpu() + + positions.detach_() + + +@pytest.mark.skipif(__tversion__ < (2, 0, 0), reason="Requires PyTorch>=2.0.0") @pytest.mark.parametrize("dtype", [torch.double]) @pytest.mark.parametrize("name1", ["LiH"]) @pytest.mark.parametrize("name2", sample_list) -def skip_test_batch(dtype: torch.dtype, name1: str, name2: str) -> None: +@pytest.mark.parametrize("chunk_size", [None, 2]) +def test_batch( + dtype: torch.dtype, name1: str, name2: str, chunk_size: int | None +) -> None: dd: DD = {"device": DEVICE, "dtype": dtype} sample1, sample2 = samples[name1], samples[name2] @@ -134,7 +183,20 @@ def skip_test_batch(dtype: torch.dtype, name1: str, name2: str) -> None: # variable to be differentiated positions.requires_grad_(True) - hess = hessian(dftd3, (numbers, positions, param), argnums=1) + def _energy(numbers: torch.Tensor, positions: torch.Tensor) -> torch.Tensor: + """ + Closure over non-tensor argument `param` for `dftd3` function. + + Returns energy as scalar, which is required for Hessian computation + to obtain the correct shape of ``(..., nat, 3, nat, 3)``. + """ + return dftd3(numbers, positions, param, chunk_size=chunk_size).sum(-1) + + hess_fn = hess_fn_rev(_energy, argnums=1) + hess_fn_batch = torch.func.vmap(hess_fn, in_dims=(0, 0)) + + hess = hess_fn_batch(numbers, positions) + assert isinstance(hess, Tensor) assert pytest.approx(ref.cpu(), abs=tol, rel=tol) == hess.detach().cpu() positions.detach_() diff --git a/test/test_grad/test_param.py b/test/test_grad/test_param.py index dfc166d..a9d6715 100644 --- a/test/test_grad/test_param.py +++ b/test/test_grad/test_param.py @@ -22,9 +22,9 @@ from tad_mctc.autograd import dgradcheck, dgradgradcheck from tad_mctc.batch import pack from tad_mctc.data.molecules import mols as samples +from tad_mctc.typing import DD, Callable, Tensor from tad_dftd3 import dftd3 -from tad_dftd3.typing import DD, Callable, Tensor from ..conftest import DEVICE, FAST_MODE @@ -136,7 +136,9 @@ def test_gradcheck_batch(dtype: torch.dtype, name1: str, name2: str) -> None: @pytest.mark.parametrize("dtype", [torch.double]) @pytest.mark.parametrize("name1", ["LiH"]) @pytest.mark.parametrize("name2", sample_list) -def test_gradgradcheck_batch(dtype: torch.dtype, name1: str, name2: str) -> None: +def test_gradgradcheck_batch( + dtype: torch.dtype, name1: str, name2: str +) -> None: """ Check a single analytical gradient of parameters against numerical gradient from `torch.autograd.gradgradcheck`. diff --git a/test/test_grad/test_pos.py b/test/test_grad/test_pos.py index e1693b4..1be8ccc 100644 --- a/test/test_grad/test_pos.py +++ b/test/test_grad/test_pos.py @@ -21,9 +21,9 @@ import torch from tad_mctc.autograd import dgradcheck, dgradgradcheck, jacrev from tad_mctc.batch import pack +from tad_mctc.typing import DD, Callable, Tensor from tad_dftd3 import dftd3 -from tad_dftd3.typing import DD, Callable, Tensor from ..conftest import DEVICE, FAST_MODE from .samples import samples @@ -137,7 +137,9 @@ def test_gradcheck_batch(dtype: torch.dtype, name1: str, name2: str) -> None: @pytest.mark.parametrize("dtype", [torch.double]) @pytest.mark.parametrize("name1", ["LiH"]) @pytest.mark.parametrize("name2", sample_list) -def test_gradgradcheck_batch(dtype: torch.dtype, name1: str, name2: str) -> None: +def test_gradgradcheck_batch( + dtype: torch.dtype, name1: str, name2: str +) -> None: """ Check a single analytical gradient of parameters against numerical gradient from `torch.autograd.gradgradcheck`. diff --git a/test/test_model/samples.py b/test/test_model/samples.py index c4bdfff..c8dd9cf 100644 --- a/test/test_model/samples.py +++ b/test/test_model/samples.py @@ -19,8 +19,7 @@ import torch from tad_mctc.data.molecules import merge_nested_dicts, mols - -from tad_dftd3.typing import Molecule, Tensor, TypedDict +from tad_mctc.typing import Molecule, Tensor, TypedDict class Refs(TypedDict): diff --git a/test/test_model/test_c6.py b/test/test_model/test_c6.py index f27092d..a66db6b 100644 --- a/test/test_model/test_c6.py +++ b/test/test_model/test_c6.py @@ -17,13 +17,16 @@ """ from __future__ import annotations +from typing import Protocol + import pytest import torch +from tad_mctc._version import __tversion__ from tad_mctc.autograd import dgradcheck, dgradgradcheck from tad_mctc.batch import pack +from tad_mctc.typing import DD, Callable, Tensor from tad_dftd3 import model, ncoord, reference -from tad_dftd3.typing import DD, Callable, Protocol, Tensor from ..conftest import DEVICE, FAST_MODE from .samples import samples @@ -109,6 +112,53 @@ def test_chunked(dtype: torch.dtype, size: int, chunk_size: int) -> None: assert pytest.approx(c6.cpu(), abs=tol, rel=tol) == c6_chunked.cpu() +@pytest.mark.skipif(__tversion__ < (2, 1, 0), reason="Requires PyTorch>=2.1.0") +def test_fail() -> None: + dd: DD = {"device": DEVICE, "dtype": torch.float64} + size = 10 + nbatch = 2 + + ref = reference.Reference(**dd) + numbers = torch.randint(1, 86, (nbatch, size), device=DEVICE) + positions = torch.rand((nbatch, size, 3), **dd) * 10 + + cn = ncoord.cn_d3(numbers, positions) + weights = model.weight_references(numbers, cn, ref) + + def _c6(nums: Tensor, ws: Tensor) -> Tensor: + return model.atomic_c6(nums, ws, ref) + + jac = torch.func.jacrev(_c6, argnums=1) + + # Correct batch dimensions + with pytest.raises(ValueError) as excinfo: + vjac = torch.func.vmap(jac, in_dims=(0, None)) + _ = vjac(numbers, weights.moveaxis(1, 0)) + + assert "Batch size mismatch" in str(excinfo.value) + assert "weights" in str(excinfo.value) + + # Correct batch dimensions + with pytest.raises(ValueError) as excinfo: + vjac = torch.func.vmap(jac, in_dims=(None, 0)) + _ = vjac(numbers.moveaxis(1, 0), weights) + + assert "Batch size mismatch" in str(excinfo.value) + assert "numbers" in str(excinfo.value) + + # Internal vmap errors + + # Batch dimensions is always 0 for numbers and weights + with pytest.raises(ValueError) as excinfo: + vjac = torch.func.vmap(jac, in_dims=(1, 0)) + _ = vjac(numbers, weights) + + # Batch dimensions is always 0 for numbers and weights + with pytest.raises(ValueError): + vjac = torch.func.vmap(jac, in_dims=(0, 1)) + _ = vjac(numbers, weights) + + ############################################################################### diff --git a/test/test_model/test_load.py b/test/test_model/test_load.py index 14637a2..c7a7dc4 100644 --- a/test/test_model/test_load.py +++ b/test/test_model/test_load.py @@ -17,7 +17,7 @@ """ import torch -from tad_dftd3 import data, defaults, reference +from tad_dftd3 import defaults, reference def test_ref() -> None: @@ -25,8 +25,3 @@ def test_ref() -> None: assert c6.shape == torch.Size( (defaults.MAX_ELEMENT, defaults.MAX_ELEMENT, 7, 7), ) - - -def test_vdw() -> None: - vdw = data.radii._load_vdw_rad_d3(dtype=torch.double) - assert vdw.shape == torch.Size((defaults.MAX_ELEMENT, defaults.MAX_ELEMENT)) diff --git a/test/test_model/test_reference.py b/test/test_model/test_reference.py index b92c7c6..f498461 100644 --- a/test/test_model/test_reference.py +++ b/test/test_model/test_reference.py @@ -15,16 +15,15 @@ """ Test the reference. """ -from typing import Optional, Union +from typing import Any, Optional, TypedDict, Union from unittest.mock import patch import pytest import torch from tad_mctc.convert import str_to_device -from tad_mctc.typing import MockTensor +from tad_mctc.typing import DD, MockTensor, Tensor from tad_dftd3 import reference -from tad_dftd3.typing import DD, Any, Tensor, TypedDict from ..conftest import DEVICE diff --git a/test/test_model/test_weights.py b/test/test_model/test_weights.py index 38fe4ba..72bfa47 100644 --- a/test/test_model/test_weights.py +++ b/test/test_model/test_weights.py @@ -18,9 +18,9 @@ import pytest import torch from tad_mctc.batch import pack +from tad_mctc.typing import DD from tad_dftd3 import model, reference -from tad_dftd3.typing import DD from ..conftest import DEVICE from .samples import samples From 1640456ddd1dbc99b8903c7c51367b8bdb46c10f Mon Sep 17 00:00:00 2001 From: Marvin Friede <51965259+marvinfriede@users.noreply.github.com> Date: Sun, 3 Aug 2025 19:23:10 +0200 Subject: [PATCH 2/2] Bump version to 0.5.0 (#90) --- src/tad_dftd3/__version__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tad_dftd3/__version__.py b/src/tad_dftd3/__version__.py index a1c88f7..b6f2a57 100644 --- a/src/tad_dftd3/__version__.py +++ b/src/tad_dftd3/__version__.py @@ -15,4 +15,4 @@ """ Version module for *tad-dftd3*. """ -__version__ = "0.4.1" +__version__ = "0.5.0"