Separate relaxation settings from hessian settings

notebooks/0_create-test-set/0_get-exact-answer.ipynb

@@ -23,7 +23,7 @@
     "from ase.vibrations import VibrationsData, Vibrations\n",
     "from ase.calculators.mopac import MOPAC\n",
     "from ase.calculators.psi4 import Psi4\n",
-    "from ase.optimize import QuasiNewton\n",
+    "from ase.optimize import BFGS\n",
     "from ase import Atoms, units\n",
     "from ase.io import write, read\n",
     "from jitterbug.utils import make_calculator\n",
@@ -58,7 +58,8 @@
    "outputs": [],
    "source": [
     "molecule_name = 'caffeine'\n",
-    "method = 'pm7'\n",
+    "relax_method = 'pm7/None'  # Method used to relax geometry \n",
+    "hess_method = None  # Method used to perform Hessian computation, None to use same\n",
     "basis = None  # Set to None for MOPAC methods\n",
     "threads = min(os.cpu_count(), 12)\n",
     "delta = 0.01"
@@ -72,18 +73,37 @@
     "Derived"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2f543fd2-ca4c-4d68-a523-14515f351c4b",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "relax_method, relax_basis = relax_method.split(\"/\")\n",
+    "if hess_method is None:\n",
+    "    hess_method, hess_basis = relax_method, relax_basis\n",
+    "else:\n",
+    "    hess_method, hess_basis = hess_method.split(\"/\")"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "aeebbc77-70e4-4709-90a0-b9aaf54d4cd9",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "outputs": [],
    "source": [
-    "run_name = f'{molecule_name}_{method}_{basis}'\n",
+    "run_name = f'{molecule_name}_{hess_method}_{hess_basis}_at_{relax_method}_{relax_basis}'\n",
     "run_name_with_delta = f'{run_name}_d={delta:.3g}'\n",
     "out_dir = Path('data') / 'exact'\n",
     "if (out_dir / f'{run_name_with_delta}-times.json').exists():\n",
-    "    raise ValueError('Already done!')"
+    "    raise ValueError('Already done!')\n",
+    "print(f'Run name: {run_name_with_delta}')"
    ]
   },
   {
@@ -160,7 +180,7 @@
    },
    "outputs": [],
    "source": [
-    "calc = make_calculator(method, basis, num_threads=threads)"
+    "calc = make_calculator(relax_method, relax_basis, num_threads=threads)"
    ]
   },
   {
@@ -180,14 +200,17 @@
    },
    "outputs": [],
    "source": [
-    "geom_path = out_dir / f'{run_name}.xyz'"
+    "geom_path = out_dir / f'{molecule_name}_{relax_method}_{relax_basis}.xyz'\n",
+    "print(f'Geometry path: {geom_path}')"
    ]
   },
   {
    "cell_type": "code",
    "execution_count": null,
    "id": "ef903a43-5d6c-47fb-a500-837599c95f91",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
    "outputs": [],
    "source": [
     "%%time\n",
@@ -196,7 +219,7 @@
     "    atoms.calc = calc\n",
     "else:\n",
     "    atoms.calc = calc\n",
-    "    dyn = QuasiNewton(atoms)\n",
+    "    dyn = BFGS(atoms)\n",
     "    with redirect_stderr(devnull):\n",
     "        dyn.run(fmax=0.01)"
    ]
@@ -242,6 +265,35 @@
     "ASE has a built-in method which uses finite displacements"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "0e70265b-cefd-4d3c-925b-8b2cf13419e4",
+   "metadata": {},
+   "source": [
+    "Make the calculator for the hessian"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a144434b-e478-42e0-a2bd-5c43beab31d0",
+   "metadata": {
+    "tags": []
+   },
+   "outputs": [],
+   "source": [
+    "calc = make_calculator(hess_method, hess_basis, num_threads=threads)\n",
+    "atoms.calc = calc"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "71949047-df5f-47a5-883f-c329b7ca12bf",
+   "metadata": {},
+   "source": [
+    "Perform the computation"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
@@ -332,11 +384,11 @@
    "outputs": [],
    "source": [
     "psi4_path = out_dir / f'{run_name}-psi4.json'\n",
-    "if isinstance(calc, Psi4) and \"cc\" not in method and not psi4_path.exists():\n",
+    "if isinstance(calc, Psi4) and \"cc\" not in hess_method and not psi4_path.exists():\n",
     "    # Compute\n",
     "    analytic_time = perf_counter()\n",
     "    calc.set_psi4(atoms)\n",
-    "    hess = calc.psi4.hessian(f'{method}/{basis}')\n",
+    "    hess = calc.psi4.hessian(f'{hess_method}/{hess_basis}')\n",
     "    analytic_time = perf_counter() - analytic_time\n",
     "\n",
     "    # Convert to ASE format\n",

notebooks/0_create-test-set/1_compare-step-sizes-and-methods.ipynb

@@ -52,7 +52,8 @@
    },
    "outputs": [],
    "source": [
-    "molecule_name = 'water'  # Which water molecule to evaluate\n",
+    "molecule_name = 'water'\n",
+    "relax_level = 'b3lyp_cc-pvtz'  # Which water molecule to evaluate. Need both the molecule name and relaxation level\n",
     "target_method = ('ccsd(t)', 'cc-pvtz')"
    ]
   },
@@ -74,7 +75,7 @@
    },
    "outputs": [],
    "source": [
-    "hessian_paths = list(Path('data/exact/').glob(f'{molecule_name}_*-ase.json'))\n",
+    "hessian_paths = list(Path('data/exact/').glob(f'{molecule_name}_*_at_{relax_level}_d=*-ase.json'))\n",
     "print(f'Found {len(hessian_paths)} hessians for {molecule_name}')"
    ]
   },
@@ -87,7 +88,7 @@
    },
    "outputs": [],
    "source": [
-    "exact_path = Path(f'data/exact/{molecule_name}_{\"_\".join(target_method)}_d=0.005-ase.json')\n",
+    "exact_path = Path(f'data/exact/{molecule_name}_{\"_\".join(target_method)}_at_{relax_level}_d=0.005-ase.json')\n",
     "assert exact_path.exists(), f'Missing reference calculation: {exact_path}'\n",
     "exact_hess = VibrationsData.read(exact_path)"
    ]
