Add pennylane and vcvqe modules

.gitignore

@@ -1,2 +1,5 @@
 __pycache__/
-*.npy
+.pytest_cache
+*.npy
+*.egg-info
+structure.xyz

README.md

@@ -0,0 +1,6 @@
+# DD2367-Quantum-Chemistry
+This is the final project of Samuil Donchev and William Galvin for Autumn 2023 DD2367 Quantum Computing for Computer Scientisits at KTH Stockholm. 
+
+In this project, we explore the variational quantum eigensolver (VQE) and its applications to computational chemistry, particularly in solving the electron structure problem. 
+
+In `/paper` we present an overview of our findings, in `/pennylane` we use an off-the-shelf library and tutorial to demonstrate the VQE with different molecules, and in `/vcvqe` we present a "very classical variational quantum eigensolver" that gives a lower level, less practical window into the VQE for quantum chemistry. See their READMEs for more information and usage instructions.

data/h2.xyz

@@ -1,4 +1,4 @@
 2
 
-H	0.0000	0.0000	-0.66
-H	0.0000	0.0000	0.66
+H	0.0000	0.0000	-0.75
+H	0.0000	0.0000	0

data/h2o2.xyz

@@ -0,0 +1,4 @@
+O	0.0000	0.7375	-0.0528
+O	0.0000	-0.7375	-0.0528
+H	0.8190	0.8170	0.4220
+H	-0.8190	-0.8170	0.4220

data/n2h4.xyz

@@ -0,0 +1,6 @@
+N	0.0000	0.7230	-0.1123
+N	0.0000	-0.7230	-0.1123
+H	-0.4470	1.0031	0.7562
+H	0.4470	-1.0031	0.7562
+H	0.9663	1.0031	0.0301
+H	-0.9663	-1.0031	0.0301

pennylane/README.md

@@ -0,0 +1,12 @@
+# PennyLane
+
+This is a slight modifaction of the code found in [this PennyLane tutorial](https://pennylane.ai/qml/demos/tutorial_vqe/).
+
+## Installation
+1. Clone this repo and `cd` into `pennylane`
+2. Run `pip install -r requirements.txt && pip install -e .`
+
+## Usage
+```
+vqe-pennylane --file ../data/[file].xyz [--max-iter <n>] [--threshold <n>]
+```

pennylane/requirements.txt

@@ -0,0 +1,2 @@
+pennylane>=0.32
+matplotlib

pennylane/setup.py

@@ -0,0 +1,11 @@
+from setuptools import setup
+
+setup(
+    name="vqe-pennylane",
+    version="0.1.0",
+    packages=["src"],
+    entry_points={
+        "console_scripts":
+        ["vqe-pennylane = src.main:main"]
+    },
+)

pennylane/src/main.py

@@ -0,0 +1,110 @@
+#!/usr/bin/env python3
+
+import argparse
+
+from pennylane import numpy as np
+import pennylane as qml
+
+
+def main():
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "--file",
+        required=True,
+        help="Path to .xyz file for input molecule."
+    )
+
+    parser.add_argument(
+        "--max-iter",
+        type=int, 
+        default=50,
+        help="Maximum number of optimization cycles"
+    )
+
+    parser.add_argument(
+        "--threshold",
+        type=float,
+        default=1e-6,
+        help="Convergence threshdold"
+    )
+
+    args = parser.parse_args()
+
+    symbols, geometry = qml.qchem.read_structure(args.file)
+
+    mol = qml.qchem.Molecule(symbols, geometry)
+    electrons = mol.n_electrons
+    H, qubits = qml.qchem.molecular_hamiltonian(symbols, geometry)
+
+    if qubits > 8:
+        print(f"{qubits} qubits: This may take awhile...")
+    else:
+        print("Number of qubits = ", qubits)
+
+
+    dev = qml.device("lightning.qubit", wires=qubits)
+
+    hf = qml.qchem.hf_state(electrons, qubits)
+    print("Hartree-fock basis state:", hf)
+
+
+    def circuit(param, wires):
+        qml.BasisState(hf, wires=wires)
+        singles, doubles = qml.qchem.excitations(electrons, qubits)
+
+        for i, s in enumerate(singles):
+            qml.SingleExcitation(param[i], wires=s)
+        for i, s in enumerate(doubles):
+            qml.DoubleExcitation(param[i + len(singles)], wires=s)
+
+
+    @qml.qnode(dev, interface="autograd")
+    def cost_fn(param):
+        circuit(param, wires=range(qubits))
+        return qml.expval(H)
+
+    opt = qml.GradientDescentOptimizer(stepsize=0.4)
+    singles, doubles = qml.qchem.excitations(electrons, qubits)
+    theta = np.random.random((len(singles) + len(doubles)), requires_grad=True)
+
+
+    energy = [cost_fn(theta)]
+    angle = [theta]
+
+    for n in range(args.max_iter):
+        theta, prev_energy = opt.step_and_cost(cost_fn, theta)
+
+        energy.append(cost_fn(theta))
+        angle.append(theta)
+
+        conv = np.abs(energy[-1] - prev_energy)
+
+        if n % 2 == 0:
+            print(f"Step = {n},  Energy = {energy[-1]:.8f} Ha")
+
+        if conv <= args.threshold:
+            break
+
+    print("\n" f"Final value of the ground-state energy = {energy[-1]:.8f} Ha")
+    print("\n" f"Optimal value of the circuit parameter:")
+    for i in theta:
+        print(f"{i: 5f}")
+
+
+    import matplotlib.pyplot as plt
+
+    plt.plot(range(n + 2), energy, "tab:red", ls="-.")
+    plt.xlabel("Optimization step", fontsize=13)
+    plt.ylabel("Energy (Hartree)", fontsize=13)
+    plt.ylabel(r"$E_\mathrm{HF}$", fontsize=15)
+    plt.title("Calculated Energy vs Optimization step")
+    plt.xticks(fontsize=12)
+    plt.yticks(fontsize=12)
+
+    plt.savefig("plot.png")
+    print("\nPlot saved at plot.png")
+
+
+if __name__ == "__main__":
+    main()

