Quantum_PSP2D_HPmodel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "3f6c669b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Enter a sequence of length at least 5 encoded as a string of 0s and 1s:11111\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Simulation failed and returned the following error message:\n",
      "ERROR:  [Experiment 0] Insufficient memory to run circuit circuit-88 using the statevector simulator. Required memory: 4096M, max memory: 3781M\n"
     ]
    },
    {
     "ename": "QiskitError",
     "evalue": "'ERROR:  [Experiment 0] Insufficient memory to run circuit circuit-88 using the statevector simulator. Required memory: 4096M, max memory: 3781M ,  ERROR: Insufficient memory to run circuit circuit-88 using the statevector simulator. Required memory: 4096M, max memory: 3781M'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mQiskitError\u001b[0m                               Traceback (most recent call last)",
      "Input \u001b[0;32mIn [1]\u001b[0m, in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m   1154\u001b[0m simulator \u001b[38;5;241m=\u001b[39m Aer\u001b[38;5;241m.\u001b[39mget_backend(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mqasm_simulator\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m   1155\u001b[0m result \u001b[38;5;241m=\u001b[39m execute(qc, backend \u001b[38;5;241m=\u001b[39m simulator, shots \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m10\u001b[39m)\u001b[38;5;241m.\u001b[39mresult()\n\u001b[0;32m-> 1156\u001b[0m counts \u001b[38;5;241m=\u001b[39m \u001b[43mresult\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget_counts\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m   1157\u001b[0m \u001b[38;5;66;03m#statevector = result.get_statevector(qc)\u001b[39;00m\n\u001b[1;32m   1158\u001b[0m \u001b[38;5;66;03m#print(statevector)\u001b[39;00m\n\u001b[1;32m   1159\u001b[0m \u001b[38;5;66;03m#from qiskit.tools.visualization import plot_histogram\u001b[39;00m\n\u001b[1;32m   1160\u001b[0m \u001b[38;5;66;03m#plot_histogram(counts)\u001b[39;00m\n\u001b[1;32m   1161\u001b[0m \u001b[38;5;28mprint\u001b[39m(counts, sep\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m'\u001b[39m)\n",
      "File \u001b[0;32m~/.local/lib/python3.8/site-packages/qiskit/result/result.py:269\u001b[0m, in \u001b[0;36mResult.get_counts\u001b[0;34m(self, experiment)\u001b[0m\n\u001b[1;32m    267\u001b[0m dict_list \u001b[38;5;241m=\u001b[39m []\n\u001b[1;32m    268\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m key \u001b[38;5;129;01min\u001b[39;00m exp_keys:\n\u001b[0;32m--> 269\u001b[0m     exp \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_get_experiment\u001b[49m\u001b[43m(\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    270\u001b[0m     \u001b[38;5;28;01mtry\u001b[39;00m:\n\u001b[1;32m    271\u001b[0m         header \u001b[38;5;241m=\u001b[39m exp\u001b[38;5;241m.\u001b[39mheader\u001b[38;5;241m.\u001b[39mto_dict()\n",
      "File \u001b[0;32m~/.local/lib/python3.8/site-packages/qiskit/result/result.py:397\u001b[0m, in \u001b[0;36mResult._get_experiment\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m    395\u001b[0m result_status \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mgetattr\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mstatus\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mResult was not successful\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m    396\u001b[0m exp_status \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mgetattr\u001b[39m(exp, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mstatus\u001b[39m\u001b[38;5;124m\"\u001b[39m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mExperiment was not successful\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m--> 397\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m QiskitError(result_status, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m, \u001b[39m\u001b[38;5;124m\"\u001b[39m, exp_status)\n",
      "\u001b[0;31mQiskitError\u001b[0m: 'ERROR:  [Experiment 0] Insufficient memory to run circuit circuit-88 using the statevector simulator. Required memory: 4096M, max memory: 3781M ,  ERROR: Insufficient memory to run circuit circuit-88 using the statevector simulator. Required memory: 4096M, max memory: 3781M'"
     ]
    }
   ],
   "source": [
    "'''\n",
    "(C) Copyright Renata Wong 2023.\n",
    "\n",
    "This code is licensed under the Apache License, Version 2.0. You may obtain a copy of this license \n",
    "at http://www.apache.org/licenses/LICENSE-2.0.\n",
    "\n",
    "Any modifications or derivative works of this code must retain this copyright notice, and modified files \n",
    "need to carry a notice indicating that they have been altered from the originals.\n",
    "\n",
    "This is the accompanying Qiskit code for the paper \n",
    "'Fast Quantum Algorithm for Protein Structure Prediction in Hydrophobic-Hydrophilic Model' (DOI: 10.1016/j.jpdc.2022.03.011)\n",
    "The code works for amino acid sequences of length > 4.\n",
    "\n",
    "NOTE: This code assumes the sequence 01001 for testing pursposes. Other sequences may require adjustments of variables\n",
    "num_solutions and j. \n",
    "'''\n",
    "\n",
    "from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "import numpy as np\n",
    "import math\n",
    "\n",
    "'''\n",
    "Taking in user input of the amino acid sequence.\n",
    "Bit 0 stands for a hydrophilic amino acid, while bit 1 stands for a hydrophobic amino acid\n",
    "'''\n",
    "\n",
    "sequence_input = list(input('Enter a sequence of length at least 5 encoded as a string of 0s and 1s:' ))\n",
    "length = len(sequence_input)\n",
    "\n",
    "'''\n",
    "Quantum register holding the sequence encoding\n",
    "'''\n",
    "sequence = QuantumRegister(length,'s')\n",
    "\n",
    "\n",
    "'''\n",
    "Quantum registers for each of the three coordinates x, y, and z. \n",
    "Coordinate order within each register: x0, x1, x2, etc., where x0 is the \n",
    "x-coordinate of the 1st amino acid, x1 is the x-coordinate of the 2nd amino acid, etc.\n",
    "Each register holds 2(length - 1) + 1 values.\n",
    "Each value can therefore be repserented with a minimum of ceil(log2(2*length - 1)) bits.\n",
    "There are a 'length' number of amino acids in a sequence.\n",
    "'''\n",
    "register_length = int(math.ceil(np.log2(2*length-1)))\n",
    "coord_reg_length = length * register_length\n",
    "\n",
    "x_coord = QuantumRegister(coord_reg_length, 'x')\n",
    "y_coord = QuantumRegister(coord_reg_length, 'y')\n",
    "\n",
    "'''\n",
    "Quantum register holding the controls w0, w1, etc.\n",
    "Each w is of length 2. \n",
    "Amino acid at position (0, 0, 0) is fixed. \n",
    "Hence, length - 1 amino acids in a sequence require directional information encoded in w. \n",
    "'''\n",
    "w = QuantumRegister(2*(length-1), 'w')\n",
    "\n",
    "'''\n",
    "Quantum register holding the binary 1, and the two's complement of 1 (= -1).\n",
    "The two's complement of 1 can be obtained by negating all the qubits in binary one but the last.\n",
    "NOTE: Encoding coordinates in two's complement notation allows one to treat subtraction as addition\n",
    "with the same ripple carry circuit. \n",
    "'''\n",
    "one_reg_length = int(math.ceil(np.log2(2*length-1)))\n",
    "binary_one = QuantumRegister(one_reg_length,'one')\n",
    "\n",
    "'''\n",
    "Quantum register holding the carry bit for ripple-carry adder\n",
    "'''\n",
    "carry = QuantumRegister(1, 'c')\n",
    "  \n",
    "    \n",
    "'''\n",
    "Register holding the Grover auxialliary qubit\n",
    "'''\n",
    "g = QuantumRegister(1, 'g')\n",
    "   \n",
    "    \n",
    "'''\n",
    "Classical register holding the output of the quantum computation\n",
    "'''\n",
    "output = ClassicalRegister(length+2*(length-1), 'out')#(length + 2*(length-1) + 2*coord_reg_length,'out')   \n",
    "\n",
    "\n",
    "'''\n",
    "The ancilla register of 3 qubits in needed for the multiplexer (controlled Toffoli gates) in the ripple-carry quantum adder.\n",
    "'''\n",
    "ancilla = QuantumRegister(3, 'anc')   \n",
    "\n",
    "\n",
    "'''\n",
    "Registers for the section FINDING ENERGY VALUES\n",
    "'''\n",
    "k_range = length-3\n",
    "j_range = int((k_range * (k_range+1))/2)\n",
    "    \n",
    "x_star = QuantumRegister(4*j_range*register_length, 'x_st')\n",
    "y_star = QuantumRegister(4*j_range*register_length, 'y_st')\n",
    "    \n",
    "r_reg = QuantumRegister(4*j_range*(2*register_length + 1), 'r')\n",
    "    \n",
    "delta = QuantumRegister(4*2*j_range, 'delta')\n",
    "    \n",
    "psi = QuantumRegister(4*j_range, 'psi')\n",
    "\n",
    "\n",
    "'''\n",
    "Registers for section SUMMING UP ENERGY VALUES. \n",
    "We create 4 index matrices that hold registers for each of 0 <= p <= 3 in the z_reg to obtain the value for kj in register z.\n",
    "'''\n",
    "\n",
    "p0_matrix = [ [ 0 for m in range(4 * (length-4) + 2) ] for n in range(length-2) ]\n",
    "\n",
    "kj = 0\n",
    "for k in range(1, length-2):\n",
    "    for j in range(4 * (k-1) + 2):\n",
    "        p0_matrix[k][j] = kj\n",
    "        kj += 1\n",
    "\n",
    "p1_matrix = [ [ 0 for m in range(4 * (length-4) + 3) ] for n in range(length-2) ]\n",
    "\n",
    "for k in range(1, length-2):\n",
    "    for j in range(4 * (k-1) + 3):\n",
    "        p1_matrix[k][j] = kj\n",
    "        kj += 1\n",
    "\n",
    "p2_matrix = [ [ 0 for m in range(4 * (length-4) + 4) ] for n in range(length-2) ]\n",
    "\n",
    "for k in range(1, length-2):\n",
    "    for j in range(4 * (k-1) + 4):\n",
    "        p2_matrix[k][j] = kj\n",
    "        kj += 1\n",
    "\n",
    "p3_matrix = [ [ 0 for m in range(4 * (length-4) + 5) ] for n in range(length-2) ]\n",
    "\n",
    "for k in range(1, length-2):\n",
    "    for j in range(4 * (k-1) + 5):\n",
    "        p3_matrix[k][j] = kj\n",
    "        kj += 1\n",
    "\n",
    "        \n",
    "        \n",
    "z = QuantumRegister(kj, 'z')      #((4 * length - 12) * (4 * length - 11) / 2, 'z') this value is too low  \n",
    "z_or = QuantumRegister(1, 'z_or')       # ancilla qubit for the multiplexer in Fig. 4   \n",
    "z_and = QuantumRegister(1, 'z_and')    # ancilla to hold the results of the AND operation\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "Defining a quantum circuit consisting of the above quantum and classical registers\n",
    "'''\n",
    "qc = QuantumCircuit(g, sequence, x_coord, y_coord, w, binary_one, carry, ancilla, x_star, y_star, r_reg, delta, psi, \n",
    "                    z, z_or, z_and, output)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "SECTION: PREPARING A UNIFORM SUPERPOSITION STATE\n",
    "\n",
    "Initialisation of sequence values based on user input\n",
    "'''\n",
    "for index, bit in enumerate(sequence_input):\n",
    "    if bit == '1':\n",
    "        qc.x(sequence[index])\n",
    "\n",
    "'''\n",
    "Initialising binary one and the two's complement of one\n",
    "''' \n",
    "qc.x(binary_one[one_reg_length-1])\n",
    "\n",
    "'''\n",
    "Setting the quantum state in a uniform superposition wrt. w\n",
    "'''\n",
    "qc.h(w[0:2*(length-1)])\n",
    "\n",
    "\n",
    "'''\n",
    "Initialization of r_reg\n",
    "'''\n",
    "multiplier = int(len(r_reg)/(4*j_range))\n",
    "for i in range(4*j_range):\n",
    "    qc.x(r_reg[i*multiplier])\n",
    "  \n",
    "\n",
    "'''\n",
    "Initialize z_or to |1> for cascading OR gates implemented with mcx gate.\n",
    "'''\n",
    "qc.x(z_or)\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "Setting the Grover qubit into the |-> state\n",
    "'''\n",
    "qc.x(g)\n",
    "qc.h(g)\n",
    "\n",
    "\n",
    "'''\n",
    "SECTION: CALCULATING COORDINATES\n",
    "Requires a doubly-controlled ripple-carry adder as a subroutine sc\n",
    "The adder below is for amino acid sequences of length > 4.\n",
    "'''\n",
    "summand_register_length = int(math.ceil(np.log2(2*length-1)))\n",
    "\n",
    "sw = QuantumRegister(2,'sw')     # control qubits (w)\n",
    "sa = QuantumRegister(3,'sa')     # ancilla qubits\n",
    "ss = QuantumRegister(1,'ss')     # carry \n",
    "sx = QuantumRegister(summand_register_length,'sx')   # summand x \n",
    "sy = QuantumRegister(summand_register_length,'sy')   # summand y\n",
    "sc = QuantumCircuit(sw, sa, ss, sx, sy, name='ccrca')\n",
    "\n",
    "sc.ccx(sw[0],sw[1],sa[0])\n",
    "\n",
    "for i in range(1, summand_register_length):\n",
    "    sc.ccx(sa[0],sx[i],sa[1])\n",
    "    sc.cx(sa[1],sy[i])\n",
    "    sc.ccx(sa[0],sx[i],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[1],sa[1])\n",
    "sc.cx(sa[1],ss[0])\n",
    "sc.ccx(sa[0],sx[1],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[0],sa[1])\n",
    "sc.ccx(sa[1],sy[0],sa[2])\n",
    "sc.cx(sa[2],ss[0])\n",
    "sc.ccx(sa[1],sy[0],sa[2])\n",
    "sc.ccx(sa[0],sx[0],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[2],sa[1])\n",
    "sc.cx(sa[1],sx[1])\n",
    "sc.ccx(sa[0],sx[2],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],ss[0],sa[1])\n",
    "sc.ccx(sa[1],sy[1],sa[2])\n",
    "sc.cx(sa[2],sx[1])\n",
    "sc.ccx(sa[1],sy[1],sa[2])\n",
    "sc.ccx(sa[0],ss[0],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[3],sa[1])\n",
    "sc.cx(sa[1],sx[2])\n",
    "sc.ccx(sa[0],sx[3],sa[1])\n",
    "\n",
    "for i in range(2, summand_register_length-2):\n",
    "    sc.ccx(sa[0],sy[i],sa[1])\n",
    "    sc.ccx(sa[1],sx[i-1],sa[2])\n",
    "    sc.cx(sa[2],sx[i])\n",
    "    sc.ccx(sa[1],sx[i-1],sa[2])\n",
    "    sc.ccx(sa[0],sy[i],sa[1])\n",
    "    sc.ccx(sa[0],sx[i+2],sa[1])\n",
    "    sc.cx(sa[1],sx[i+1])\n",
    "    sc.ccx(sa[0],sx[i+2],sa[1])\n",
    "    \n",
    "sc.ccx(sa[0],sy[summand_register_length-2],sa[1])\n",
    "sc.ccx(sa[1],sx[summand_register_length-3],sa[2])\n",
    "sc.cx(sa[2],sx[summand_register_length-2])\n",
    "sc.ccx(sa[1],sx[summand_register_length-3],sa[2])\n",
    "sc.ccx(sa[0],sy[summand_register_length-2],sa[1])\n",
    "\n",
    "for i in range(1, summand_register_length-1):\n",
    "    sc.cx(sa[0],sy[i])\n",
    "    \n",
    "sc.ccx(sa[0],ss[0],sa[1])\n",
    "sc.cx(sa[1],sy[1])\n",
    "sc.ccx(sa[0],ss[0],sa[1])\n",
    "\n",
    "for i in range(2, summand_register_length):\n",
    "    sc.ccx(sa[0],sx[i-1],sa[1])\n",
    "    sc.cx(sa[1],sy[i])\n",
    "    sc.ccx(sa[0],sx[i-1],sa[1])\n",
    "    \n",
    "sc.ccx(sa[0],sx[summand_register_length-3],sa[1])\n",
    "sc.ccx(sa[1],sy[summand_register_length-2],sa[2])\n",
    "sc.cx(sa[2],sx[summand_register_length-2])\n",
    "sc.ccx(sa[1],sy[summand_register_length-2],sa[2])\n",
    "sc.ccx(sa[0],sx[summand_register_length-3],sa[1])\n",
    "\n",
    "for i in range(summand_register_length-3, 1, -1):\n",
    "    sc.ccx(sa[0],sx[i-1],sa[1])\n",
    "    sc.ccx(sa[1],sy[i],sa[2])\n",
    "    sc.cx(sa[2],sx[i])\n",
    "    sc.ccx(sa[1],sy[i],sa[2])\n",
    "    sc.ccx(sa[0],sx[i-1],sa[1])\n",
    "    sc.ccx(sa[0],sx[i+2],sa[1])\n",
    "    sc.cx(sa[1],sx[i+1])\n",
    "    sc.ccx(sa[0],sx[i+2],sa[1])\n",
    "    sc.cx(sa[0],sy[i+1])\n",
    "\n",
    "sc.ccx(sa[0],ss[0],sa[1])\n",
    "sc.ccx(sa[1],sy[1],sa[2])\n",
    "sc.cx(sa[2],sx[1])\n",
    "sc.ccx(sa[1],sy[1],sa[2])\n",
    "sc.ccx(sa[0],ss[0],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[3],sa[1])\n",
    "sc.cx(sa[1],sx[2])\n",
    "sc.ccx(sa[0],sx[3],sa[1])\n",
    "\n",
    "sc.cx(sa[0],sy[2])\n",
    "\n",
    "sc.ccx(sa[0],sx[0],sa[1])\n",
    "sc.ccx(sa[1],sy[0],sa[2])\n",
    "sc.cx(sa[2],ss[0])\n",
    "sc.ccx(sa[1],sy[0],sa[2])\n",
    "sc.ccx(sa[0],sx[0],sa[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[2],sa[1])\n",
    "sc.cx(sa[1],sx[1])\n",
    "sc.ccx(sa[0],sx[2],sa[1])\n",
    "\n",
    "sc.cx(sa[0],sy[1])\n",
    "\n",
    "sc.ccx(sa[0],sx[1],sa[1])\n",
    "sc.cx(sa[1],ss[0])\n",
    "sc.ccx(sa[0],sx[1],sa[1])\n",
    "\n",
    "for i in range(summand_register_length):\n",
    "    sc.ccx(sa[0],sx[i],sa[1])\n",
    "    sc.cx(sa[1],sy[i])\n",
    "    sc.ccx(sa[0],sx[i],sa[1])\n",
    "    \n",
    "sc.ccx(sw[0],sw[1],sa[0])\n",
    "\n",
    "subinst = sc.to_instruction()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "Ripple-carry adder without control qubits |w>\n",
    "'''\n",
    "ssr = QuantumRegister(1,'ssr')     # carry \n",
    "sxr = QuantumRegister(summand_register_length,'sxr')   # summand x \n",
    "syr = QuantumRegister(summand_register_length,'syr')   # summand y\n",
    "rca = QuantumCircuit(ssr, sxr, syr, name='rca')\n",
    "\n",
    "for i in range(1, summand_register_length):\n",
    "    rca.cx(sxr[1],syr[i])\n",
    "\n",
    "rca.cx(sxr[1],ssr[0])\n",
    "rca.ccx(sxr[0],syr[0],ssr[0])\n",
    "rca.cx(sxr[2],sxr[1])\n",
    "rca.ccx(ssr[0],syr[1],sxr[1])\n",
    "rca.cx(sxr[3],sxr[2])\n",
    "\n",
    "for i in range(2, summand_register_length-2):\n",
    "    rca.ccx(syr[i],sxr[i-1],sxr[i])\n",
    "    rca.cx(sxr[i+2],sxr[i+1])\n",
    "    \n",
    "rca.ccx(syr[summand_register_length-2],sxr[summand_register_length-3],sxr[summand_register_length-2])\n",
    "\n",
    "for i in range(1, summand_register_length-1):\n",
    "    rca.x(syr[i])\n",
    "    \n",
    "rca.cx(ssr[0],syr[1])\n",
    "\n",
    "for i in range(2, summand_register_length):\n",
    "    rca.cx(sxr[i-1],syr[i])\n",
    "    \n",
    "rca.ccx(sxr[summand_register_length-3],syr[summand_register_length-2],sxr[summand_register_length-2])\n",
    "\n",
    "for i in range(summand_register_length-3, 1, -1):\n",
    "    rca.ccx(sxr[i-1],syr[i],sxr[i])\n",
    "    rca.cx(sxr[i+2],sxr[i+1])\n",
    "    rca.cx(sxr[i+2],syr[i+1])\n",
    "\n",
    "rca.ccx(ssr[0],syr[1],sxr[1])\n",
    "rca.cx(sxr[3],sxr[2])\n",
    "rca.x(syr[2])\n",
    "rca.ccx(sxr[0],syr[0],ssr[0])\n",
    "rca.cx(sxr[2],sxr[1])\n",
    "rca.x(syr[1])\n",
    "rca.cx(sxr[1],ssr[0])\n",
    "\n",
    "for i in range(summand_register_length):\n",
    "    rca.cx(sxr[i],syr[i])\n",
    "    \n",
    "subinst_rca = rca.to_instruction()\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "AMPLITUDE AMPLIFICATION PROCEDURE:\n",
    "\n",
    "Calculation of coordinates using the ripple-carry adder\n",
    "b is a global variable used to navigate among the values of vector w.\n",
    "d stands for the position of an amino acid in a sequence.\n",
    "\n",
    "MODIFY NUMBER OF SOLUTIONS num_solutions FOR DIFFERENT SEQUENCES.\n",
    "'''\n",
    "num_solutions = 8\n",
    "num_iter = int(math.ceil(np.pi*(np.sqrt(2**(2*(length-1)/num_solutions)))/4))    # Number of iterations \n",
    "\n",
    "\n",
    "def calculate_neighbour(register, b, d):\n",
    "    \n",
    "    arglist = []\n",
    "    \n",
    "    arglist.append(w[b])\n",
    "    arglist.append(w[b+1])\n",
    "    \n",
    "    for i in range(len(ancilla)):\n",
    "        arglist.append(ancilla[i])\n",
    "        \n",
    "    arglist.append(carry[0])\n",
    "\n",
    "    for i in range(summand_register_length-1, -1, -1): # digit order reversed for adder\n",
    "        arglist.append(binary_one[i])\n",
    "        \n",
    "    for i in range(d*summand_register_length-1, (d-1)*summand_register_length-1, -1): # digit order reversed for adder \n",
    "        arglist.append(register[i])\n",
    "    \n",
    "    return arglist\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "def update_coordinates(register_1, register_2):\n",
    "    \n",
    "    arglist = []\n",
    "        \n",
    "    arglist.append(carry[0])\n",
    "\n",
    "    for i in range(summand_register_length-1, -1, -1): # digit order reversed for adder\n",
    "        arglist.append(register_1[i])\n",
    "        \n",
    "    reg = len(register_2)\n",
    "    for i in reversed(range(reg)): # digit order reversed for adder \n",
    "        arglist.append(register_2[i])\n",
    "    \n",
    "    return arglist\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "THIS IS THE MAIN PART OF THE CODE. \n",
    "'''\n",
    "for grover_iteration in range(num_iter): \n",
    "    \n",
    "    b = 0\n",
    "\n",
    "    # Subroutine 1: Generating conformational space by calculating coordinates\n",
    "    # This part has been tested and performs correctly\n",
    "\n",
    "    for d in range(2, length+1):\n",
    "        for q in range(summand_register_length):\n",
    "            qc.cx(x_coord[(d-2)*summand_register_length+q],x_coord[(d-1)*summand_register_length+q])\n",
    "            qc.cx(y_coord[(d-2)*summand_register_length+q],y_coord[(d-1)*summand_register_length+q])\n",
    "          \n",
    "            \n",
    "        '''\n",
    "        Calculating the western neighbour of site d-1 (this is the case when |w> = |11>)\n",
    "        '''\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "        neighbour_west = calculate_neighbour(x_coord, b, d)\n",
    "        qc.append(subinst, neighbour_west)\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "\n",
    "        '''\n",
    "        Calculating the eastern neighbour of site d-1 (this is the case when |w> = |01>)\n",
    "        '''\n",
    "        qc.x(w[b])\n",
    "        neighbour_east = calculate_neighbour(x_coord, b, d)\n",
    "        qc.append(subinst, neighbour_east)\n",
    "        qc.x(w[b])\n",
    "        \n",
    "        '''\n",
    "        Calculating the northern neighbour of site d-1 (this is the case when |w> = |00>)\n",
    "        '''\n",
    "        qc.x(w[b])\n",
    "        qc.x(w[b+1])\n",
    "        neighbour_north = calculate_neighbour(y_coord, b, d)\n",
    "        qc.append(subinst, neighbour_north)\n",
    "        qc.x(w[b])\n",
    "        qc.x(w[b+1])\n",
    "        \n",
    "        '''\n",
    "        Calculating the southern neighbour of site d-1 (this is the case when |w> = |10>)\n",
    "        '''\n",
    "        qc.x(w[b+1])\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "        neighbour_south = calculate_neighbour(x_coord, b, d)\n",
    "        qc.append(subinst, neighbour_south)\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "        qc.x(w[b+1])\n",
    "\n",
    "        b = b+2\n",
    "        \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "    '''\n",
    "    SECTION: FINDING ENERGY VALUES\n",
    "    For the conformation that corresponds to the ground energy of the protein, the energy value will be e = 1, \n",
    "    while e = 0 for all other confornations. \n",
    "    The energy is calculated based on the number of hydrophobic-hydrophobic contacts in the lattice. \n",
    "    For each such contact, a 1 is added to the energy value for the given conformation.\n",
    "    \n",
    "    New registers for calculating energy values based on the coordinates x and y: |x_star> and |y_star>.\n",
    "    Register r_reg stores the result of coordinate comparison between sites k and j. \n",
    "    Register delta stores the information as to whether the k-th and j-th lattice sites are both hydrophobic and form a \n",
    "    loose contact (i.e., are adjacent in the lattice but not in the sequence). \n",
    "    \n",
    "    At the end of the computation, vector psi contains all the information about hydrophobic-hydrophobic contacts:\n",
    "    psi(k, j, p) = |1> iff amino acid k and amino acid j, where j is p-neighbour of k, are both hydrophobic and have \n",
    "    a contact in the lattice but not in the sequence.\n",
    "    '''\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    # star_matrix: 3D array to hold the |x_star> and y_star indices\n",
    "    # r_matrix: 3D array to hold the indices for r_reg\n",
    "    # delta_matrix: 3D array to hold pairs of delta indices\n",
    "    # psi_matrix: 3D array to hold psi indices. \n",
    "    # All indexed in the coding convention, i.e. from 0 to n-1 instead of from 1 to n. \n",
    "     \n",
    "    star_matrix = [[[0 for k in range(length)] for j in range(length)] for p in range(4)]\n",
    "    r_matrix = [[[0 for k in range(length)] for j in range(length)] for p in range(4)]\n",
    "    delta_matrix = [[[0 for k in range(length)] for j in range(length)] for p in range(4)]\n",
    "    psi_matrix = [[[0 for k in range(length)] for j in range(length)] for p in range(4)]\n",
    "    \n",
    "    star_index, r_index, delta_index, psi_index = 0, 0, 0, 0\n",
    "    \n",
    "    for k in range(length-3):\n",
    "        for j in range(k+3, length):\n",
    "            for p in range(4):\n",
    "                \n",
    "                star_matrix[k][j][p] = [m for m in range(star_index, star_index + register_length)]\n",
    "                star_index += register_length\n",
    "                \n",
    "                r_matrix[k][j][p] = [m for m in range(r_index, r_index + 2*register_length + 1)]\n",
    "                r_index += 2 * register_length + 1\n",
    "                \n",
    "                delta_matrix[k][j][p] = (delta_index, delta_index+1)\n",
    "                delta_index += 2\n",
    "                \n",
    "                psi_matrix[k][j][p] = psi_index\n",
    "                psi_index += 1\n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    #The code below has not been tested within the algorithm but the indices are assigned correctly and so it should work\n",
    "   \n",
    "    \n",
    "    for k in range(length-3):\n",
    "        for j in range(k+3, length):\n",
    "            for p in range(4):\n",
    "                \n",
    "                for r in range(register_length):\n",
    "                    qc.cx(x_coord[k+r], x_star[star_matrix[k][j][p][r]])\n",
    "                    qc.cx(y_coord[k+r], y_star[star_matrix[k][j][p][r]])\n",
    "                    \n",
    "                if p == 0:\n",
    "                    neighbour_north = update_coordinates(binary_one, [y_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_north)\n",
    "                \n",
    "                                                                      \n",
    "                if p == 1:\n",
    "                    neighbour_east = update_coordinates(binary_one, [x_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_east)\n",
    "                                                                      \n",
    "                if p == 2:\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                    neighbour_south = update_coordinates(binary_one, [y_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_south)\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                     \n",
    "                if p == 3:\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                    neighbour_west = update_coordinates(binary_one, [x_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_west)\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                \n",
    "                for r in range(register_length):\n",
    "                    qc.cx(x_coord[j+r], x_star[star_matrix[k][j][p][r]])\n",
    "                    qc.cx(y_coord[j+r], y_star[star_matrix[k][j][p][r]])\n",
    "                  \n",
    "                \n",
    "                # rows 13 to 19\n",
    "                for r in range(register_length):\n",
    "                    qc.x(y_star[star_matrix[k][j][p][r]])\n",
    "                    qc.x(x_star[star_matrix[k][j][p][r]])\n",
    "                \n",
    "                qc.ccx(y_star[star_matrix[k][j][p][0]], r_reg[r_matrix[k][j][p][0]], r_reg[r_matrix[k][j][p][1]])\n",
    "                    \n",
    "                for r in reversed(range(1, register_length)):\n",
    "                    qc.ccx(y_star[star_matrix[k][j][p][r]], r_reg[r_matrix[k][j][p][r]], r_reg[r_matrix[k][j][p][r+1]])\n",
    "                    \n",
    "                for r in reversed(range(register_length, 2*register_length)):\n",
    "                    qc.ccx(x_star[star_matrix[k][j][p][r-register_length]], r_reg[r_matrix[k][j][p][r]], \n",
    "                           r_reg[r_matrix[k][j][p][r+1]])\n",
    "                    \n",
    "                for r in range(register_length):\n",
    "                    qc.x(y_star[star_matrix[k][j][p][r]])\n",
    "                    qc.x(x_star[star_matrix[k][j][p][r]]) \n",
    "                     \n",
    "                \n",
    "                # row 20 onwards\n",
    "                qc.ccx(sequence[k], r_reg[r_matrix[k][j][p][2*register_length]], delta[delta_matrix[k][j][p][0]])\n",
    "                qc.ccx(sequence[j], delta[delta_matrix[k][j][p][0]], delta[delta_matrix[k][j][p][1]])\n",
    "                \n",
    "                qc.cx(delta[delta_matrix[k][j][p][1]], psi[psi_matrix[k][j][p]])\n",
    "                \n",
    "\n",
    "    \n",
    "\n",
    " \n",
    "\n",
    "\n",
    "            \n",
    "            \n",
    "    '''\n",
    "    Uncomputing calculation of H-H contacts (stored in vector psi) to free up qubits\n",
    "    '''\n",
    "    #This part has not been tested due to running out of memory. But should be correct. \n",
    "    \n",
    "    \n",
    "    for k in reversed(range(length-3)):\n",
    "        for j in reversed(range(k+3, length)):\n",
    "            for p in reversed(range(3)):\n",
    "                \n",
    "                qc.ccx(sequence[j], delta[delta_matrix[k][j][p][0]], delta[delta_matrix[k][j][p][1]])\n",
    "                qc.ccx(sequence[k], r_reg[r_matrix[k][j][p][2*register_length]], delta[delta_matrix[k][j][p][0]])\n",
    "                \n",
    "                for r in range(register_length):\n",
    "                    qc.x(y_star[star_matrix[k][j][p][r]])\n",
    "                    qc.x(x_star[star_matrix[k][j][p][r]])  \n",
    "                \n",
    "                for r in range(register_length, 2*register_length):\n",
    "                    qc.ccx(x_star[star_matrix[k][j][p][r-register_length]], r_reg[r_matrix[k][j][p][r]], \n",
    "                           r_reg[r_matrix[k][j][p][r+1]])\n",
    "                  \n",
    "                for r in range(1, register_length):\n",
    "                    qc.ccx(y_star[star_matrix[k][j][p][r]], r_reg[r_matrix[k][j][p][r]], r_reg[r_matrix[k][j][p][r+1]]) \n",
    "                \n",
    "                qc.ccx(y_star[star_matrix[k][j][p][0]], r_reg[r_matrix[k][j][p][0]], r_reg[r_matrix[k][j][p][1]]) \n",
    "                \n",
    "                for r in range(register_length):\n",
    "                    qc.x(y_star[star_matrix[k][j][p][r]])\n",
    "                    qc.x(x_star[star_matrix[k][j][p][r]])\n",
    "                    \n",
    "                for r in reversed(range(register_length)):\n",
    "                    qc.cx(x_coord[j+r], x_star[star_matrix[k][j][p][r]])\n",
    "                    qc.cx(y_coord[j+r], y_star[star_matrix[k][j][p][r]])    \n",
    "                \n",
    "                if p == 0:\n",
    "                    neighbour_north = update_coordinates(binary_one, [y_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_north)\n",
    "                \n",
    "                                                                      \n",
    "                if p == 1:\n",
    "                    neighbour_east = update_coordinates(binary_one, [x_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_east)\n",
    "                                                                      \n",
    "                if p == 2:\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                    neighbour_south = update_coordinates(binary_one, [y_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_south)\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                     \n",
    "                if p == 3:\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                    neighbour_west = update_coordinates(binary_one, [x_star[star_matrix[k][j][p][r]] \n",
    "                                                                      for r in range(register_length)])\n",
    "                    qc.append(subinst_rca, neighbour_west)\n",
    "                    for i in range(summand_register_length-1):\n",
    "                        qc.x(binary_one[i])\n",
    "                \n",
    "                for r in reversed(range(register_length)):\n",
    "                    qc.cx(x_coord[k+r], x_star[star_matrix[k][j][p][r]])\n",
    "                    qc.cx(y_coord[k+r], y_star[star_matrix[k][j][p][r]])\n",
    "\n",
    "                \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    # This part has been tested and performs correctly\n",
    "    '''\n",
    "    # Reversing coordinate calculation \n",
    "    '''\n",
    "    b = b-2\n",
    "\n",
    "    for d in range (length, 1, -1):\n",
    "\n",
    "        \n",
    "        # Uncomputing the western neighbour of site d-1 (|w> = |11>)\n",
    "        \n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "        neighbour_west = calculate_neighbour(x_coord, b, d)\n",
    "        qc.append(subinst.inverse(), neighbour_west)\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "              \n",
    "\n",
    "        \n",
    "        # Uncomputing the eastern neighbour of site d-1 (|w> = |01>)\n",
    "        \n",
    "        qc.x(w[b])\n",
    "        neighbour_east = calculate_neighbour(x_coord, b, d)\n",
    "        qc.append(subinst.inverse(), neighbour_east)\n",
    "        qc.x(w[b])\n",
    "\n",
    "        \n",
    "        # Uncomputing the northern neighbour of site d-1 (|w> = |00>)\n",
    "        \n",
    "        qc.x(w[b])\n",
    "        qc.x(w[b+1])\n",
    "        neighbour_north = calculate_neighbour(y_coord, b, d)\n",
    "        qc.append(subinst.inverse(), neighbour_north)\n",
    "        qc.x(w[b])\n",
    "        qc.x(w[b+1])\n",
    "        \n",
    "        \n",
    "        \n",
    "        # Uncomputing the southern neighbour of site d-1 (|w> = |10>)\n",
    "        \n",
    "        qc.x(w[b+1])\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "        neighbour_south = calculate_neighbour(x_coord, b, d)\n",
    "        qc.append(subinst.inverse(), neighbour_south)\n",
    "        for i in range(summand_register_length-1):\n",
    "            qc.x(binary_one[i])\n",
    "        qc.x(w[b+1])\n",
    "\n",
    "        b = b-2\n",
    "\n",
    "        for q in range(summand_register_length-1, -1, -1):\n",
    "            qc.cx(x_coord[(d-2)*summand_register_length+q],x_coord[(d-1)*summand_register_length+q])\n",
    "            qc.cx(y_coord[(d-2)*summand_register_length+q],y_coord[(d-1)*summand_register_length+q])\n",
    "\n",
    "     \n",
    "\n",
    "\n",
    "\n",
    "            \n",
    "            \n",
    "            \n",
    "            \n",
    "            \n",
    "    '''\n",
    "    SUMMING UP ENERGY VALUES IN PSI = THIS IS THE CRUCIAL ORACLE PART\n",
    "    STORING THE SUM OF ENERGIES IN VECTOR Z_{length-3, 3, j}\n",
    "    FLIPPING THE BIT ON GROVER QUBIT g FOR THE LARGEST ENERGY VALUE \n",
    "    \n",
    "    \n",
    "    p0_matrix = [ [ 0 for m in range(4 * (length-4) + 2) ] for n in range(length-2) ]\n",
    "\n",
    "    kj = 0\n",
    "    for k in range(1, length-2):\n",
    "        for j in range(4 * (k-1) + 2):\n",
    "            p0_matrix[k][j] = kj\n",
    "            kj += 1\n",
    "            \n",
    "    p1_matrix = [ [ 0 for m in range(4 * (length-4) + 3) ] for n in range(length-2) ]\n",
    "\n",
    "    for k in range(1, length-2):\n",
    "        for j in range(4 * (k-1) + 3):\n",
    "            p1_matrix[k][j] = kj\n",
    "            kj += 1\n",
    "    \n",
    "    p2_matrix = [ [ 0 for m in range(4 * (length-4) + 4) ] for n in range(length-2) ]\n",
    "\n",
    "    for k in range(1, length-2):\n",
    "        for j in range(4 * (k-1) + 4):\n",
    "            p2_matrix[k][j] = kj\n",
    "            kj += 1\n",
    "    \n",
    "    p3_matrix = [ [ 0 for m in range(4 * (length-4) + 5) ] for n in range(length-2) ]\n",
    "\n",
    "    for k in range(1, length-2):\n",
    "        for j in range(4 * (k-1) + 5):\n",
    "            p3_matrix[k][j] = kj\n",
    "            kj += 1\n",
    "    \n",
    "    ''' \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "    \n",
    "    # Adding up the energies. \n",
    "    # The energies for each conformation are stored in qubits z_{length-3, 3, j}.\n",
    "    # If j = 0, there are no HH contacts and the energy is 0.\n",
    "    # If j = 1, there is 1 HH contact and, hence, the energy is 1. And so on. \n",
    "    \n",
    "    \n",
    "    \n",
    "    for k in range(length-3):\n",
    "        for p in range(4):\n",
    "            \n",
    "            multiplexer = [psi[psi_matrix[k][i][p]] for i in range(k+3, length)]\n",
    "                \n",
    "            for j in reversed(range(4*k + p + 1)):\n",
    "                \n",
    "                \n",
    "                if p == 0:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.cx(z_and, z[p0_matrix[k+1][j+1]])\n",
    "                    # uncomputing to reuse the z_or and z_and ancillas\n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.cx(z_and, z[p0_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                elif p == 1:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p1_matrix[k+1][j+1]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p1_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "\n",
    "                elif p == 2:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p2_matrix[k+1][j+1]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p2_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                elif p == 3:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p3_matrix[k+1][j+1]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p3_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    # This is the crucial part of the oracle\n",
    "    # WE ASSUME j = 1 SINCE WE ASSUME THAT OUR SEQUENCE IS 01001, \n",
    "    # i.e. second and fifth amino acids are hydrophobic, first, third and fourth are hydrophilic.\n",
    "    # PARAMETER j IS ADJUSTABLE.\n",
    "    \n",
    "    j = 1\n",
    "    qc.cx(z[p3_matrix[length-3][j]], g)\n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    # Uncomputing the summation of energy\n",
    "    \n",
    "    for k in reversed(range(length-3)):\n",
    "        for p in reversed(range(4)):\n",
    "            \n",
    "            multiplexer = [psi[psi_matrix[k][i][p]] for i in range(k+3, length)]\n",
    "                \n",
    "            for j in range(4*k + p + 1):\n",
    "                \n",
    "                \n",
    "                if p == 0:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.cx(z_and, z[p0_matrix[k+1][j+1]])\n",
    "                    # uncomputing to reuse the z_or and z_and ancillas\n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.cx(z_and, z[p0_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p3_matrix[k][3]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                elif p == 1:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p1_matrix[k+1][j+1]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p1_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p0_matrix[k+1][j]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "\n",
    "                elif p == 2:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p2_matrix[k+1][j+1]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p2_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p1_matrix[k+1][j]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                elif p == 3:\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p3_matrix[k+1][j+1]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    \n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    qc.x(z_or)\n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.cx(z_and, z[p3_matrix[k+1][j]])\n",
    "                    \n",
    "                    qc.ccx(z_or, z[p2_matrix[k+1][j]], z_and)\n",
    "                    qc.x(z_or)\n",
    "                    qc.mcx(multiplexer, z_or)\n",
    "                    for index, qubit in enumerate(multiplexer):\n",
    "                        qc.x(multiplexer[index])    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "    '''\n",
    "    SECTION: IDENTIFYING THE CONFORMATION WITH THE MINIMAL ENERGY\n",
    "    \n",
    "    Using Grover's diffusion operation on register w\n",
    "    '''\n",
    "    \n",
    "    for i in range(len(w)):\n",
    "        qc.h(w[i])\n",
    "        qc.x(w[i])\n",
    "        \n",
    "    \n",
    "    multiplexer = [w[i] for i in range(len(w) - 1)]\n",
    "    \n",
    "    qc.h(w[len(w) - 1]) \n",
    "    qc.mcx(multiplexer, w[len(w) - 1])\n",
    "    qc.h(w[len(w) - 1])\n",
    "\n",
    "\n",
    "    for i in range(len(w)):\n",
    "        qc.x(w[i])\n",
    "        qc.h(w[i])  \n",
    "    \n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "'''\n",
    "MEASUREMENT.\n",
    "WE ASSUME THAT THE AMINO ACID SEQUENCE IS PHPPH = 01001. \n",
    "HENCE, WE ARE EXPECTING TO SEE THE FOLLOWING VALUES IN THE CHART:\n",
    "\n",
    "01011011\n",
    "01101111\n",
    "00000110\n",
    "00001110\n",
    "11111001\n",
    "11110001\n",
    "10100100\n",
    "10101100\n",
    "'''\n",
    "\n",
    "conformations = []\n",
    "\n",
    "for i in range(length):\n",
    "    conformations.append(sequence[i])\n",
    "for i in range(2*(length-1)):\n",
    "    conformations.append(w[i])\n",
    "\n",
    "# Reverse the order in which the output is shown so that it can be read from left to right.\n",
    "conformations.reverse()\n",
    "\n",
    "qc.measure(conformations, output)\n",
    "\n",
    "from qiskit import Aer, execute\n",
    "#from qiskit.providers.aer import QasmSimulator, StatevectorSimulator\n",
    "simulator = Aer.get_backend('qasm_simulator')\n",
    "result = execute(qc, backend = simulator, shots = 10).result()\n",
    "counts = result.get_counts()\n",
    "#statevector = result.get_statevector(qc)\n",
    "#print(statevector)\n",
    "#from qiskit.tools.visualization import plot_histogram\n",
    "#plot_histogram(counts)\n",
    "print(counts, sep='\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb38490a",
   "metadata": {},
   "outputs": [],
   "source": [
    "qc.draw('mpl')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}