diff --git a/docs/ase.ipynb b/docs/ase.ipynb index fafb6db..94facd0 100644 --- a/docs/ase.ipynb +++ b/docs/ase.ipynb @@ -12,19 +12,10 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "id": "81f0cb23", "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/giulialuise/miniconda3/envs/skala/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html\n", - " from .autonotebook import tqdm as notebook_tqdm\n" - ] - } - ], + "outputs": [], "source": [ "import numpy as np\n", "from ase.build import molecule\n", @@ -47,31 +38,17 @@ "- `basis`: Basis set for the calculation (e.g., \"def2-svp\", \"def2-tzvp\")\n", "- `with_density_fit`: Enable density fitting for faster calculations\n", "- `with_dftd3`: Include DFT-D3 dispersion correction (default: True)\n", - "- `verbose`: Control output verbosity (0-4)" + "- `verbose`: Control output verbosity (0-4)\n", + "- `with_retry`: Enable retry mechanism for SCF convergence (default: True)\n", + "- `ks_config`: Additional configuration options for the KS solver (e.g., convergence criteria)" ] }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "id": "37d38241", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Calculator parameters:\n", - " xc: skala\n", - " basis: def2-svp\n", - " with_density_fit: False\n", - " with_newton: False\n", - " with_dftd3: True\n", - " charge: None\n", - " multiplicity: None\n", - " verbose: 4\n" - ] - } - ], + "outputs": [], "source": [ "# Create a water molecule\n", "atoms = molecule(\"H2O\")\n", @@ -97,120 +74,10 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": null, "id": "3de5e444", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "System: uname_result(system='Linux', node='GCRAZGDL1498', release='6.8.0-1034-azure', version='#39~22.04.1-Ubuntu SMP Wed Aug 13 22:25:47 UTC 2025', machine='x86_64') Threads 24\n", - "Python 3.13.7 | packaged by conda-forge | (main, Sep 3 2025, 14:30:35) [GCC 14.3.0]\n", - "numpy 2.3.3 scipy 1.16.2 h5py 3.14.0\n", - "Date: Tue Oct 7 15:43:37 2025\n", - "PySCF version 2.9.0\n", - "PySCF path /home/giulialuise/miniconda3/envs/skala/lib/python3.13/site-packages/pyscf\n", - "\n", - "[CONFIG] conf_file None\n", - "[INPUT] verbose = 4\n", - "[INPUT] num. atoms = 3\n", - "[INPUT] num. electrons = 10\n", - "[INPUT] charge = 0\n", - "[INPUT] spin (= nelec alpha-beta = 2S) = 0\n", - "[INPUT] symmetry False subgroup None\n", - "[INPUT] Mole.unit = Angstrom\n", - "[INPUT] Symbol X Y Z unit X Y Z unit Magmom\n", - "[INPUT] 1 O 0.000000000000 0.000000000000 0.119262000000 AA 0.000000000000 0.000000000000 0.225372517068 Bohr 0.0\n", - "[INPUT] 2 H 0.000000000000 0.763239000000 -0.477047000000 AA 0.000000000000 1.442312677587 -0.901488178545 Bohr 0.0\n", - "[INPUT] 3 H 0.000000000000 -0.763239000000 -0.477047000000 AA 0.000000000000 -1.442312677587 -0.901488178545 Bohr 0.0\n", - "\n", - "nuclear repulsion = 9.08829376913928\n", - "number of shells = 12\n", - "number of NR pGTOs = 38\n", - "number of NR cGTOs = 24\n", - "basis = def2-svp\n", - "ecp = {}\n", - "CPU time: 1.87\n", - "/home/giulialuise/.cache/huggingface/hub/models--microsoft--skala/snapshots/9eedaca41deeaebe197bec77d1898f777770d3ef/skala-1.0.fun\n", - "Downloaded in 0.16 seconds\n", - "\n", - "\n", - "******** ********\n", - "method = SkalaRKS\n", - "initial guess = minao\n", - "damping factor = 0\n", - "level_shift factor = 0\n", - "DIIS = \n", - "diis_start_cycle = 1\n", - "diis_space = 8\n", - "diis_damp = 0\n", - "SCF conv_tol = 1e-09\n", - "SCF conv_tol_grad = None\n", - "SCF max_cycles = 50\n", - "direct_scf = True\n", - "direct_scf_tol = 1e-13\n", - "chkfile to save SCF result = /tmp/tmpcw83w2ua\n", - "max_memory 4000 MB (current use 736 MB)\n", - "XC library libxc version None\n", - " None\n", - "XC functionals = custom\n", - "radial grids: \n", - "Treutler-Ahlrichs [JCP 102, 346 (1995); DOI:10.1063/1.469408] (M4) radial grids\n", - "\n", - "becke partition: Becke, JCP 88, 2547 (1988); DOI:10.1063/1.454033\n", - "pruning grids: \n", - "grids dens level: 3\n", - "symmetrized grids: False\n", - "atomic radii adjust function: \n", - "small_rho_cutoff = 1e-07\n", - "Set gradient conv threshold to 3.16228e-05\n", - "Initial guess from minao.\n", - "tot grids = 33704\n", - "init E= -76.2302482032812\n", - " HOMO = -0.430502394424911 LUMO = -0.00557962018270088\n", - "cycle= 1 E= -76.1844346685567 delta_E= 0.0458 |g|= 0.818 |ddm|= 1.41\n", - " HOMO = -0.0571497729893296 LUMO = 0.092111968345332\n", - "cycle= 2 E= -75.9921484773828 delta_E= 0.192 |g|= 1.28 |ddm|= 1.4\n", - " HOMO = -0.266575634970147 LUMO = 0.0399392427556757\n", - "cycle= 3 E= -76.3081751739264 delta_E= -0.316 |g|= 0.0337 |ddm|= 0.956\n", - " HOMO = -0.25652282325852 LUMO = 0.0420512491571007\n", - "cycle= 4 E= -76.30841420647 delta_E= -0.000239 |g|= 0.0064 |ddm|= 0.0228\n", - " HOMO = -0.257181682273231 LUMO = 0.0406822667985362\n", - "cycle= 5 E= -76.3084203471352 delta_E= -6.14e-06 |g|= 0.00122 |ddm|= 0.00337\n", - " HOMO = -0.257376647610882 LUMO = 0.0405720409497939\n", - "cycle= 6 E= -76.3084205799498 delta_E= -2.33e-07 |g|= 3.1e-05 |ddm|= 0.000715\n", - " HOMO = -0.257386810929883 LUMO = 0.0405644826339652\n", - "cycle= 7 E= -76.3084206025217 delta_E= -2.26e-08 |g|= 3.89e-06 |ddm|= 3.08e-05\n", - " HOMO = -0.257387527784119 LUMO = 0.0405648211980371\n", - "cycle= 8 E= -76.3084205960069 delta_E= 6.51e-09 |g|= 3.05e-07 |ddm|= 3.24e-06\n", - " HOMO = -0.25738749260452 LUMO = 0.0405648168862294\n", - "cycle= 9 E= -76.3084206012037 delta_E= -5.2e-09 |g|= 3.17e-08 |ddm|= 2.17e-07\n", - " HOMO = -0.257387492441066 LUMO = 0.0405648161651567\n", - "cycle= 10 E= -76.3084206024703 delta_E= -1.27e-09 |g|= 2.04e-08 |ddm|= 1.73e-08\n", - " HOMO = -0.257387491334681 LUMO = 0.0405648167913778\n", - "cycle= 11 E= -76.3084206001009 delta_E= 2.37e-09 |g|= 1.44e-08 |ddm|= 7.84e-09\n", - " HOMO = -0.257387491019448 LUMO = 0.0405648174191108\n", - "cycle= 12 E= -76.3084205990203 delta_E= 1.08e-09 |g|= 1.38e-08 |ddm|= 4.07e-09\n", - " HOMO = -0.257387490637195 LUMO = 0.0405648176151493\n", - "cycle= 13 E= -76.3084206006191 delta_E= -1.6e-09 |g|= 2.06e-08 |ddm|= 7.24e-09\n", - " HOMO = -0.257387491245932 LUMO = 0.0405648173875793\n", - "cycle= 14 E= -76.3084206022113 delta_E= -1.59e-09 |g|= 1.76e-08 |ddm|= 4.61e-09\n", - " HOMO = -0.257387491753683 LUMO = 0.0405648174842908\n", - "cycle= 15 E= -76.3084205999321 delta_E= 2.28e-09 |g|= 2.05e-08 |ddm|= 4.88e-09\n", - " HOMO = -0.25738749162096 LUMO = 0.0405648180109603\n", - "cycle= 16 E= -76.3084206010219 delta_E= -1.09e-09 |g|= 1.75e-08 |ddm|= 3.15e-09\n", - " HOMO = -0.257387492013501 LUMO = 0.0405648177311237\n", - "cycle= 17 E= -76.3084206001081 delta_E= 9.14e-10 |g|= 1.46e-08 |ddm|= 3.2e-09\n", - " HOMO = -0.257387493925225 LUMO = 0.0405648195794535\n", - "Extra cycle E= -76.3084206016336 delta_E= -1.53e-09 |g|= 1.48e-08 |ddm|= 1.74e-08\n", - "converged SCF energy = -76.3084206016336\n", - "Dipole moment(X, Y, Z, Debye): -0.00000, 0.00000, -2.00440\n", - "Total energy: -2076.457890 eV\n", - "Total energy: -76.308421 Hartree\n" - ] - } - ], + "outputs": [], "source": [ "# Calculate the total energy\n", "energy_eV = atoms.get_potential_energy() # eV\n", @@ -232,33 +99,15 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": null, "id": "b2cc79b6", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "{'with_density_fit': True, 'verbose': 0}\n", - "Changed parameters: {'with_density_fit': True, 'verbose': 0}\n", - "\n", - "Current calculator parameters:\n", - " xc: skala\n", - " basis: def2-svp\n", - " with_density_fit: True\n", - " with_newton: False\n", - " with_dftd3: True\n", - " charge: None\n", - " multiplicity: None\n", - " verbose: 0\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "# Update calculator settings\n", - "changed_params = atoms.calc.set(with_density_fit=True, verbose=0)\n", + "changed_params = atoms.calc.set(\n", + " with_density_fit=True, verbose=0, ks_config={\"conv_tol\": 1e-6}\n", + ")\n", "print(changed_params)\n", "print(f\"Changed parameters: {changed_params}\")\n", "\n", @@ -282,25 +131,10 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": null, "id": "be609051", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "/home/giulialuise/.cache/huggingface/hub/models--microsoft--skala/snapshots/9eedaca41deeaebe197bec77d1898f777770d3ef/skala-1.0.fun\n", - "Downloaded in 0.10 seconds\n", - "Forces on atoms (eV/Å):\n", - " Atom 1 (O): [ 0.0000, 0.0000, -0.5867]\n", - " Atom 2 (H): [ -0.0000, -0.6148, 0.2933]\n", - " Atom 3 (H): [ 0.0000, 0.6148, 0.2933]\n", - "\n", - "Maximum force component: 0.6148 eV/Å\n" - ] - } - ], + "outputs": [], "source": [ "# Calculate forces\n", "forces = atoms.get_forces() # eV/Å\n", @@ -328,27 +162,10 @@ }, { "cell_type": "code", - "execution_count": 7, + "execution_count": null, "id": "cb74929f", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - " Step Time Energy fmax\n", - "LBFGSLineSearch: 0 15:43:56 -2076.458638 0.681230\n", - "LBFGSLineSearch: 1 15:44:02 -2076.466242 0.308524\n", - "LBFGSLineSearch: 2 15:44:13 -2076.467743 0.153213\n", - "LBFGSLineSearch: 3 15:44:26 -2076.468485 0.002282\n", - "Optimization completed in 3 steps\n", - "Initial energy: -2076.458638 eV\n", - "Final energy: -2076.468485 eV\n", - "Energy change: -0.009847 eV\n", - "Maximum final force: 0.0023 eV/Å\n" - ] - } - ], + "outputs": [], "source": [ "# Store initial geometry for comparison\n", "initial_positions = atoms.positions.copy()\n", @@ -382,19 +199,10 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": null, "id": "192bdbfa", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Dipole moment vector (Debye): [ -0.0000, 0.0000, -0.4200]\n", - "Dipole moment magnitude: 0.4200 Debye\n" - ] - } - ], + "outputs": [], "source": [ "# Calculate dipole moment\n", "dipole = atoms.get_dipole_moment() # Debye\n", @@ -417,26 +225,10 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "id": "11a8bbd0", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Methane molecule:\n", - "Chemical formula: CH4\n", - "Number of atoms: 5\n", - "/home/giulialuise/.cache/huggingface/hub/models--microsoft--skala/snapshots/9eedaca41deeaebe197bec77d1898f777770d3ef/skala-1.0.fun\n", - "Downloaded in 0.11 seconds\n", - "\n", - "Calculation results:\n", - "Total energy: -1102.070344 eV\n", - "Maximum force: 0.351115 eV/Å\n" - ] - } - ], + "outputs": [], "source": [ "# Create a methane molecule\n", "atoms = molecule(\"CH4\")\n", @@ -497,7 +289,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.7" + "version": "3.13.5" } }, "nbformat": 4, diff --git a/src/skala/ase/calculator.py b/src/skala/ase/calculator.py index c9dee44..4b623de 100644 --- a/src/skala/ase/calculator.py +++ b/src/skala/ase/calculator.py @@ -16,6 +16,7 @@ from skala.functional.base import ExcFunctionalBase from skala.pyscf import SkalaKS +from skala.pyscf.retry import retry_scf class Skala(Calculator): # type: ignore[misc] @@ -43,9 +44,11 @@ class Skala(Calculator): # type: ignore[misc] "auxbasis": None, "with_newton": False, "with_dftd3": True, + "with_retry": True, "charge": None, "multiplicity": None, "verbose": 0, + "ks_config": None, } _mol: gto.Mole | None = None @@ -72,6 +75,10 @@ def set(self, **kwargs: Any) -> dict[str, Any]: self._ks.verbose = verbose self._ks.base.verbose = verbose + if "ks_config" in changed_parameters and self.parameters.ks_config is not None: + if self._ks is not None: + self._ks.base(**self.parameters.ks_config) + if ( "charge" in changed_parameters or "multiplicity" in changed_parameters @@ -165,12 +172,17 @@ def calculate( auxbasis=self.parameters.auxbasis, with_newton=bool(self.parameters.with_newton), with_dftd3=bool(self.parameters.with_dftd3), + ks_config=self.parameters.ks_config, ).nuc_grad_method() self._ks = grad_method else: self._ks.reset(self._mol) - energy = self._ks.base.kernel() + if self.parameters.with_retry: + self._ks.base, _ = retry_scf(self._ks.base) + energy = self._ks.base.e_tot + else: + energy = self._ks.base.kernel() gradient = self._ks.kernel() self.results["energy"] = float(energy) * Hartree diff --git a/tests/test_ase.py b/tests/test_ase.py index 78a09d5..a397c1c 100644 --- a/tests/test_ase.py +++ b/tests/test_ase.py @@ -53,3 +53,26 @@ def test_missing_basis() -> None: calculator.InputError, match="Basis set must be specified in the parameters." ): atoms.get_potential_energy() + + +@pytest.mark.skipif(ase is None, reason="ASE is not installed") +def test_ks_config() -> None: + atoms = molecule("H2O") + atoms.calc = Skala( + xc="pbe", + basis="def2-svp", + with_density_fit=True, + auxbasis="def2-svp-jkfit", + ks_config={"conv_tol": 1e-6}, + ) + + energy = atoms.get_potential_energy() + + assert ( + atoms.calc._ks.base.conv_tol == 1e-6 + ), "KS solver convergence tolerance not set correctly" + + reference_energy = -2075.4896490374904 + assert ( + pytest.approx(energy, rel=1e-3) == reference_energy + ), f"Energy mismatch with custom KS config: {energy} vs {reference_energy}"