From def786d1ac841e8ebdb0bc1357c09e441b759f98 Mon Sep 17 00:00:00 2001 From: Grant Hill Date: Mon, 26 Jun 2023 16:29:44 +0100 Subject: [PATCH 1/6] Adds a well tempered Strategy and an example --- basisopt/basis/atomic.py | 46 +++++++++- basisopt/basis/basis.py | 23 ++++- basisopt/data.py | 18 ++++ basisopt/opt/welltemper.py | 177 +++++++++++++++++++++++++++++++++++++ examples/wt_strategy.py | 15 ++++ 5 files changed, 277 insertions(+), 2 deletions(-) create mode 100644 basisopt/opt/welltemper.py create mode 100644 examples/wt_strategy.py diff --git a/basisopt/basis/atomic.py b/basisopt/basis/atomic.py index 8559b75..7bff3f2 100644 --- a/basisopt/basis/atomic.py +++ b/basisopt/basis/atomic.py @@ -13,10 +13,11 @@ from basisopt.opt.eventemper import EvenTemperedStrategy from basisopt.opt.optimizers import optimize from basisopt.opt.strategies import Strategy +from basisopt.opt.welltemper import WellTemperedStrategy from basisopt.util import bo_logger from . import zetatools as zt -from .basis import Basis, even_temper_expansion +from .basis import Basis, even_temper_expansion, well_temper_expansion def needs_element(func: Callable) -> Callable: @@ -40,6 +41,7 @@ class AtomicBasis(Basis): Attributes: et_params (data.ETParams): even tempered expansion parameters + wt_params (data.WTParams): well tempered expansion parameters charge (int): net charge on atom multiplicity (int): spin multiplicity of atom config (dict): configuration of basis, (k, v) pairs @@ -65,6 +67,7 @@ def __init__(self, name: str = 'H', charge: int = 0, mult: int = 1): self.element = name self._done_setup = False self.et_params = None + self.wt_params = None if self._element is not None: self.charge = charge @@ -83,6 +86,7 @@ def as_dict(self) -> dict[str, Any]: d["@module"] = type(self).__module__ d["@class"] = type(self).__name__ d["et_params"] = self.et_params + d["wt_params"] = self.wt_params if hasattr(self, 'strategy'): if isinstance(self.strategy, Strategy): @@ -105,6 +109,7 @@ def from_dict(cls, d: dict[str, Any]) -> object: instance._tests = basis._tests instance._molecule = basis._molecule instance.et_params = d.get("et_params", None) + instance.wt_params = d.get("wt_params", None) instance.strategy = d.get("strategy", None) if instance.strategy: instance._done_setup = d.get("done_setup", False) @@ -301,6 +306,45 @@ def set_even_tempered( else: self._molecule.basis[self._symbol] = even_temper_expansion(self.et_params) + @needs_element + def set_well_tempered( + self, + method: str = 'hf', + accuracy: float = 1e-5, + max_n: int = 18, + max_l: int = -1, + exact_ref: bool = True, + params: dict[str, Any] = {}, + ): + """Looks up or computes a well tempered basis expansion for the atom + + Arguments: + method (str): method to use; possibilities can be found through Wrapper object + accuracy (float): the tolerance to optimize to, compared to reference value + max_n (int): max number of primitives per shell + max_l (int): angular momentum to go up to; if -1, will use max l in minimal config + exact_ref (bool): uses exact numerical HF energy if True, + calculates cc-pV5Z reference value if False + params (dict): dictionary of parameters to pass to the backend - + see the relevant Wrapper object for options + + Sets: + self.wt_params + """ + self.wt_params = data.get_well_temper_params(atom=self._symbol.title(), accuracy=accuracy) + if len(self.wt_params) == 0: + # optimize new params + if exact_ref: + reference = ('exact', data._ATOMIC_HF_ENERGIES[self._element.atomic_number]) + else: + reference = ('cc-pV5Z', None) + strategy = WellTemperedStrategy(max_n=max_n, max_l=max_l) + self.setup(method=method, strategy=strategy, reference=reference, params=params) + self.optimize(algorithm='Nelder-Mead', params=params) + self.wt_params = strategy.shells + else: + self._molecule.basis[self._symbol] = well_temper_expansion(self.wt_params) + @needs_element def optimize(self, algorithm: str = 'Nelder-Mead', params: dict[str, Any] = {}) -> OptResult: """Runs the basis optimization diff --git a/basisopt/basis/basis.py b/basisopt/basis/basis.py index 7a09add..af1227d 100644 --- a/basisopt/basis/basis.py +++ b/basisopt/basis/basis.py @@ -7,7 +7,7 @@ from basisopt import data from basisopt.containers import InternalBasis, Result, Shell -from basisopt.data import ETParams +from basisopt.data import ETParams, WTParams from basisopt.testing import Test from basisopt.util import bo_logger, dict_decode @@ -67,6 +67,27 @@ def even_temper_expansion(params: ETParams) -> list[Shell]: return el_basis +def well_temper_expansion(params: WTParams) -> list[Shell]: + """Forms a basis for an element from well tempered expansion parameters + + Arguments: + params (list): list of tuples corresponding to shells + e.g. [(c_s, x_s, g_s, d_s, n_s), (c_p, x_p, g_p, d_p, n_p), ...] where each shell + is expanded as c_l * (x_l**k)*(1 + g_l*((k+1)/n_l)**d_l) for k=0,...,n_l + + Returns: + list of Shell objects for the well tempered expansion + """ + el_basis = [] + for ix, (c, x, g, d, n) in enumerate(params): + new_shell = Shell() + new_shell.l = data.INV_AM_DICT[ix] + new_shell.exps = np.array([c * (x**p) * (1.0 + g * ((p + 1) / n) ** d) for p in range(n)]) + uncontract_shell(new_shell) + el_basis.append(new_shell) + return el_basis + + def fix_ratio(exps: np.ndarray, ratio: float = 1.4) -> np.ndarray: """Returns a sorted numpy array of exponents where x_{i+1}/x_i >= ratio diff --git a/basisopt/data.py b/basisopt/data.py index 25d7b17..cc36894 100644 --- a/basisopt/data.py +++ b/basisopt/data.py @@ -29,6 +29,11 @@ def atomic_number(element: str) -> int: ETParams = list[tuple[float, float, int]] +"""Dictionary with pre-optimised well-tempered expansions for atoms""" +_WELL_TEMPERED_DATA = {} + +WTParams = list[tuple[float, float, float, float, int]] + def get_even_temper_params(atom: str = 'H', accuracy: float = 1e-5) -> ETParams: """Searches for the relevant even tempered expansion @@ -43,6 +48,19 @@ def get_even_temper_params(atom: str = 'H', accuracy: float = 1e-5) -> ETParams: return [] +def get_well_temper_params(atom: str = 'H', accuracy: float = 1e-5) -> WTParams: + """Searches for the relevant well tempered expansion + from _WELL_TEMPERED_DATA + """ + if atom in _WELL_TEMPERED_DATA: + log_acc = -np.log10(accuracy) + index = max(4, log_acc) - 4 + index = int(min(index, 3)) + return _WELL_TEMPERED_DATA[atom][index] + else: + return [] + + """Essentially exact numerical Hartree-Fock energies for all atoms in Hartree. Ref: Saito 2009, doi.org/10.1016/j.adt.2009.06.001""" _ATOMIC_HF_ENERGIES = { diff --git a/basisopt/opt/welltemper.py b/basisopt/opt/welltemper.py new file mode 100644 index 0000000..e1de34d --- /dev/null +++ b/basisopt/opt/welltemper.py @@ -0,0 +1,177 @@ +from typing import Any + +import numpy as np +from mendeleev import element as md_element + +from basisopt.basis.basis import well_temper_expansion +from basisopt.basis.guesses import null_guess +from basisopt.containers import InternalBasis + +from .preconditioners import unit +from .strategies import Strategy + +_INITIAL_GUESS = (0.1, 2.2, 20.0, 10.0, 8) + + +class WellTemperedStrategy(Strategy): + """Implements a strategy for a well tempered basis set, where each angular + momentum shell is described by five parameters: (c, x, g, d, n) + Each exponent in that shell is then given by + y_k = c*(x**k)*(1 + g_l*((k+1)/n_l)**d_l) for k=0,...,n + + Algorithm: + Evaluate: energy (can change to any RMSE-compatible property) + Loss: root-mean-square error + Guess: null, uses _INITIAL_GUESS above + Pre-conditioner: None + + Initialisation: + - Find minimum no. of shells needed + - max_l >= min_l + - generate initial parameters for each shell + + First run: + - optimize parameters for each shell once, sequentially + + Next shell in list not marked finished: + - re-optimise + - below threshold or n=max_n: mark finished + - above threshold: increment n + Repeat until all shells are marked finished. + + Uses iteration, limited by two parameters: + max_n: max number of exponents in shell + target: threshold for objective function + + Additional attributes: + shells (list): list of (c, x, g, d, n) parameter tuples + shell_done (list): list of flags for whether shell is finished (0) or not (1) + target (float): threshold for optimization delta + max_n (int): maximum number of primitives in shell expansion + max_l (int): maximum angular momentum shell to do; + if -1, does minimal configuration + """ + + def __init__( + self, eval_type: str = 'energy', target: float = 1e-5, max_n: int = 18, max_l: int = -1 + ): + super().__init__(eval_type=eval_type, pre=unit) + self.name = 'WellTemper' + self.shells = [] + self.shell_done = [] + self.target = target + self.guess = null_guess + self.guess_params = {} + self.max_n = max_n + self.max_l = max_l + + def as_dict(self) -> dict[str, Any]: + """Returns MSONable dictionary of object""" + d = super().as_dict() + d["@module"] = type(self).__module__ + d["@class"] = type(self).__name__ + d["shells"] = self.shells + d["shell_done"] = self.shell_done + d["target"] = self.target + d["max_n"] = self.max_n + d["max_l"] = self.max_l + return d + + @classmethod + def from_dict(cls, d: dict[str, Any]) -> object: + """Creates WellTemperedStrategy from MSONable dictionary""" + strategy = Strategy.from_dict(d) + instance = cls( + eval_type=d.get("eval_type", 'energy'), + target=d.get("target", 1e-5), + max_n=d.get("max_n", 18), + max_l=d.get("max_l", -1), + ) + instance.name = strategy.name + instance.params = strategy.params + instance.first_run = strategy.first_run + instance._step = strategy._step + instance.last_objective = strategy.last_objective + instance.delta_objective = strategy.delta_objective + instance.shells = d.get("shells", []) + instance.shell_done = d.get("shell_done", []) + return instance + + def set_basis_shells(self, basis: InternalBasis, element: str): + """Expands well tempered parameters into a basis set + + Arguments: + basis (InternalBasis): the basis set to expand + element (str): the atom type + """ + basis[element] = well_temper_expansion(self.shells) + + def initialise(self, basis: InternalBasis, element: str): + """Initialises the strategy by determing the initial + parameters for each angular momentum shell for the + given element. + + Arguments: + basis (InternalBasis): the basis set being optimized + element (str): the atom type of interest + """ + if self.max_l < 0: + el = md_element(element.title()) + l_list = [l for (n, l) in el.ec.conf.keys()] + min_l = len(set(l_list)) + + self.max_l = max(min_l, self.max_l) + self.shells = [_INITIAL_GUESS] * self.max_l + self.shell_done = [1] * self.max_l + self.set_basis_shells(basis, element) + self.last_objective = 0.0 + self.delta_objective = 0.0 + self.first_run = True + + def get_active(self, basis: InternalBasis, element: str) -> np.ndarray: + """Returns the well temper params for the current shell""" + (c, x, g, d, _) = self.shells[self._step] + return np.array([c, x, g, d]) + + def set_active(self, values: np.ndarray, basis: InternalBasis, element: str): + """Given the well temper params for a shell, expands the basis + Checks that the smallest exponent is >= 1e-5 + that the ratio is >= 1.01, to prevent impossible exponents + that the gamma spacing parameter is >= 1.01 to add well tempering + and that the delta parameter is also >= 1.01 to aid in deviation + from a rigid geometric series. + """ + (c, x, g, d, n) = self.shells[self._step] + c = max(values[0], 1e-5) + x = max(values[1], 1.01) + g = max(values[2], 1.01) + d = max(values[3], 1.01) + self.shells[self._step] = (c, x, g, d, n) + self.set_basis_shells(basis, element) + + def next(self, basis: InternalBasis, element: str, objective: float) -> bool: + self.delta_objective = np.abs(self.last_objective - objective) + self.last_objective = objective + + carry_on = True + if self.first_run: + self._step = self._step + 1 + if self._step == self.max_l: + self.first_run = False + self._step = 0 + (c, x, g, d, n) = self.shells[self._step] + self.shells[self._step] = (c, x, g, d, min(n + 1, self.max_n)) + else: + if self.delta_objective < self.target: + self.shell_done[self._step] = 0 + + self._step = (self._step + 1) % self.max_l + (c, x, g, d, n) = self.shells[self._step] + if n == self.max_n: + self.shell_done[self._step] = 0 + elif self.shell_done[self._step] != 0: + self.shells[self._step] = (c, x, g, d, n + 1) + + carry_on = np.sum(self.shell_done) != 0 + + return carry_on diff --git a/examples/wt_strategy.py b/examples/wt_strategy.py new file mode 100644 index 0000000..84df14b --- /dev/null +++ b/examples/wt_strategy.py @@ -0,0 +1,15 @@ +import basisopt as bo + +bo.set_backend('psi4') +bo.set_tmp_dir('tmp/') + +from basisopt.basis.atomic import AtomicBasis +from basisopt.bse_wrapper import fetch_basis + +ne = AtomicBasis('Ne') +ne.set_well_tempered(method='scf', accuracy=1e-4) +print(ne.wt_params) +print(ne._molecule.get_result('energy'), ne._molecule.get_reference('energy')) +for i, p in enumerate(ne.get_basis()['ne']): + print(f"l = {i}:", p.exps) +ne.save('tmp/neon-wt.obj') From b3b7146ba935626271903676b396d04f91ecdb99 Mon Sep 17 00:00:00 2001 From: Grant Hill Date: Tue, 27 Jun 2023 15:39:13 +0100 Subject: [PATCH 2/6] Adds ability to use a well-tempered guess at exponents, and removes unused import in WT example --- basisopt/basis/guesses.py | 18 +++++++++++++++++- examples/wt_strategy.py | 1 - 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/basisopt/basis/guesses.py b/basisopt/basis/guesses.py index b66e95c..9a2fb3a 100644 --- a/basisopt/basis/guesses.py +++ b/basisopt/basis/guesses.py @@ -5,7 +5,12 @@ from basisopt.bse_wrapper import fetch_basis from basisopt.containers import Shell -from .basis import even_temper_expansion, fix_ratio, uncontract_shell +from .basis import ( + even_temper_expansion, + fix_ratio, + uncontract_shell, + well_temper_expansion, +) # All guess functions need this signature # func(atomic, params={}), where atomic is an AtomicBasis object @@ -57,3 +62,14 @@ def even_tempered_guess(atomic, params={}): if atomic.et_params is None: atomic.set_even_tempered(**params) return even_temper_expansion(atomic.et_params) + + +def well_tempered_guess(atomic, params={}): + """Takes guess from a well-tempered expansion + + Params: + see signature for AtomicBasis.set_well_tempered + """ + if atomic.wt_params is None: + atomic.set_well_tempered(**params) + return well_temper_expansion(atomic.wt_params) diff --git a/examples/wt_strategy.py b/examples/wt_strategy.py index 84df14b..62203a0 100644 --- a/examples/wt_strategy.py +++ b/examples/wt_strategy.py @@ -4,7 +4,6 @@ bo.set_tmp_dir('tmp/') from basisopt.basis.atomic import AtomicBasis -from basisopt.bse_wrapper import fetch_basis ne = AtomicBasis('Ne') ne.set_well_tempered(method='scf', accuracy=1e-4) From 7916843fba77fb387dc73f990df44cc9d6eb2625 Mon Sep 17 00:00:00 2001 From: Grant Hill Date: Wed, 28 Jun 2023 15:36:04 +0100 Subject: [PATCH 3/6] Adds a tutorial and example comparing even- and well-tempered basis sets --- doc/src/_tutorials/compare.rst | 114 +++++++++++++++++++++++++ doc/src/_tutorials/comparison_plot.png | Bin 0 -> 17933 bytes doc/src/tutorials.rst | 1 + examples/et_wt_compare.py | 33 +++++++ 4 files changed, 148 insertions(+) create mode 100644 doc/src/_tutorials/compare.rst create mode 100644 doc/src/_tutorials/comparison_plot.png create mode 100644 examples/et_wt_compare.py diff --git a/doc/src/_tutorials/compare.rst b/doc/src/_tutorials/compare.rst new file mode 100644 index 0000000..2e36fe4 --- /dev/null +++ b/doc/src/_tutorials/compare.rst @@ -0,0 +1,114 @@ +:orphan: +.. _`sec:cpmareetwt`: + +==================================================== +Comparing even-tempered and well-tempered basis sets +==================================================== + +This tutorial compares an even-tempered basis set for the neon atom with an equivalent well-tempered basis set, both in terms of +the SCF energy and the exponents that arrise from the expansions. An example script can also be found in +``examples/et_wt_compare.py``. + +An even-tempered basis set +-------------------------- + +We use the same procedure for optimizing an even-tempered basis set as in the :ref:`sec:eventemper` example, +please see that tutorial for further details of how to perform this in BasisOpt. An even-tempered expansion of :math:`n` exponents +is defined as: + +.. math:: + + \alpha_n = \alpha_0 k^{n-1} + +where :math:`\alpha_0` and :math:`k` are the two variables to be optimized. These can be thought of as the initial (most diffuse) +exponent and the spacing between all subsequent exponents, respectively. By following the even-tempered example +and running the following optimization with the Psi4 backend + +.. code-block:: python + + ne_et.set_even_tempered(method='scf', accuracy=1e-4) + +We should get results that look very similar to:: + + l = 0: (0.24303219077706045, 2.295167745669027, 18) + l = 1: (0.12199362316666561, 2.180056874861101, 13) + +where the first value in parentheses corresponds to the optimized :math:`\alpha_0`, the second value is :math:`k`, +and the final integer is the number of functions that were required to reach the desired accuracy. + +A well-tempered basis set +------------------------- + +A well-tempered basis set builds upon the even-tempered basis by introducing two additional parameters, such that +the expansion is defined as: + +.. math:: + + \alpha_n = \alpha_0 k^{n-1} [1 + \gamma(n / N)^{\delta}] + +where :math:`N` is the total number of primitive exponents in the expansion, and :math:`\gamma` and :math:`\delta` are the +additional variables to be optimized. + +A well-tempered Atomic Basis can be optimized in a similar fashion as the even-tempered example. A full example input is +provided in ``examples/wt_strategy.py``, but perhaps the most important part is: + +.. code-block:: python + + ne_wt.set_well_tempered(method='scf', accuracy=1e-4) + + +The results should look similar to:: + + l = 0: (0.21397184263141658, 2.129617983476232, 16.3786689342292, 11.79947009544172, 18) + l = 1: (0.11831666016132977, 2.130099585140868, 2.7477960324041435, 10.480981519305843, 12) + +A comparison of the above with the optimized even-tempered parameters shows that there are +some minor changes in the initial exponent and first spacing parameter. Perhaps more significant is that only 12 p-type +functions were required to reach the desired accuracy, rather than 13 in the even-tempered case. The effect of the third +and fourth parameters (those introduced in the well-tempered expansion) are perhaps better visualised. + +Visualising the exponents +------------------------- + +In this tutorial example, we create a Internal Basis object as a data structure for storing and comparing the even- and +well-tempered basis sets: + +.. code-block:: python + + viz_basis = InternalBasis() + viz_basis['even'] = ne_et.get_basis()['ne'] + viz_basis['well'] = ne_wt.get_basis()['ne'] + +We then plot the exponents using a similar process as for the :ref:`sec:visualize` example. + +.. code-block:: python + + fig, ax = plot_exponents(viz_basis, atoms=['even', 'well'], split_by_shell=True) + ax[0].set_title('Even-tempered') + ax[1].set_title('Well-tempered') + +Once the plot is saved, the resulting image should look like the following. + +.. image:: comparison_plot.png + :scale: 75 % + :alt: An event plot comparing even- and well-tempered basis sets. The well-tempered version has increased spacing between the lines for larger exponent values. + +The effect of using a well-tempered expansion is particularly clear for the s shell - it can be seen that as the exponent value +increases (becomes tighter), the spacing between the exponents increases. For the p shell the exponents from the well-tempered expansion +still span a similar range as from even-tempered, even though there is one fewer function. + +Effect on energy +---------------- + +The resulting neon SCF energies are:: + + Even-tempered energy: -128.54704832 + Well-tempered energy: -128.54708446 + Numerical SCF energy: -128.54709811 + +We can see that switching to a well-tempered scheme lowers the energy by roughly 36 :math:`\mu E_h` and is significantly closer to the +numerical SCF result, despite having one fewer p function. + +.. toctree:: + :hidden: + \ No newline at end of file diff --git a/doc/src/_tutorials/comparison_plot.png b/doc/src/_tutorials/comparison_plot.png new file mode 100644 index 0000000000000000000000000000000000000000..9fd27dcb65acc7adcd8050c2d0d19aa5d0e10472 GIT binary patch literal 17933 zcmeHvc|4VEyY?-WQW=_%DU>;)GEbEu4Jbu2WlAz-$THItB^k zGPgXLg%*~1SZkjb>KWeszWaIi_wMi8`}^%bnwERr!*!p-aUREUbzMhGV*~v*dJMxh zXdY5Mfnl@(=+~On@X3r8RUZD3byYj~M!&L5VK7e8VszR_X*!2@@XfSMJ8q*34+j01RzwrOvD*1;#bshL97=2g*_r0Bv zid~l(xF2a=!HKhTd2&*?EYGIF`F%q)TbE|sfi=nX?JrzXPYU}CRdj_bC*evr%=U=; z4cuSxt8M|gRvvy*ra0A}OSJGBex@IFIPLV)v#&kF?MoIWB|OLbomKai)3Rr35NT2} zZcA9G59Kvz4RUfzS{`DsQECzQo2e@3yV%ez&~h}CqaQvXc2^xnOt zE_I45ik{;^8r(j68I-1fZLEk9H|45|6%X$Ulg|=m3mtt`eqzntDM>TMlsS?H-lS+OY=tZHmDL-96Wv;c;Ka?n(-z^1TMc5$>NU8Fh!%rNruzH?tmW?0 zL!xYfT9dEoxJ$Ph+A#|xdyjdB@7udqwJu4qtKFr}Ua_&cnU%%Z_4s-=$sLa~GE$#E zm#UQ~5zgD$*^!^6oe-*3nj4yKsq@`>@ESwtw{PE?;X<8o8`RRH0*S6Cj9hB=Wp-G9 zsgT>XYnM|`w0=34*ATO!*JMa$K2fe#rYpp~FA>-3oy>`Amk3iBd*|MN_nLWW5OaNo zKE4@Nq3iB7uHc}cm!x9(Nw}f@mR#%g+)4^t#Z3q%VZ}akXVb#*r-&$aiNf4r}m?c7tTJx`*)Z$Vv{3NO3CQsG4^YVsU{6%19J zZ00Ku=jIMqrntvv_n?nXJlbqg9&$Y>h?_cHA>lne(*HR_-!8jst5dsEMcG(SmHdu) z2_h?F{0>^B$asfBp&LO%u*!04Z@JWD9GmDf*?a%Dz~NH>YT_Riw?W z%|7>rc(KgOu}YELL>y6e@XJfwVBLaD31tGG`<6jo6iyst=wP#Oy6=;&b`6B=Er=GI zmIPUs9FljRyg{;(PtGkqjdRD29m}k1GrPu_V~qQiw z6L!b%vDR|yye@DKqEEgr!R1Uzaa!6=`H@$&xei|nt*8qn)TyCw8TeWT+3uTm%}>n? zR4t{yxaz*7+F#ylg!^&{{tYKSujuxU=bB}Dvajc4r3l`V`PRZ@eu=%09jR!l)~&;5 zO0+6xDB_6N`@?A`+=j9%&QbfQ^T~Fz;{z@JG_*M=UcDf6&y<{lU}<5HANtVSr>^Q& zug3B$zJ}E?jOmD37vX~tDQNUH@h{j7hn!w|V`;Oj_GW9BD*06*`avAK)Z+AY@ zm}!_lbW4NFhD<;Y&U{O8wm&tqf|obr@#9JpkDlm8i;{(5$K0TIS;e#8w5xp@4Xmx$ zXX8Kcg)2m;sk#kj=0BJ9g-FHYysvS(z2oJbgIN0FY4^o;7y@y!LYb01F@04liT~K0 zC^c0nkFoAX3rNY9fVG=)9zT9O1nWNkDeij-l_FJFG~Ib47xKwcylDJHd{PW0uO(BkyWb8SRW?$TqTc;F1(O;ih zGpcR4CbyckBus)!DQ;)#hyCcw?!R89YxQo&b@(tdJHZPqg_qG8Cjb4s4+U0p6Qjb2 z3WwKFc4md|SY&eP_m?EDm&4Bu2@{QmncUYnoDZp%dg3zK<=yPob14)|kj>s+RijcC zx<2chE_}jw2`JBKe){yO6r#P6|Lw@vx2cXD#e-H4Gcs-*PgQ^BTioG0_pG>gIKRWU zsEfVQZL-BWpF+BH#4P^Kolbb(`IB+zkyxc2+CA}u{OvPhTUbm?j7;xsi`EdeHoM$oLpfDq>;!om4asu$uCI3 zN{)n4_4qofr=*EBxklRYi%b2dPn#V}7Tn@1=omOyp_EXYn?}c~#5H_;eC(cm)=L&X z2z5eL^;2V`!^m{Uce(+H=ph2R);`SV>pIh*AjS@sZw;2MAHJBF`p)!Bef-sZsWGdF2Zp3Zvtz3?@j_Tx$6X!G&@`q7cMj2cId98rK=$o{Nn_>5=(WvYJk zEjML1@4aWNC|aFbtYXz-tg@l7*p19NJa0M2p7@%TnjA7ZG9VP7!wOy^H3 zv-B5IW*63Q`;6C5_ez(v%r>J;fg;D{*?;$R-udGkPNmGuqhC))Ss5VQlls7d;#reP z+w3$t8?+DNr`T1lPJHSiK}m{U<$zYU;Y%o<@gf1AKYz}yQ<^K?*JQK_K#IezaMOiv z`l()1-zA{hO=LZ=GS~{O*i>q^kA3;6^#Q4&*E;%Uo%j=+H^gw}oIn_$hyqytu6sAc zf&oeHK}|=hE>XtGyd9!pG@Lpeo_+T9FRXjc9YwXn#Qkd^0-lB)xJP^l1M%_}kA^yy z7C3d?gNrZXj1KBTM+KX80(0)UnqccB4}J#a`H^P3+?vjwSktiN`ECKH&!_aul;__n zyB@k4wD@QcO5BnB$AT9o+r*)!6cCY3 zuT^~eoYB={og51JDiR}Ryp4AC`XF(ifg=*_vA(!LPrpGl%axeLAk8R^q$R z39sdb*Y_6!TFMLKQlw*ZHf&v?qfhIdUzP~DG7lFGtIpT#dSZCyOipWSEA*u!VlOPO zbl@L75(0cz96Ov{!DVpi(xvA!z46v>uQAp=&`wVJ+|-n>B+i?YTPtT(6J$|DZHhqpd!H$C#z0C07V@4~DV+={gu(WB2wMDNLYaxVYT z+fRv8M4jhQPj&&On(LtY4rasItWfWxj9gRBE2eff8agH#yVk4O6RRHcNkSVff=b`A zipA(A0!gUJuU^l3Gn7kF?j=<5RzlrpQ_5JZU1#Za187gVckbLle{bn5yYVH0<8&SR zB6M+Bz6q|Ea3;V7vlD}{5dQfdL+Q#3HqWV)IcszSfL9`iq1?3?SSLoK#-3R>BkCw` z79w5#wp`e%CB1oJ-3J;ixq!g{|wZWb?tlnU#l@W@3@w4a+^D z10qIHNuHd}=7oD;0aQkk!4ie5a@?oB;4D6*?61!>Jmc4}mq;t5=|K>e4t)A16OIKS zP7b5d(L*y0+@KhV)yj%McOt6O-PD2JEO4+1Z>s7&UG=fVkDAbcxn^Cxdi7yurakA} z33E6sVr=O!!+g7kepHGUlfUKt0?Uu?jQ0i0MvL1;->v!|FzhRiF<{tzt7Ez@L}vO- zTy(>p6z>ojOvT{yvBX$hzyU)!uLDG8|A5nyK;YaI)L%p6pL0gi!idO>O$Pyy5f!cI zjFp+a7yQRGMs>>|zJ$wkuC>X$YFB=UNs(t+lGk@OZN-U`#BAU`Ru+(4!XTg%QKjJk zK(N>=Ap^Ql(Lj}#bRMp-V)GjgACgzc4M?+sfd}zt!iS-@%kl$Pa>D%m`7mD!Ik_T; znXePtTv+DwR^ehO@8uh$+i&7Tfo)-?Hw!cInfP2V*=lPMYwUI%rn9m*>m;_0>2m4L z-grk!Gv3C|E(G4Q6gTFoSnh>Y?&qmK4xu#!kWj*xx3_28M5`mV{JnaFDj@Cc!-k9JzGnk=g@f0L6>R@DIB2L0!!)r)g67@j zQ7(Vyi@I?fp=^8jO}t0Tq$*4Rg*HP$P@1$w^KM~bVQwnT%{^U;SkJHU2*nebj%C0H zaP6_K4beC8*y#ZHvVbx>*!3j0wo-l@yeKSiw+5d^SO}_Flpa&c&Z`WZ+Yl%=0vzx3 zX01Rs=LU#Q|DCC_@~Ei3;|s5r2SV=Vc6h`8=PO|wtG&vXpyH(!1nk$}zR5o}$sP#C6p93IxI`=5ANm%a9_dBAenD*;~2Cs;sK zT^}iTlx;)g42&su3j_^DYAl?)=dMrp6OC}L>>8=IZcbVV%bO6EUNb$hsB&im+^AS9 zkib+Gv#UyeBZTN0{WUImFa(tWH)hdfj0;s>7|=$zN9I?5+Rhh!NjIQH+K#^EL~8G9^&!gl!dui z5q=f1ZQHipd+>n&x&&S%eClTiivIdbRQVi%&57O;`rGI=<#`@{t%T|9x<^@AVZd2Q zK3K7x(b$84kCziCtqFYYIvGIV!`5b7?$zi^@^YAOD_JlUH7?xiyD(^4pFjFdX`{T$1;R4wb*r`OCs#& zQY$OqP&Cm9%x<5$Pt6m9T*%-D6F@XSkJMdplMXdr9T0P>C4 zKROutvtjHQABcO0uX( zPS@g!b4}Pnf*vl_1ZTCC+Xa4L5ST*D`#;>Z51SZh5P071Ts08}iSN4iN6eVLR)Db^ zS{GDUkgu&@S^w4{-Ke9Tm<@hY4mBJ|d4nt2|<@{Z4+xO{OF22m~wz;*oP9 zRq`t^?&Aph+TKI~q*I!(a)4{`Zj09)qZY`RD!YrvKzhLN7hsYOTHY!ztD}D84~Q&6;&T zDtR4-?t6Z6QTanT^EL_EouGK&%{NceA-^lD{$7!;G*-Z{hoZ+a$XYsfr+mh=R93;; zf=~9)^(3Y=Klt32=uaQv%`^Y&6~;~&q+KAw`BwvpW@wkKM!bbRJlnek6V{9f9QlIl z2t_=))7)@=#U^b9X27)|Wq`EQ3?Q=;$guDRZDmFEFpfCfNle(^V^!tQn=Kv!{$hMD z4!D!D=yjo-D$wi5)MW+j8<@SQwRMP@JP^U^&!rYO1NG6YC~a&pPp?hF?rw` zp9LiDg$oz%6&EYbPLYyv9R6CTN+s|{T*)pek2KpZM%)YcNR-8IUU)|9y2HX`X2*h3fc$WtE&&S6ASV^|QS?uS>)1m}v_!`~ z244Jfm0I%b!so|G2Kwy>&eGA>gKkS&9adF^#hU((D`5w^+QG+exb<1~`TBd!z$3XF zrc+6iz+*a(B18jju!x5RF;XLdl&e!eb?Vfh1^k8x!7@NF@%jqCHk?0`i2KFdMWa~D z%7iB6#7=`XETX5Gn2|!}Y;@`KDXGT?^mn`H5Z4CVpYb8CE!<0{$D)PbA1J6jCjBJ} zdTSvQ`hVTDSXW+78;+#Z46+C+Awpg?gph7j;S963J6 zB(-l}8BpdhP^++@vweiK!F>0r04D3P-##p|90ynu4&$H_7tDFAQg)~9fHAe_DgpY_ z`S5C$oszIZHW6Y0Ajy!-K?w;93IYg}Rb%0>5_|RF$4>`|hh)TIx$KQV`^*MKH7E#$ ziTSrrk0gBOr;S`_Mh3bS@$ikqo{VNh^eneQ6PgFI@OpZZr!Z1xrUCeI)G1TF0UW^$ zy2Er0w@3G_>DH#=t*Ppv!A+)qioCoIFM4857hE(;rw%Q~G>9gpVlRg_nohrA>HK&^ zJQOK_FpUuhl-d1jZ}kyYVDZ8OJ5v!VwE)5dO^V!}WC-8^e82;~QM7uBv+Q;6YjQn7~Z| zld}8$I|tu+qyap=q}$g6Jv^p1=_YRFb{@@SvtpoUoIu~m!Vv((>SY-j#HytqzAK%NzUW zE&ND2YIdsNY~wJ6WkiI$)q9;&r?wotx*oJjMCCyPWYagkbPAE=?A{b$mY8%O9JN^> zpzOA_t%}&iIJOkeG{N*V+6%mKXGVb?6N0K`4#MDQBFS3-l)I;)6s`>$HX!wPVq&5m z7$>cx>u#1DW+Tg4m^_MVV)(!u;?Wa1zcX%z+*dI_L+%QdCI*4_A{7pato-d$vJZUX z_HBwTWqvFH-OA3tNoID$iygvF?p`hnGC*J4md809O^}fw31Y4^izP0fK->KYm`n9= zu0-5=tZ>cp(eB?Kr6c7g=64_NOeK|2eU1DUCPtyPls`3lYYtSmEqQv2zR@jRNlf_G zVs&uknFq0j1xg@JRaIfC-Nl~Eg%87|qL+Va^z`h0Ae8zKREuCI@7$k49q_*BfJgu{ zvxt6H2)iRe!k0SNp`og3Y;2t1Tqy+06zTHg^e2{(nMso8w?}B48~^rMWD*)#P>}oB zdx5T^DdG*pjw8anJ7FXeSqq@60VS*fv?e)>cE)D0&}Mv##Qnl>d|U#*r!V?miJO0z zWwl81{dHey=s1{Atg<^4@OX|-kC@wl-k1Ov8wHBK{4Jw$g3BB2-q=-2k~N=W-Hv6{ z$?v}=ID5M(j(oTtnvjMhG@*}2)RW;GQju<|svwkiLK57765D-yz5ij$fPPPx6jyr3D? z+F%b|Q3R=>2pTSDcMQ3E@R~ZfIFMc<;oARUy3dc&2TXv@UG$zK?&F!X66N+~f$oLU zEM(h0VP0OaW@Jw{<`-{4E5WwsTz3MK|6%jRCn>si5tPAEkk+Xtq*_Pnm^gJKsFq>* zL{E`)OZK|i3)^?9ag*g-OoV1_>G}_q6UH=W>A%u(dp~_GP0yv|?LPSBD|Xf*PCw=7 z9lbsPMuaV!;2>&uu0= z1X6y1y2Gc`K2OasT@a>xCsD+}HM8r?W^}gzmn0q>jr%xDC`oLmIeT#FA=iCHpa;6V zwyNsjoZ1G~-P?g?XR1@2G8!beOLF*!nweadd~s2YxHhHgP$84p>8IOacyi_G-6FWj zt2^Z_bR|Y`-)vBF`*NwN41|^Ay{a~b4c%AQ4G%zfz>Az*Nv`r4|gvi zB>gpfyl7Yh>4?gT+6C0$QV{Jpz^PzTW;abb?5CI=@X#lQX}?v7N8&#G$v3tgr_1Li z*+_S7Lp{`D=Yz=5R?Bp?6}Lc{HAbX9qD4>Vo;x-^V3F(mA`ZKp^@=d|9{3!n-x2iy z#3~*hPhC_#;U2}3+&izICb&DT2g?SECqfsEj@eRaS@-pnzVR02xt6+4NZv3ct9`h8 zuUd~ zXJ7fdci+S0h&K~Ib!62)ZGX;VI;T0-)8KedmUyXJ?Q&>J172 zsZ>?vS^>(aLs95iW7>~d!NdZopE7bAnUj774)1OVj~gATA~IyJCgh;GtDr? zcP3K?k~@R=!rsYCVyjT)jo9g>8wXh;{v8?+!jJQ#;eH|L1C;)z^PMaf?FCM#EE4z} zHnTfij@s-|Jvii%_72aPlFCM|6lgind}`W?26U#|v17J$A1$CAbuUhbBRqd)3Pv+( za3^q@*HQZ$g|eohx@b>>>n<`V3CqoAV*yME(nNO)6uhXK!Fz@wVdIUhyMmEX0FcIQ zoB}2c<=+RUc%M8njgHpRkk)OF2w23k00?a8orX7Tp6oMu9f|qrmdVIt^XKtQ$il?B zS|oa`Rk}Q^uOt4a#vo{=#@&LQ=apu@_dRi<_iy` z?S}b4D8Wp<)GF3^s~>gBj~n^^kUk8pkVK`BM(4?O6`&g>xfHrCliOq#i-#P&(UZG= z5}+YLU|W(zpiWz1R23B&aJI~skc#Kjr@nt#Zi{(!!F2weyLtWYBWWD~`b3!cnLI%l zWlK0Vrd61RT+S+mK;NmeM)@^3@HiJaYW0otWvquTL*(%~t>$8`;Bcv49d2!~6?s?U z8Zey!REp>Wck`n8XR%6v!N)!fj$l-tzo_*Hdw3FsU?98{I!;%e@xdl>92ylc=yTdOWrZTk9WO7^q!8zfd@G~ibU>^Z zR26kVWVJTZ)jh7iYrfZVaxC+FVnTA(Mid_?_>Rv*-nT*JJjgnkiBGzVZpH>J`CXJY zzd6UzOh+G)-kSJWmoq8GS`4jCtF#z|Vf89YZ?ER$2N0M{o$u=c7s+W7J`W9Rh&({i zym*7xxP9G`vr0D2F1aM>mBxe5-O$(z-}(M|w0EtNS4o2D3*qL|q-n4=bCFcfz`Xwk7l)H(J%l&LR)fsXxpou4W{|=G z4^dEyS}vpC^g&fsFb&FhPJK6?To2}&bI8(@eevrXM9(QrdNe#Zc4MF+9gk}VG`uQ0 zC;dKO@0Kpbn>mNsNH2fsrvDlxq$@T1Lf^t z%-~Ey#+rc^tXT%Ioommi?+;z*lykbZN>~KD(^k>bH*o=ah_`PC-g^{JoyFJVO-rW9 z$yBgH0VY7s>x2fHrxL+?@h5}-NKW6pPm;Wo56&_vl+7tu7=j^So=+f4+Pg*V+rK-D zX{)Arn>5XKUfyMyNokQsj>@Y`+L)5DJBeoVHRd>ARqy*uMq${8`k$VX%~ut7VwjEz zJx`#S_Sz*u;va6x$iK9*@alulXh3+Zkx$?Vj1q0CE z&^WM6lK*1>&|@!vfluQHhtxB^UXYuQ0aT?w5k5cl87$?=Z!5I|v$6{3_^Eg=U8>Pe zht=UAQR&DNL{cMlJfFugx$Pz`EZ=W}@m2JJiic1DETl79$fYa%VINEUbN{TnHkSbD zWj@v>uEqRUoLF-?%d9sAz++q9EMXIdt+MVK&9y{s6?;$b1jj~ zb36-p0v&^PgF+v@nx|{g9va~r@<;OP64P2%{`O>xA0E86M2Jtm+kY*4>D#~8NLs*V z6!rPegI+KPoFNHAoX)U1!3sY9gG2#p?9c03J+mw?0EKph8J0&-gTue06b&tC`0IR_8dfR46;sK6k?T6wZgA(= zS@7I|Ii3@3KSA3pz_$guQ8$Ru2v~Jhgx7(oM8(R<+NOUuTJPvG#P5T(3b|~`%gY<< z*tczSK1vHrJrITr@Pjz$9Y>Ru6iKX@YYxze_}^<)+CpGqCov(0J+!c@|0>1Y$k34e z2VbVO_9h-Vv>rXuG_M}4>{|E}yc7zZ@*oR5q}PS{7|N7K1(p)YMxqg<0~IhsFEbdN zG-(rLWT8uf>8IKJ#x8T{hjlP-xhDOf9R&kqfyPPi1h&gz!Mid_2#6A!kM|La2EZ-F zrG<}9heN%agVXiHezSE@4sOOecz{reZRh@dy!~CBeF~KEWnqgK~Pd=gTZ<168Dkfw{$ zjs5!mzn594C{$kfQT+d_jAj2&_nvqsk3GqK@JJRoj>eajcoA68s=u_L(NF$Vppg4z z`sg`^O8NcSh~lHcbX1mW%ucox+r?=U79wcVd3WN=o_ukOHU8Bx zD_FxteHBt7v^3>!>-&}4ypMVe%{aYWr^9PGDc;On-4Ehr8y-tc=Ei;27071IXT-3FT>GLvD?o$Cd4Y^)w+j(x{! zJTwINI=)>Q70I0bruxo45#$*L;zQ0rpNlm3d!vQQ0a3P^df0=4*rri#>L|B;7`4K0 zUIRSWQtb18a2%S1J;b8oes~{%PC?tmkTA$`>Ca) z-srD@aNG}P9KBgB(3~5QUmK}}U~ObaYYh8D}qwSS#l&VuOJ6cip;Nm;P}a{u5hdQ3+7Fwj`w@ zhJOW?;=lPoitTM#L70}MGL=XxF%=rXE9F1Ce@PbHh^Y1L&FWu1dB9u9OUnhjHlWAQ5go~g-hbKA<4w=%SJSV~SaKgLzgixb1XLlD{HUZ_cT}aMnTT%J`hba-6YvBVyhQEw3{spV%+0BWi9Ne2(3w7?) z#QViMN6+(NUo}S*h|%$S$;&HV<`Y--!7)F%8S!mJ}PA#XkH0fo&8(g zwS76iZ!>l_@^3vtI_vcURhB!*(%=4F?}scbU}fbo=tqqARb<%^G;`hrwe&ixj@B3C zV8&FGjKCL;j=_S_HtWhtE?NK1r+*9e{%-62&n+JELSg+(|6$2PT#WgtLvv5)%{Cz* zyVP`b#PkC5>K-$p{@}WS0-NNYzF)$zr+7Y?f!P`DM{@!Vyc{-aZ9)bC*bz$Lqro;O zqXVEChCyLkbUVWVdqBj)RDOXXSPx_zAFuiM_7^j&k&EU&liqEBTWH7ZqNrN27g7zm zngGBpLN3^Fmgfsqua|SJzt&kab#8vc;D&_JM<*{#bp^7ROn>U2{<2Oi2soTZGzjC) zdSbzEkTdeMWb1261$CO+>7%B=v$+K1qSW=H0^tBYCdwyMMir~wOX!Y{9gK+UIb?PN zy2G&KLvW)xR#lydaj>e2LDjJf+}I!r^BlaY3Nki%13)KN4}g;ov=X;TcX-22MZ}!K zqrLuR)goCR1^%!@f$^2y=iy;*S+#g_S8j*7gW+G$?fz1gfqe0=sVZfSO)A)SdCNKNA;D{crasZVXO-rFer>!cb`@9o_g7_f%KiW%&NMtVK=cMS}3Qa((Z*&5#qke|`@+#K1EdNgNe!N8KN1A9Cc3xtnjxWM(P1ka73XlzNm>Tj5jz Date: Wed, 28 Jun 2023 16:10:54 +0100 Subject: [PATCH 4/6] Forgot to run black --- examples/et_wt_compare.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/examples/et_wt_compare.py b/examples/et_wt_compare.py index 22403ee..9a1305b 100644 --- a/examples/et_wt_compare.py +++ b/examples/et_wt_compare.py @@ -20,6 +20,7 @@ viz_basis['well'] = ne_wt.get_basis()['ne'] import matplotlib.pyplot as plt + from basisopt.viz.basis import plot_exponents fig, ax = plot_exponents(viz_basis, atoms=['even', 'well'], split_by_shell=True) @@ -29,5 +30,3 @@ filename = "comparison_plot.png" bo_logger.info("Writing exponents plot to %s", filename) plt.savefig(filename) - - From b8dd34e74d83538ac7a62f0c8c112542cb1b5188 Mon Sep 17 00:00:00 2001 From: Grant Hill Date: Tue, 4 Jul 2023 14:46:31 +0100 Subject: [PATCH 5/6] Adds a test for a well-tempered expansion --- tests/test_basis.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/tests/test_basis.py b/tests/test_basis.py index 773c3d8..7054f3d 100644 --- a/tests/test_basis.py +++ b/tests/test_basis.py @@ -40,6 +40,20 @@ def test_even_temper_expansion(): assert almost_equal(p_shell.exps[11], 474.989023, thresh=1e-6) +def test_well_temper_expansion(): + wt_params = [(0.1, 2.3, 12.2, 8.6, 13), (0.2, 2.1, 32.8, 9.9, 10)] + wt_basis = basis.well_temper_expansion(wt_params) + assert len(wt_basis) == 2 + + s_shell = wt_basis[0] + p_shell = wt_basis[1] + assert len(s_shell.exps) == 13 + assert almost_equal(s_shell.exps[0], 0.1, thresh=1e-6) + assert almost_equal(s_shell.exps[10], 1615.706126, thresh=1e-6) + assert len(p_shell.exps) == 10 + assert almost_equal(p_shell.exps[6], 33.623095, thresh=1e-6) + + def test_fix_ratio(): exps = np.array([0.1, 0.2, 0.58, 1.3, 3.0, 8.2]) new_exps = basis.fix_ratio(exps) From 49b059c181cbd3a16d5c8946ccd2ed54b584d9a8 Mon Sep 17 00:00:00 2001 From: Robert Shaw Date: Fri, 27 Jun 2025 01:51:23 +0100 Subject: [PATCH 6/6] Moved minimum well-temper ratio into a constant --- basisopt/opt/welltemper.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/basisopt/opt/welltemper.py b/basisopt/opt/welltemper.py index e1de34d..61a8808 100644 --- a/basisopt/opt/welltemper.py +++ b/basisopt/opt/welltemper.py @@ -1,4 +1,4 @@ -from typing import Any +from typing import Any, Final import numpy as np from mendeleev import element as md_element @@ -10,7 +10,8 @@ from .preconditioners import unit from .strategies import Strategy -_INITIAL_GUESS = (0.1, 2.2, 20.0, 10.0, 8) +_INITIAL_GUESS: Final = (0.1, 2.2, 20.0, 10.0, 8) +_MIN_RATIO: Final = 1.01 class WellTemperedStrategy(Strategy): @@ -143,9 +144,9 @@ def set_active(self, values: np.ndarray, basis: InternalBasis, element: str): """ (c, x, g, d, n) = self.shells[self._step] c = max(values[0], 1e-5) - x = max(values[1], 1.01) - g = max(values[2], 1.01) - d = max(values[3], 1.01) + x = max(values[1], _MIN_RATIO) + g = max(values[2], _MIN_RATIO) + d = max(values[3], _MIN_RATIO) self.shells[self._step] = (c, x, g, d, n) self.set_basis_shells(basis, element)