vcvqe/README.md

@@ -0,0 +1,39 @@
+# Very Classical Variational Quantum Eigensolver (VCVQE)
+
+## Background
+
+This is a classical implementation of the Variational Quantum Eigensolver for quantum chemistry described in [this PennyLane tutorial](https://pennylane.ai/qml/demos/tutorial_vqe/). See the accompanying write-up in `../paper` for more details. 
+
+This is considered a "classical" implementation because it takes the ideas from the VQE and implements them in classical matrix algebra, an approach that quickly becomes unmanageable as molecules grow large. For example, with a H2O molecule, this requires 14 qubits and therefor a Hamiltonian $H \in \mathbb{C}^{2^{14} \times 2^{14}}$, which is too large to even write down on a laptop. Nonetheless, it provides an instructive window into the inner workings of the VQE
+
+The energies calculated by this program account for only electron-electron and electron-nucleus interactions and electron kinetic energy, under the Born-Oppenheimer and Hartree-Fock approoximations. It does not incorporate nucleus-nucleus interactions, nuclear kinetic energies, or pairwise electron-electron energies. As such, **the energies calculated here do not match the energies calculatd by PySCF and PennyLane**.
+
+This program takes as input a molecular geometry (the 3D coordinates of the nuclei) and attempts to solve the electron structure problem. It prints a brief report of its findings.
+
+## Installation
+1. Clone this repo and `cd` to the root directory `vcvqe`
+2. Install with `pip install -r requirements.txt` and `pip install -e .`
+
+## Usage
+Run with 
+```
+$ vcvqe [-h] [--random] [--no-random] [--num-H NUM_H] [--file FILE] [--max-iter MAX_ITER] [--min-iter MIN_ITER] [--threshold THRESHOLD] [--step-size STEP_SIZE]
+```
+For example
+```
+vcvqe --file data/h4.xyz
+```
+or
+```
+vcvqe --random --num-H 2 --max-iter 25 --min-iter 5 --threshold 0.001 --step-size 0.25
+```
+
+This program can be used without `pip install -e .` by running `python src/main.py [args]`.
+
+### args
+- `random`: Use a hydrogen molecules with atoms randomly placed in 3D space. Almost certainly will not be a physically observable molecule. If `random` is used, `file` should not be used and `num-H` should be used.
+- `num-H`: Number of hydrogen atoms to use in random hydrogen molecule. H2 and H4 are supported, but note that H3 is not (PySCF doesn't support it). Defaults to 4
+- `file`: path to file with a molecule's `.xyz` coords. Theoretically may be any molecule, works best if it's a very small molecule (H2 or H4).
+- `[max|min]-iter` max/min iterations of optimization loop to perform. Defaults to 25 and 5
+- `threshold` If two consecutuve energies are calculated with a difference not greater than threshold, break optimization loop. Defaults to 0.001.
+- `step-size` Learning rate in optimization step. Defaults to 1

vcvqe/setup.py

@@ -0,0 +1,11 @@
+from setuptools import setup
+
+setup(
+    name="vcvqe",
+    version="0.1.0",
+    packages=["src"],
+    entry_points={
+        "console_scripts":
+        ["vcvqe = src.main:main"]
+    },
+)

src/main.py → vcvqe/src/main.py

@@ -11,14 +11,14 @@
 from jax import numpy as jnp
 from jax import value_and_grad
 
-from givens import givens
-from givens import hf_ground
-from givens import full_wavefunction
-from givens import excitations
-from hamiltonian import expectation
-from hamiltonian import hamiltonian
-from utils import get_electrons_qubits
-from utils import random_xyz
+from src.givens import givens
+from src.givens import hf_ground
+from src.givens import full_wavefunction
+from src.givens import excitations
+from src.hamiltonian import expectation
+from src.hamiltonian import hamiltonian
+from src.utils import get_electrons_qubits
+from src.utils import random_xyz
 
 
 def energy(H: np.ndarray, electrons: int, qubits: int, thetas: np.ndarray) -> float:
@@ -47,10 +47,15 @@ def get_args() -> argparse.Namespace:
 
     parser.add_argument(
         "--random",
-        action=argparse.BooleanOptionalAction,
+        action="store_true",
         default=False,
         help="Specify whether to use a random molecule or supply your own. If --random is set, --num-H must be too. If not, --file must be set"
     )
+    parser.add_argument(
+        "--no-random",
+        dest="random",
+        action="store_false"
+    )
 
     parser.add_argument(
         "--num-H",
@@ -92,13 +97,16 @@ def get_args() -> argparse.Namespace:
         help="Optimization step size aka learning rate"
     )
 
-    return parser.parse_args()
+    return parser.parse_args(), parser
 
 
-if __name__ == "__main__":
+def main():
     linebreak = "\n" + 100 * "=" + "\n"
 
-    args = get_args()
+    args, parser = get_args()
+    if not args.random and args.file is None:
+        parser.print_help()
+        exit(-1)
 
     with tempfile.NamedTemporaryFile() as f:
         if args.random:
@@ -116,11 +124,11 @@ def get_args() -> argparse.Namespace:
 
         try:
             electrons, qubits = get_electrons_qubits(file)
-        except RuntimeError:
-            raise RuntimeError("Molecule not valid for PySCF")
+        except RuntimeError as e:
+            raise RuntimeError(f"Molecule not valid for PySCF:\n {str(e)}")
 
+        print(f"\nCreating Hamiltonian for {electrons} electrons and {qubits} spin orbitals")
         H = hamiltonian(file)
-        print(f"\nHamiltonian created for {electrons} electrons and {qubits} spin orbitals")
         print(linebreak)
 
     singles, doubles = excitations(electrons, qubits)
@@ -142,11 +150,14 @@ def get_args() -> argparse.Namespace:
 
         print(f"\tstep: {i: 2}\t\tE: {E: .3f}\t\tDifference: {E - prev: .3f}")
 
-        if abs(E - prev) < args.threshold and i >= args.min_iter:
+        if abs(E - prev) < args.threshold and i >= args.min_iter - 1:
             break
 
         prev = E
 
     print(linebreak)
-    print(f"Final E converged to {E: 5f} in {i} steps\n")
-
+    print(f"Final E converged to {E: 5f} in {i + 1} steps\n")
+
+
+if __name__ == "__main__":
+    main()

tests/test_hamiltonian.py → vcvqe/tests/test_hamiltonian.py

@@ -1,7 +1,11 @@
 import numpy as np
+import os
 
 from src import hamiltonian
 
+root = "." if os.path.isdir("data") else ".."
+print(root)
+
 def test_creation_annihilation_shape():
     for n in range(8):
         for i in range(n):
@@ -22,6 +26,6 @@ def test_creation_annihilation_anti_symmetry():
 
 
 def test_hamiltonian_hermitian():
-    file = "data/h2.xyz"
+    file = f"{root}/data/h2.xyz"
     H = hamiltonian.hamiltonian(file)
     assert np.allclose(H, H.conj().T)

tests/test_utils.py → vcvqe/tests/test_utils.py

@@ -1,7 +1,10 @@
 import tempfile
+import os
 
 from src import utils
 
+root = "." if os.path.isdir("data") else ".."
+
 def test_random_xyz():
     atoms = {"H": 12, "O": 89, "Cl": 4}
     with tempfile.NamedTemporaryFile() as f:
@@ -12,11 +15,11 @@ def test_random_xyz():
 
 
 def test_get_electrons_qubits():
-    file = "data/h2.xyz"
+    file = f"{root}/data/h2.xyz"
     assert utils.get_electrons_qubits(file) == (2, 4)
 
-    file = "data/h4.xyz"
+    file = f"{root}/data/h4.xyz"
     assert utils.get_electrons_qubits(file) == (4, 8)
 
-    file = "data/h2o.xyz"
+    file = f"{root}/data/h2o.xyz"
     assert utils.get_electrons_qubits(file) == (10, 14)