added free projection without open ended walks

examples/13-free_projection/fp_afqmc.ipynb

@@ -0,0 +1,180 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from pyscf import fci, gto, scf\n",
+    "\n",
+    "np.set_printoptions(precision=5, suppress=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "r = 1.6\n",
+    "nH = 10\n",
+    "# geom = f\"H {-3*r/2} 0 0; H {-r/2} 0 0; H {r/2} 0 0; H {3*r/2} 0 0\"\n",
+    "geom = \"\"\n",
+    "for i in range(nH):\n",
+    "    geom += \"H 0 0 %g\\n\" % (i * r)\n",
+    "mol = gto.M(atom=geom, basis=\"sto-6g\", verbose=3, unit=\"bohr\")\n",
+    "\n",
+    "mf = scf.RHF(mol)\n",
+    "mf.kernel()\n",
+    "\n",
+    "umf = scf.UHF(mol)\n",
+    "umf.kernel()\n",
+    "mo1 = umf.stability(external=True)[0]\n",
+    "umf = umf.newton().run(mo1, umf.mo_occ)\n",
+    "mo1 = umf.stability(external=True)[0]\n",
+    "umf = umf.newton().run(mo1, umf.mo_occ)\n",
+    "\n",
+    "# fci\n",
+    "cisolver = fci.FCI(mol, mf.mo_coeff)\n",
+    "e, ci = cisolver.kernel()\n",
+    "print(\"FCI energy: \", e)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Free projection AFQMC\n",
+    "\n",
+    "Evaluates the quantity\n",
+    "\n",
+    "$$ E(\\tau) = \\frac{\\langle \\Psi_l | H e^{-\\tau H} | \\Psi_r \\rangle}{\\langle \\Psi_l | e^{-\\tau H} | \\Psi_r \\rangle} $$\n",
+    "\n",
+    "where $|\\Psi_l\\rangle$ is a trial wave function and $|\\Psi_r\\rangle$ is an initial state. The propagator is sampled using Monte Carlo. $E(\\tau)$ converges to the ground state energy at long $\\tau$, but the energies get noisier at long $\\tau$ due to the sign problem.\n",
+    "\n",
+    "In the following, energy evaluations are performed after a block consisting of `num_steps` steps of duration `dt`. In one iteration, energy samples are collected at `num_blocks` different $\\tau$ values. Multiple walkers are used to batch operations together for computational efficiency. The total number of samples at a given $\\tau$ is given by `num_walkers` $\\times$ `num_iterations_fp`. The energy is then averaged over walkers and iterations.\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ipie.addons.free_projection.qmc.calc import build_fpafqmc_driver\n",
+    "from ipie.config import MPI\n",
+    "from ipie.utils.from_pyscf import gen_ipie_input_from_pyscf_chk\n",
+    "\n",
+    "comm = MPI.COMM_WORLD\n",
+    "\n",
+    "gen_ipie_input_from_pyscf_chk(umf.chkfile, verbose=0)\n",
+    "qmc_options = {\n",
+    "    \"num_iterations_fp\": 100,\n",
+    "    \"num_blocks\": 4,\n",
+    "    \"num_steps\": 30,\n",
+    "    \"num_walkers\": 50,\n",
+    "    \"dt\": 0.05,\n",
+    "}\n",
+    "afqmc = build_fpafqmc_driver(\n",
+    "    comm,\n",
+    "    nelec=mol.nelec,\n",
+    "    seed=212503,\n",
+    "    qmc_options=qmc_options,\n",
+    ")\n",
+    "afqmc.run()\n",
+    "\n",
+    "# analysis\n",
+    "from ipie.addons.free_projection.analysis.extraction import extract_observable\n",
+    "from ipie.addons.free_projection.analysis.jackknife import jackknife_ratios\n",
+    "\n",
+    "for i in range(afqmc.params.num_blocks):\n",
+    "    print(\n",
+    "        f\"\\nEnergy statistics at time {(i+1) * afqmc.params.num_steps_per_block * afqmc.params.timestep}:\"\n",
+    "    )\n",
+    "    qmc_data = extract_observable(afqmc.estimators[i].filename, \"energy\")\n",
+    "    mean_energy, energy_err = jackknife_ratios(qmc_data[\"ENumer\"], qmc_data[\"EDenom\"])\n",
+    "    print(f\"  Energy: {mean_energy:.8e} +- {energy_err:.8e}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Phaseless AFQMC"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from ipie.qmc.calc import build_afqmc_driver\n",
+    "from ipie.config import MPI\n",
+    "from ipie.utils.from_pyscf import gen_ipie_input_from_pyscf_chk\n",
+    "\n",
+    "comm = MPI.COMM_WORLD\n",
+    "\n",
+    "gen_ipie_input_from_pyscf_chk(mf.chkfile, verbose=0)\n",
+    "\n",
+    "# fixing random seed for reproducibility\n",
+    "afqmc = build_afqmc_driver(comm, nelec=mol.nelec, num_walkers_per_task=100, seed=41100801)\n",
+    "if comm.rank == 0:\n",
+    "    print(afqmc.params)  # Inspect the default qmc options\n",
+    "\n",
+    "# Let us override the number of blocks to keep it short\n",
+    "afqmc.params.num_blocks = 400\n",
+    "afqmc.run()\n",
+    "\n",
+    "if comm.rank == 0:\n",
+    "    # We can extract the qmc data as as a pandas data frame like so\n",
+    "    from ipie.analysis.extraction import extract_observable\n",
+    "\n",
+    "    qmc_data = extract_observable(afqmc.estimators.filename, \"energy\")\n",
+    "    y = qmc_data[\"ETotal\"]\n",
+    "    y = y[50:]  # discard first 50 blocks\n",
+    "\n",
+    "    from ipie.analysis.autocorr import reblock_by_autocorr\n",
+    "\n",
+    "    df = reblock_by_autocorr(y, verbose=1)\n",
+    "    print(df.to_csv(index=False))\n",
+    "    # assert np.isclose(df.at[0,'ETotal_ac'], -5.325611614468466)\n",
+    "    # assert np.isclose(df.at[0,'ETotal_error_ac'], 0.00938082351500978)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "py39",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.5"
+  },
+  "orig_nbformat": 4
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}