@@ -105,15 +106,15 @@
     "    \"\"\"Load the Hessian and parse the metadata from the filename\n",
     "    \n",
     "    Args:\n",
-    "        path: Path to the Hessia\n",
+    "        path: Path to the Hessian\n",
     "    Returns:\n",
     "        Dictionary the includes the metadata for the calculation and errors wrt true Hessian\n",
     "    \"\"\"\n",
     "    \n",
     "    # Get some of the basic information\n",
     "    method_name, delta = path.name[:-9].rsplit(\"_d=\", 1)\n",
     "    delta = float(delta)\n",
-    "    _, method, basis = method_name.split(\"_\")\n",
+    "    _, method, basis = method_name.split(\"_\")[:3]\n",
     "    \n",
     "    # Compare to reference\n",
     "    approx_hess = VibrationsData.read(path)\n",

notebooks/0_create-test-set/run-all-methods.sh

@@ -1,18 +1,17 @@
 #! /bin/bash
 
-molecule=butanol
-methods="pm7//None xtb//None hf//cc-pvtz b3lyp//cc-pvtz wb97x-d//cc-pvtz m062x//cc-pvtz ccsd(t)//cc-pvdz"
+molecule=water
+relax_method="b3lyp/cc-pvtz"
+
+hess_methods="pm7/None xtb/None hf/cc-pvtz b3lyp/cc-pvtz wb97x-d/cc-pvtz m062x/cc-pvtz ccsd(t)/cc-pvdz"
 deltas="0.04 0.02 0.01 0.005 0.0025"
 
-#methods="ccsd(t)//cc-pvtz"
+#hess_methods="ccsd(t)/cc-pvtz"
 #deltas=0.005
 
 notebook=0_get-exact-answer.ipynb
-for name in $methods; do
-    echo $name
+for method in $hess_methods; do
     for delta in $deltas; do
-        method=$(echo $name | cut -d "/" -f 1)
-        basis=$(echo $name | cut -d "/" -f 3)
-        papermill -p method $method -p basis $basis -p delta $delta -p molecule_name $molecule $notebook live.ipynb
+        papermill -p hess_method $method -p relax_method $relax_method -p delta $delta -p molecule_name $molecule $notebook live.ipynb
     done
 done

notebooks/1_explore-sampling-methods/0_random-directions-same-distance.ipynb

@@ -47,9 +47,8 @@
    },
    "outputs": [],
    "source": [
-    "molecule_name = 'caffeine'\n",
-    "method = 'hf'\n",
-    "basis = 'def2-svpd'\n",
+    "starting_geometry = '../data/exact/caffeine_pm7_None.xyz'\n",
+    "method = 'hf/def2-svpd'\n",
     "threads = min(os.cpu_count(), 12)\n",
     "step_size: float = 0.005 # Perturbation amount, used as maximum L2 norm"
    ]
