vp_h4_square.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <center> Variational Post-selection\n",
    "## <center> H4 (square)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since `PyScf` is not portable in Windows, we have to calculate the following codebox in `Linux`. The MolecularData file will be saved in `filename.hdf5`. Use command `cp filename.hdf5 /mnt/d` in WSL to copy it to your Windows `D:`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# in Linux\n",
    "from openfermion.chem import MolecularData\n",
    "from openfermionpyscf import run_pyscf\n",
    "\n",
    "d = 1.7\n",
    "geometry = [('H', (0, 0, 0)), ('H', (0, 0, d)), ('H', (0, d, d)), ('H', (0, d, 0))]\n",
    "basis = \"sto-3g\"\n",
    "multiplicity = 1\n",
    "description = \"h4_square\"\n",
    "molecule = MolecularData(geometry, basis, multiplicity, description=description, filename=\"h4_square_17\")\n",
    "molecule = run_pyscf(molecule, run_mp2=True, run_cisd=True, run_ccsd=True, run_fci=True)\n",
    "molecule.save()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can further process in Windows. (Don't forget to move the MolecularData file to the file where your `.ipynb` is.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import tensorflow as tf\n",
    "import tensorcircuit as tc\n",
    "import cotengra\n",
    "import quimb\n",
    "\n",
    "optc = cotengra.ReusableHyperOptimizer(\n",
    "    methods=[\"greedy\"],\n",
    "    parallel=\"ray\",\n",
    "    minimize=\"combo\",\n",
    "    max_time=30,\n",
    "    max_repeats=1024,\n",
    "    progbar=True,\n",
    ")\n",
    "tc.set_contractor(\"custom\", optimizer=optc, preprocessing=True)\n",
    "\n",
    "K = tc.set_backend(\"tensorflow\")\n",
    "tc.set_dtype(\"complex128\")\n",
    "\n",
    "from openfermion.chem import MolecularData\n",
    "from openfermion.linalg import LinearQubitOperator\n",
    "from openfermion.transforms import (\n",
    "    get_fermion_operator,\n",
    "    jordan_wigner\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Molecule"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# check the informations below\n",
    "d = 1.7\n",
    "geometry = [('H', (0, 0, 0)), ('H', (0, 0, d)), ('H', (0, d, d)), ('H', (0, d, 0))]\n",
    "basis = \"sto-3g\"\n",
    "multiplicity = 1\n",
    "description = \"h4_square\"\n",
    "molecule = MolecularData(geometry, basis, multiplicity, description=description, filename=\"models/h4_square_17\")   # check the filename\n",
    "molecule.load()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get Fermionic Hamiltonian\n",
    "mh = molecule.get_molecular_hamiltonian()\n",
    "fh = get_fermion_operator(mh)\n",
    "\n",
    "# Transform into qubit Hamiltonian\n",
    "b = jordan_wigner(fh) \n",
    "ns = LinearQubitOperator(b).n_qubits\n",
    "print(\"system qubit number:\", ns)\n",
    "\n",
    "# Transform from OpenFermion to TensorCircuit\n",
    "lsb, wb = tc.templates.chems.get_ps(b, ns)\n",
    "mb = tc.quantum.PauliStringSum2COO_numpy(lsb, wb)\n",
    "mbd = mb.todense()\n",
    "mbd_tf = tc.array_to_tensor(mbd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e_anal = np.min(quimb.eigvalsh(mbd))\n",
    "print(\"fci energy: {:.10f} \\nccsd energy: {:.10f} \\nhf energy: {:.10f} \\nanalytical energy: {:.10f} \\n\".format(molecule.fci_energy, molecule.ccsd_energy, molecule.hf_energy, e_anal))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Conventional VQE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def U_conventional(ns, depth, params):\n",
    "    params = K.cast(params, \"complex128\")\n",
    "    c = tc.Circuit(ns)\n",
    "    idx = 0\n",
    "    \n",
    "    c.x(0)\n",
    "    c.x(1)\n",
    "    c.x(2)\n",
    "    c.x(3)\n",
    "\n",
    "    for _ in range(depth):\n",
    "        for i in range(ns-1):\n",
    "            c.rzz(i, (i + 1) % ns, theta=params[idx+2*i])\n",
    "            c.iswap(i, (i + 1) % ns, theta=params[idx+2*i+1])\n",
    "        idx+=2*ns-2\n",
    "\n",
    "        c.iswap(ns-1, 0, theta=params[idx])\n",
    "        idx += 1\n",
    "\n",
    "        for i in range(ns):\n",
    "            c.rz(i, theta=params[idx+i])\n",
    "        idx+=ns\n",
    "\n",
    "    return c, idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# circuit visulization, optional\n",
    "ns_ = ns\n",
    "p = 1\n",
    "cirq, idx = U_conventional(ns_, p, np.zeros([1000]))\n",
    "print(\"The number of parameters is\", idx)\n",
    "cirq.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def train(ns, p, e_anal, maxiter=10000, lr=0.01, stddev=1.0):\n",
    "    _, idx = U_conventional(ns, p, np.zeros([1000]))\n",
    "    params = tf.Variable(\n",
    "        initial_value=tf.random.normal(\n",
    "            shape=[idx], stddev=stddev, dtype=getattr(tf, tc.rdtypestr)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    exp_lr = tf.keras.optimizers.schedules.ExponentialDecay(\n",
    "        initial_learning_rate=lr, decay_steps=150, decay_rate=0.7\n",
    "    )\n",
    "    opt = tf.keras.optimizers.Adam(exp_lr)\n",
    "\n",
    "    e_list = []\n",
    "    re_list = []\n",
    "    params_list = []\n",
    "    for i in range(maxiter):\n",
    "        with tf.GradientTape() as tape:\n",
    "            c, _ = U_conventional(ns, p, params)\n",
    "            e = tc.templates.measurements.operator_expectation(c, mbd_tf)\n",
    "        e_list.append(K.real(e).numpy())\n",
    "        re_list.append(np.abs((e_list[-1] - e_anal) / e_anal))\n",
    "        params_list.append(params.numpy())\n",
    "        grads = tape.gradient(e, params)\n",
    "        opt.apply_gradients(zip([grads], [params]))\n",
    "        if (i + 1) % 50 == 0 or i==0:\n",
    "            print(\"epoch{:>4}, e: {:.10f}, re: {:.10f}\".format(i, e_list[-1], re_list[-1]))\n",
    "    \n",
    "    params_list.append(params.numpy())\n",
    "    print(params.numpy())\n",
    "    \n",
    "    return e_list, re_list, params_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "\n",
    "p = 8\n",
    "maxiter = 500\n",
    "lr = 0.1\n",
    "stddev = np.pi\n",
    "\n",
    "t = time.gmtime()\n",
    "\n",
    "# train\n",
    "with tf.device(\"/cpu:0\"):\n",
    "    with open('data/h4_square_17_p{}_{}.csv'.format(p, time.strftime(\"%Y_%m%d_%H_%M_%S\",t)), 'w', newline='') as fp1:   # check the filename\n",
    "        with open('data/h4_square_17_p{}_params_{}.csv'.format(p, time.strftime(\"%Y_%m%d_%H_%M_%S\",t)), 'w', newline='') as fp2:   # check the filename\n",
    "            writer1 = csv.writer(fp1)\n",
    "            writer2 = csv.writer(fp2)\n",
    "            for j in range(50):\n",
    "                print(j+1)\n",
    "                e_list, re_list, params_list = train(ns, p, e_anal, maxiter=maxiter, lr=lr, stddev=stddev)\n",
    "                writer1.writerow([j])\n",
    "                writer1.writerow(e_list)\n",
    "                writer1.writerow(re_list)\n",
    "                writer2.writerow([j])\n",
    "                for k in range(maxiter+1):\n",
    "                    writer2.writerow(params_list[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <center> Post-selected VQE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def U_postselection(n, na, depth, params):\n",
    "    params = K.cast(params, \"complex128\")\n",
    "    c = tc.Circuit(n)\n",
    "    idx = 0\n",
    "    \n",
    "    c.x(1)\n",
    "    c.x(2)\n",
    "    c.x(3)\n",
    "    c.x(4)\n",
    "    \n",
    "    for _ in range(depth):\n",
    "        for i in range(n-1):\n",
    "            c.rzz(i, (i + 1) % n, theta=params[idx+2*i])\n",
    "            c.iswap(i, (i + 1) % n, theta=params[idx+2*i+1])\n",
    "        idx+=2*n-2\n",
    "\n",
    "        c.iswap(n-1, 0, theta=params[idx])\n",
    "        idx += 1\n",
    "\n",
    "        for i in range(n):\n",
    "            c.rz(i, theta=params[idx+i])\n",
    "        idx+=n\n",
    "\n",
    "    return c, idx"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# circuit visulization, optional\n",
    "ns_ = ns\n",
    "na = 1\n",
    "n = ns_ + na\n",
    "p = 1\n",
    "cirq, idx = U_postselection(n, na, p, np.zeros([1000]))\n",
    "print(\"The number of parameters is\", idx)\n",
    "cirq.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def postselect(n, na, ns, p, params):\n",
    "    c, _ = U_postselection(n, na, p, params)\n",
    "\n",
    "    w = c.wavefunction()[ : 2**ns]\n",
    "    norm = tf.linalg.norm(w)\n",
    "    w = w / norm\n",
    "    \n",
    "    return tc.Circuit(ns, inputs=w), tf.abs(norm)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only energy as loss\n",
    "def train(n, ns, na, p, e_anal, maxiter=10000, lr=0.01, stddev=1.0):\n",
    "    _, idx = U_postselection(n, na, p, np.zeros([1000]))\n",
    "    params = tf.Variable(\n",
    "        initial_value=tf.random.normal(\n",
    "            shape=[idx], stddev=stddev, dtype=getattr(tf, tc.rdtypestr)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    exp_lr = tf.keras.optimizers.schedules.ExponentialDecay(\n",
    "        initial_learning_rate=lr, decay_steps=150, decay_rate=0.7\n",
    "    )\n",
    "    opt = tf.keras.optimizers.Adam(exp_lr)\n",
    "\n",
    "    e_list = []\n",
    "    re_list = []\n",
    "    params_list = []\n",
    "    ps_list = []\n",
    "    for i in range(maxiter):\n",
    "        with tf.GradientTape() as tape:\n",
    "            c_s, norm = postselect(n, na, ns, p, params)\n",
    "            e = tc.templates.measurements.operator_expectation(c_s, mbd_tf)\n",
    "            ps_list.append(norm.numpy())\n",
    "        e_list.append(K.real(e).numpy())\n",
    "        re_list.append(np.abs((e_list[-1] - e_anal) / e_anal))\n",
    "        params_list.append(params.numpy())\n",
    "        grads = tape.gradient(e, params)\n",
    "        opt.apply_gradients(zip([grads], [params]))\n",
    "        if (i + 1) % 50 == 0 or i==0:\n",
    "            print(\"epoch{:>4}, e: {:.10f}, re: {:.10f}\".format(i, e_list[-1], re_list[-1]))\n",
    "    \n",
    "    params_list.append(params.numpy())\n",
    "    print(params.numpy())\n",
    "\n",
    "    return e_list, re_list, params_list, ps_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# consider ancilla accuracy in loss\n",
    "def train(n, ns, na, p, e_anal, maxiter=10000, lr=0.01, stddev=1.0):\n",
    "    _, idx = U_postselection(n, na, p, np.zeros([1000]))\n",
    "    params = tf.Variable(\n",
    "        initial_value=tf.random.normal(\n",
    "            shape=[idx], stddev=stddev, dtype=getattr(tf, tc.rdtypestr)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    exp_lr = tf.keras.optimizers.schedules.ExponentialDecay(\n",
    "        initial_learning_rate=lr, decay_steps=150, decay_rate=0.7\n",
    "    )\n",
    "    opt = tf.keras.optimizers.Adam(exp_lr)\n",
    "\n",
    "    e_list = []\n",
    "    re_list = []\n",
    "    params_list = []\n",
    "    ps_list = []\n",
    "    for i in range(maxiter):\n",
    "        with tf.GradientTape() as tape:\n",
    "            c_s, norm = postselect(n, na, ns, p, params)\n",
    "            e = tc.templates.measurements.operator_expectation(c_s, mbd_tf)\n",
    "            loss = e - 0.01 * norm\n",
    "            tape.watch(loss)\n",
    "    \n",
    "        e_list.append(K.real(e).numpy())\n",
    "        re_list.append(np.abs((e_list[-1] - e_anal) / e_anal))\n",
    "        params_list.append(params.numpy())\n",
    "        ps_list.append(norm.numpy())\n",
    "        if (i + 1) % 50 == 0 or i==0:\n",
    "            print(\"epoch{:>4}, e: {:.10f}, re: {:.10f}, ps: {:.2f}\".format(i, e_list[-1], re_list[-1], ps_list[-1]))\n",
    " \n",
    "        grads = tape.gradient(loss, params)\n",
    "        opt.apply_gradients(zip([grads], [params]))\n",
    "        \n",
    "    params_list.append(params.numpy())\n",
    "    print(params.numpy())\n",
    "\n",
    "    return e_list, re_list, params_list, ps_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# consider ancilla accuracy with sigmoid function in loss\n",
    "def train(n, ns, na, p, e_anal, maxiter=10000, lr=0.01, stddev=1.0):\n",
    "    _, idx = U_postselection(n, na, p, np.zeros([1000]))\n",
    "    params = tf.Variable(\n",
    "        initial_value=tf.random.normal(\n",
    "            shape=[idx], stddev=stddev, dtype=getattr(tf, tc.rdtypestr)\n",
    "        )\n",
    "    )\n",
    "\n",
    "    exp_lr = tf.keras.optimizers.schedules.ExponentialDecay(\n",
    "        initial_learning_rate=lr, decay_steps=150, decay_rate=0.7\n",
    "    )\n",
    "    opt = tf.keras.optimizers.Adam(exp_lr)\n",
    "\n",
    "    e_list = []\n",
    "    re_list = []\n",
    "    params_list = []\n",
    "    ps_list = []\n",
    "    for i in range(maxiter):\n",
    "        with tf.GradientTape() as tape:\n",
    "            c_s, norm = postselect(n, na, ns, p, params)\n",
    "            e = tc.templates.measurements.operator_expectation(c_s, mbd_tf)\n",
    "            loss = e - 0.005 * tf.sigmoid((norm-0.70) * 30)\n",
    "            tape.watch(loss)\n",
    "    \n",
    "        e_list.append(K.real(e).numpy())\n",
    "        re_list.append(np.abs((e_list[-1] - e_anal) / e_anal))\n",
    "        params_list.append(params.numpy())\n",
    "        ps_list.append(norm.numpy())\n",
    "        if (i + 1) % 50 == 0 or i==0:\n",
    "            print(\"epoch{:>4}, e: {:.10f}, re: {:.10f}, ps: {:.2f}, loss: {:.10f}\".format(i, e_list[-1], re_list[-1], ps_list[-1], loss))\n",
    " \n",
    "        grads = tape.gradient(loss, params)\n",
    "        opt.apply_gradients(zip([grads], [params]))\n",
    "        \n",
    "    params_list.append(params.numpy())\n",
    "    print(params.numpy())\n",
    "\n",
    "    return e_list, re_list, params_list, ps_list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import csv\n",
    "na = 1\n",
    "n = ns + na\n",
    "p = 8\n",
    "maxiter = 500\n",
    "lr = 0.1\n",
    "stddev = np.pi\n",
    "\n",
    "t = time.gmtime()\n",
    "\n",
    "# train\n",
    "with tf.device(\"/cpu:0\"):\n",
    "    with open('data/h4_square_17_ps_p{}_sigmoid70_{}.csv'.format(p, time.strftime(\"%Y_%m%d_%H_%M_%S\",t)), 'w', newline='') as fp1:   # check the filename\n",
    "        with open('data/h4_square_17_ps_p{}_sigmoid70_params_{}.csv'.format(p, time.strftime(\"%Y_%m%d_%H_%M_%S\",t)), 'w', newline='') as fp2:   # check the filename\n",
    "            writer1 = csv.writer(fp1)\n",
    "            writer2 = csv.writer(fp2)\n",
    "            for j in range(50):\n",
    "                print(j+1)\n",
    "                e_list, re_list, params_list, ps_list = train(n, ns, na, p, e_anal, maxiter=maxiter, lr=lr, stddev=stddev)\n",
    "                writer1.writerow([j])\n",
    "                writer1.writerow(e_list)\n",
    "                writer1.writerow(re_list)\n",
    "                writer1.writerow(ps_list)\n",
    "\n",
    "                writer2.writerow([j])\n",
    "                writer2.writerows(params_list)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}