examples/13-free_projection/run_afqmc.py

@@ -0,0 +1,54 @@
+from pyscf import cc, gto, scf
+
+from ipie.config import MPI
+from ipie.utils.from_pyscf import gen_ipie_input_from_pyscf_chk
+
+comm = MPI.COMM_WORLD
+mol = gto.M(
+    atom=[("H", 1.6 * i, 0, 0) for i in range(0, 10)],
+    basis="sto-6g",
+    verbose=4,
+    unit="Bohr",
+)
+if comm.rank == 0:
+    mf = scf.UHF(mol)
+    mf.chkfile = "scf.chk"
+    mf.kernel()
+    mycc = cc.UCCSD(mf).run()
+    et = mycc.ccsd_t()
+    print("UCCSD(T) energy {}".format(mf.e_tot + mycc.e_corr + et))
+
+    gen_ipie_input_from_pyscf_chk(mf.chkfile, verbose=0)
+comm.barrier()
+
+from ipie.addons.free_projection.qmc.calc import build_fpafqmc_driver
+
+qmc_options = {
+    "num_iterations_fp": 100,
+    "num_blocks": 5,
+    "num_steps": 20,
+    "num_walkers": 10,
+    "dt": 0.05,
+}
+afqmc = build_fpafqmc_driver(
+    comm,
+    nelec=mol.nelec,
+    seed=41100801,
+    qmc_options=qmc_options,
+)
+if comm.rank == 0:
+    print(afqmc.params)  # Inspect the default qmc options
+afqmc.run()
+
+# analysis
+if comm.rank == 0:
+    from ipie.addons.free_projection.analysis.extraction import extract_observable
+    from ipie.addons.free_projection.analysis.jackknife import jackknife_ratios
+
+    for i in range(afqmc.params.num_blocks):
+        print(
+            f"\nEnergy statistics at time {(i+1) * afqmc.params.num_steps_per_block * afqmc.params.timestep}:"
+        )
+        qmc_data = extract_observable(afqmc.estimators[i].filename, "energy")
+        energy_mean, energy_err = jackknife_ratios(qmc_data["ENumer"], qmc_data["EDenom"])
+        print(f"Energy: {energy_mean:.8e} +/- {energy_err:.8e}")

ipie/addons/free_projection/__init__.py

@@ -0,0 +1,13 @@
+# Copyright 2022 The ipie Developers. All Rights Reserved.
+#
+# 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.

ipie/addons/free_projection/analysis/__init__.py

@@ -0,0 +1,13 @@
+# Copyright 2022 The ipie Developers. All Rights Reserved.
+#
+# 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.

ipie/addons/free_projection/analysis/extraction.py