@@ -69,8 +68,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "run_name = Path(starting_geometry).name[:-4]\n",
-    "name, method, basis = run_name.split(\"_\")"
+    "relax_name = Path(starting_geometry).name[:-4]\n",
+    "name, relax_method, relax_basis = relax_name.split(\"_\")\n",
+    "method, basis = method.split(\"/\")\n",
+    "run_name = f'{name}_{method}_{basis}_at_{relax_method}_{relax_basis}'\n",
+    "print(f'Run name: {run_name}')"
    ]
   },
   {

notebooks/1_explore-sampling-methods/1_random-directions-variable-distance.ipynb

@@ -39,7 +39,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "52252ee2-315c-48bb-8cba-07620e6e2faa",
+   "id": "c61794ce-ca24-470f-903f-4fc5118af1d3",
    "metadata": {
     "tags": [
      "parameters"
@@ -48,13 +48,14 @@
    "outputs": [],
    "source": [
     "starting_geometry = '../data/exact/caffeine_pm7_None.xyz'\n",
+    "method = 'pm7/None'\n",
     "threads = min(os.cpu_count(), 12)\n",
-    "step_size: float = 0.005 # Lambda parameter for an expontential distribution for the Perturbation amount"
+    "step_size: float = 0.005 # Perturbation amount, used as maximum L2 norm"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "7010df09-73b2-4d58-be03-15a5f0d04b4c",
+   "id": "7ebb8a2a-b2f8-4647-9cd4-d9a05efc4790",
    "metadata": {},
    "source": [
     "Derived"
@@ -63,12 +64,15 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0b6794cd-477f-45a1-b96f-2332804ddb20",
+   "id": "3177eaf6-7af7-4e7c-8440-e32c172f8669",
    "metadata": {},
    "outputs": [],
    "source": [
-    "run_name = Path(starting_geometry).name[:-4]\n",
-    "name, method, basis = run_name.split(\"_\")"
+    "relax_name = Path(starting_geometry).name[:-4]\n",
+    "name, relax_method, relax_basis = relax_name.split(\"_\")\n",
+    "method, basis = method.split(\"/\")\n",
+    "run_name = f'{name}_{method}_{basis}_at_{relax_method}_{relax_basis}'\n",
+    "print(f'Run name: {run_name}')"
    ]
   },
   {
@@ -226,7 +230,8 @@
     "    # Sample a perturbation\n",
     "    disp = np.random.normal(0, 1, size=(n_atoms * 3))\n",
     "    disp /= np.linalg.norm(disp)\n",
-    "    my_step_dist = np.random.exponential(scale=step_size)\n",
+    "    my_step_dist = np.random.uniform(0, step_size)\n",
+    "    pbar.set_description(f'd={my_step_dist:.3e}')\n",
     "    disp *= my_step_dist * len(atoms)\n",
     "    disp = disp.reshape((-1, 3))\n",
     "    \n",

notebooks/1_explore-sampling-methods/2_displace-along-axes.ipynb

@@ -52,14 +52,15 @@
    "outputs": [],
    "source": [
     "starting_geometry = '../data/exact/caffeine_pm7_None.xyz'\n",
+    "method = 'hf/def2-svpd'\n",
     "threads = min(os.cpu_count(), 12)\n",
     "step_size: float = 0.005 # Lambda parameter for an expontential distribution for the Perturbation amount\n",
     "perturbs_per_evaluation: int = 2  # Number of perturbations to perform at once"
    ]
   },
   {
    "cell_type": "markdown",
-   "id": "7010df09-73b2-4d58-be03-15a5f0d04b4c",
+   "id": "134b0aa4-f7ef-415f-8334-7039bdf66152",
    "metadata": {},
    "source": [
     "Derived"
@@ -68,12 +69,15 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0b6794cd-477f-45a1-b96f-2332804ddb20",
+   "id": "f47df53a-1b81-4504-a9db-2fcc583d7096",
    "metadata": {},
    "outputs": [],
    "source": [
-    "run_name = Path(starting_geometry).name[:-4]\n",
-    "name, method, basis = run_name.split(\"_\")"
+    "relax_name = Path(starting_geometry).name[:-4]\n",
+    "name, relax_method, relax_basis = relax_name.split(\"_\")\n",
+    "method, basis = method.split(\"/\")\n",
+    "run_name = f'{name}_{method}_{basis}_at_{relax_method}_{relax_basis}'\n",
+    "print(f'Run name: {run_name}')"
    ]
   },
   {

notebooks/1_explore-sampling-methods/3_displace-along-vibrational-modes.ipynb

@@ -29,6 +29,7 @@
     "from pathlib import Path\n",
     "from tqdm import tqdm \n",
     "import numpy as np\n",
+    "import shutil\n",
     "import os"
    ]
   },
@@ -52,6 +53,7 @@
    "outputs": [],
    "source": [
     "starting_geometry = '../data/exact/caffeine_pm7_None.xyz'\n",
+    "method = 'pm7/None'\n",
     "threads = min(os.cpu_count(), 12)\n",
     "step_size: float = 0.002 # Target energy increase (units: eV)\n",
     "perturbs_per_evaluation: int = 16  # Number of perturbations to perform at once\n",
@@ -70,14 +72,17 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "0b6794cd-477f-45a1-b96f-2332804ddb20",
+   "id": "91cc7cb8-a620-4395-84fc-533c041c652e",
    "metadata": {
     "tags": []
    },
    "outputs": [],
    "source": [
-    "run_name = Path(starting_geometry).name[:-4]\n",
-    "name, method, basis = run_name.split(\"_\")"
+    "relax_name = Path(starting_geometry).name[:-4]\n",
+    "name, relax_method, relax_basis = relax_name.split(\"_\")\n",
+    "method, basis = method.split(\"/\")\n",
+    "run_name = f'{name}_{method}_{basis}_at_{relax_method}_{relax_basis}'\n",
+    "print(f'Run name: {run_name}')"
    ]
   },
   {
@@ -142,6 +147,8 @@
    "source": [
     "%%time\n",
     "atoms.calc = lower_calc\n",
+    "if Path('vib').exists():\n",
+    "    shutil.rmtree('vib')\n",
     "vib = Vibrations(atoms)\n",
     "vib.run()"
    ]