@@ -0,0 +1,65 @@
+# Copyright 2022 The ipie Developers. All Rights Reserved.
+#
+# 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.
+#
+# Authors: Fionn Malone <[email protected]>
+#          Joonho Lee <[email protected]>
+#
+
+
+import h5py
+import numpy
+import pandas as pd
+
+
+def extract_hdf5_data(filename, block_idx=1):
+    shapes = {}
+    with h5py.File(filename, "r") as fh5:
+        keys = fh5[f"block_size_{block_idx}/data/"].keys()
+        shape_keys = fh5[f"block_size_{block_idx}/shape/"].keys()
+        data = numpy.concatenate([fh5[f"block_size_{block_idx}/data/{d}"][:] for d in keys])
+        for k in shape_keys:
+            shapes[k] = {
+                "names": fh5[f"block_size_{block_idx}/names/{k}"][()],
+                "shape": fh5[f"block_size_{block_idx}/shape/{k}"][:],
+                "offset": fh5[f"block_size_{block_idx}/offset/{k}"][()],
+                "size": fh5[f"block_size_{block_idx}/size/{k}"][()],
+                "scalar": bool(fh5[f"block_size_{block_idx}/scalar/{k}"][()]),
+                "num_walker_props": fh5[f"block_size_{block_idx}/num_walker_props"][()],
+                "walker_header": fh5[f"block_size_{block_idx}/walker_prop_header"][()],
+            }
+        size_keys = fh5[f"block_size_{block_idx}/max_block"].keys()
+        max_block = sum(fh5[f"block_size_{block_idx}/max_block/{d}"][()] for d in size_keys)
+
+    return data[: max_block + 1], shapes
+
+
+def extract_observable(filename, name="energy", block_idx=1):
+    data, info = extract_hdf5_data(filename, block_idx=block_idx)
+    obs_info = info.get(name)
+    if obs_info is None:
+        raise RuntimeError(f"Unknown value for name={name}")
+    obs_slice = slice(obs_info["offset"], obs_info["offset"] + obs_info["size"])
+    if obs_info["scalar"]:
+        obs_data = data[:, obs_slice].reshape((-1,) + tuple(obs_info["shape"]))
+        nwalk_prop = obs_info["num_walker_props"]
+        weight_data = data[:, :nwalk_prop].reshape((-1, nwalk_prop))
+        results = pd.DataFrame(numpy.hstack([weight_data, obs_data]))
+        header = list(obs_info["walker_header"]) + obs_info["names"].split()
+        results.columns = [n.decode("utf-8") for n in header]
+        return results
+    else:
+        obs_data = data[:, obs_slice]
+        nsamp = data.shape[0]
+        walker_averaged = obs_data[:, :-1] / obs_data[:, -1].reshape((nsamp, -1))
+        return walker_averaged.reshape((nsamp,) + tuple(obs_info["shape"]))

ipie/addons/free_projection/analysis/jackknife.py

@@ -0,0 +1,52 @@
+# Copyright 2022 The ipie Developers. All Rights Reserved.
+#
+# 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.
+#
+# Authors: Fionn Malone <[email protected]>
+#          Joonho Lee
+#
+
+#!/usr/bin/env python
+
+import numpy
+
+
+def jackknife_ratios(num: numpy.ndarray, denom: numpy.ndarray):
+    r"""Jackknife estimation of standard deviation of the ratio of means.
+
+    Parameters
+    ----------
+    num : :class:`np.ndarray
+        Numerator samples.
+    denom : :class:`np.ndarray`
+        Denominator samples.
+
+    Returns
+    -------
+    mean : :class:`np.ndarray`
+        Ratio of means.
+    sigma : :class:`np.ndarray`
+        Standard deviation of the ratio of means.
+    """
+    n_samples = num.size
+    num_mean = numpy.mean(num)
+    denom_mean = numpy.mean(denom)
+    mean = num_mean / denom_mean
+    jackknife_estimates = numpy.zeros(n_samples, dtype=num.dtype)
+    for i in range(n_samples):
+        mean_num_i = (num_mean * n_samples - num[i]) / (n_samples - 1)
+        mean_denom_i = (denom_mean * n_samples - denom[i]) / (n_samples - 1)
+        jackknife_estimates[i] = (mean_num_i / mean_denom_i).real
+    mean = numpy.mean(jackknife_estimates)
+    sigma = numpy.sqrt((n_samples - 1) * numpy.var(jackknife_estimates))
+    return mean, sigma

ipie/addons/free_projection/analysis/tests/__init__.py

@@ -0,0 +1,13 @@
+# Copyright 2022 The ipie Developers. All Rights Reserved.
+#
+# 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.