noisy_simulation.py

import numpy as np
import pandas as pd
import csv
import matplotlib.pyplot as plt
from copy import deepcopy
import cmath

num_rot = 2

df1 = pd.read_csv("energy_files/one_body_terms.csv")
q = df1['E_ii'].values
num = len(q)
N_res = int(num/num_rot)

df = pd.read_csv("energy_files/two_body_terms.csv")
v = df['E_ij'].values
numm = len(v)
Q = np.zeros((num,num))
n = 0

for i in range(num-2):
    if i%2 == 0:
        Q[i][i+2] = deepcopy(v[n])
        Q[i+2][i] = deepcopy(v[n])
        Q[i][i+3] = deepcopy(v[n+1])
        Q[i+3][i] = deepcopy(v[n+1])
        n += 2
    elif i%2 != 0:
        Q[i][i+1] = deepcopy(v[n])
        Q[i+1][i] = deepcopy(v[n])
        Q[i][i+2] = deepcopy(v[n+1])
        Q[i+2][i] = deepcopy(v[n+1])
        n += 2

print("q: \n", q)
print("Q: \n", Q)
H = np.zeros((num,num))

for i in range(num):
    for j in range(num):
        if i != j:
            H[i][j] = np.multiply(0.25, Q[i][j])

for i in range(num):
    H[i][i] = -(0.5 * q[i] + sum(0.25 * Q[i][j] for j in range(num) if j != i))

# add penalty terms to the matrix so as to discourage the selection of two rotamers on the same residue - implementation of the Hammings constraint
def add_penalty_term(M, penalty_constant, residue_pairs):
    for i, j in residue_pairs:
        M[i][j] += penalty_constant
    return M

def generate_pairs(N):
    pairs = [(i, i+1) for i in range(0, 2*N, 2)]
    return pairs

P = 6
pairs = generate_pairs(N_res)

M = deepcopy(H)
M = add_penalty_term(M, P, pairs)


## Classical optimisation:
from scipy.sparse.linalg import eigsh
from scipy.linalg import eig

Z_matrix = np.array([[1, 0], [0, -1]])
identity = np.eye(2)

def construct_operator(qubit_indices, num):
    operator = np.eye(1)
    for qubit in range(num):
        if qubit in qubit_indices:
            operator = np.kron(operator, Z_matrix)
        else:
            operator = np.kron(operator, identity)
    return operator

C = np.zeros((2**num, 2**num))

for i in range(num):
    operator = construct_operator([i], num)
    C += H[i][i] * operator

for i in range(num):
    for j in range(i+1, num):
        operator = construct_operator([i, j], num)
        C += H[i][j] * operator

def create_hamiltonian(pairs, P, num):
    H_pen = np.zeros((2**num, 2**num))
    def tensor_term(term_indices):
        term = [Z_matrix if i in term_indices else identity for i in range(num)]
        result = term[0]
        for t in term[1:]:
            result = np.kron(result, t)
        return result
    
    for pair in pairs:
        term = tensor_term(pair)
        H_pen += P * term

    return H_pen

H_penalty = create_hamiltonian(pairs, P, num)
H_tot = C + H_penalty

eigenvalues, eigenvectors = eigsh(H_tot, k=num, which='SA')

## Quantum optimisation
#  Find minimum value using optimisation technique of QAOA
from qiskit import Aer, QuantumCircuit, execute
from qiskit_algorithms.minimum_eigensolvers import QAOA
from qiskit.quantum_info.operators import Operator, Pauli, SparsePauliOp
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler


def X_op(i, num):
    op_list = ['I'] * num
    op_list[i] = 'X'
    return SparsePauliOp(Pauli(''.join(op_list)))

def generate_pauli_zij(n, i, j):
    if i<0 or i >= n or j<0 or j>=n:
        raise ValueError(f"Indices out of bounds for n={n} qubits. ")
   
    pauli_str = ['I']*n

    if i == j:
        pauli_str[i] = 'Z'
    else:
        pauli_str[i] = 'Z'
        pauli_str[j] = 'Z'

    return Pauli(''.join(pauli_str))

q_hamiltonian = SparsePauliOp(Pauli('I'*num), coeffs=[0])

for i in range(num):
    for j in range(i+1, num):
        if M[i][j] != 0:
            pauli = generate_pauli_zij(num, i, j)
            op = SparsePauliOp(pauli, coeffs=[M[i][j]])
            q_hamiltonian += op

for i in range(num):
    pauli = generate_pauli_zij(num, i, i)
    Z_i = SparsePauliOp(pauli, coeffs=[M[i][i]])
    q_hamiltonian += Z_i

def format_sparsepauliop(op):
    terms = []
    labels = [pauli.to_label() for pauli in op.paulis]
    coeffs = op.coeffs
    for label, coeff in zip(labels, coeffs):
        terms.append(f"{coeff:.10f} * {label}")
    return '\n'.join(terms)

print(f"\nThe hamiltonian constructed using Pauli operators is: \n", format_sparsepauliop(q_hamiltonian))

mixer_op = sum(X_op(i,num) for i in range(num))

p = 1  # Number of QAOA layers
initial_point = np.ones(2 * p)
qaoa = QAOA(sampler=Sampler(), optimizer=COBYLA(), reps=p, mixer=mixer_op, initial_point=initial_point)
result = qaoa.compute_minimum_eigenvalue(q_hamiltonian)
print("\n\nThe result of the quantum optimisation using QAOA is: \n")
print('best measurement', result.best_measurement)
print('Optimal parameters: ', result.optimal_parameters)

# from qiskit_optimization.applications import max_cut
# optimal_params = result.optimal_parameters
# graph = [(0, 1, 1.0), (1, 2, 1.0), (2, 3, 1.0), (3, 0, 1.0)]
# qubit_op, offset = max_cut.get_operator(graph)
# circuit = qaoa.construct_circuit(optimal_params, qubit_op).decompose()
# depth = circuit.depth()
# print("Depth of the QAOA circuit:", depth)

k = 0
for i in range(num):
    k += 0.5 * q[i]

for i in range(num):
    for j in range(num):
        if i != j:
            k += 0.5 * 0.25 * Q[i][j]

print('\nThe ground state energy classically is: ', eigenvalues[0] + N_res*P + k)
print('The ground state energy with QAOA is: ', np.real(result.best_measurement['value']) + N_res*P + k)

from qiskit_aer.noise import NoiseModel
from qiskit_ibm_provider import IBMProvider
from qiskit_aer import AerSimulator
from qiskit.providers.fake_provider import FakeKolkata, FakeVigo
from qiskit_ibm_runtime import QiskitRuntimeService, Options, Session, Sampler

IBMProvider.save_account('25a4f69c2395dfbc9990a6261b523fe99e820aa498647f92552992afb1bd6b0bbfcada97ec31a81a221c16be85104beb653845e23eeac2fe4c0cb435ec7fc6b4', overwrite=True)
provider = IBMProvider()
available_backends = provider.backends()
print([backend.name for backend in available_backends])
service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.backend("ibmq_qasm_simulator")
noise_model = NoiseModel.from_backend(backend)
simulator = AerSimulator(noise_model = noise_model)
fake_backend = FakeKolkata()
noise_model = NoiseModel.from_backend(fake_backend)
options = Options()
options.simulator = {
    "noise_model": noise_model,
    "basis_gates": fake_backend.configuration().basis_gates,
    "coupling_map": fake_backend.configuration().coupling_map,
    "seed_simulator": 42
}
options.execution.shots = 1000
options.optimization_level = 0
options.resilience_level = 0

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    qaoa = QAOA(sampler=sampler, optimizer=COBYLA(), reps=p, mixer=mixer_op, initial_point=initial_point)
    result1 = qaoa.compute_minimum_eigenvalue(q_hamiltonian)

print("\n\nThe result of the noisy quantum optimisation using QAOA is: \n")
print('best measurement', result1.best_measurement)
print('Optimal parameters: ', result1.optimal_parameters)
print('The ground state energy with noisy QAOA is: ', np.real(result1.best_measurement['value']) + N_res*P + k)



## Creation and Annihilation operators Hamiltonian
num_qubits = N_res

## First Classically

H_s = np.zeros((2**num_qubits, 2**num_qubits), dtype=complex)
H_i = np.zeros((2**num_qubits, 2**num_qubits), dtype=complex)

X = np.array([[0, 1], [1, 0]])
Y = np.array([[0, -1j], [1j, 0]])

a = 0.5*(X + 1j*Y)
a_dagger = 0.5 *(X - 1j*Y) 

aad = a@a_dagger
ada = a_dagger@a

def extended_operator(n, qubit, op):
    ops = [identity if i != qubit else op for i in range(n)]
    extended_op = ops[0]
    for op in ops[1:]:
        extended_op = np.kron(extended_op, op)
    return extended_op

s = 0
for i in range(N_res):
    aad_extended = extended_operator(num_qubits, i, aad)
    ada_extended = extended_operator(num_qubits, i, ada)
    H_s += q[s] * aad_extended + q[s+1] * ada_extended 
    s += num_rot
    if s >= num:
        break

k = 0
for i in range(N_res-1):
    aad_extended = extended_operator(num_qubits, i, aad)
    ada_extended = extended_operator(num_qubits, i, ada)
    aad_extended1 = extended_operator(num_qubits, i+1, aad)
    ada_extended1 = extended_operator(num_qubits, i+1, ada)
    H_i += v[k] * aad_extended @ aad_extended1 + \
                v[k+1] * aad_extended @ ada_extended1 + \
                v[k+2] * ada_extended @ aad_extended1 + \
                v[k+3] * ada_extended @ ada_extended1
    k += num_rot**2
    if k >= numm:
        break

H_tt = H_i + H_s 
eigenvalue, eigenvector = eigsh(H_tt, k=num_qubits, which='SA')
print('\n\nThe ground state with the number operator classically is: ', eigenvalue[0])
# print('The classical eigenstate is: ', eigenvalue)

# ground_state = eig(H_tt)
# print('eig result:', ground_state)


## Mapping to qubits
H_self = SparsePauliOp(Pauli('I'* num_qubits), coeffs=[0])
H_int = SparsePauliOp(Pauli('I'* num_qubits), coeffs=[0]) 

def N_0(i, n):
    pauli_str = ['I'] * n
    pauli_str[i] = 'Z'
    z_op = SparsePauliOp(Pauli(''.join(pauli_str)), coeffs=[0.5])
    i_op = SparsePauliOp(Pauli('I'*n), coeffs=[0.5])
    return z_op + i_op

def N_1(i, n):
    pauli_str = ['I'] * n
    pauli_str[i] = 'Z'
    z_op = SparsePauliOp(Pauli(''.join(pauli_str)), coeffs=[-0.5])
    i_op = SparsePauliOp(Pauli('I'*n), coeffs=[0.5])
    return z_op + i_op

l = 0
for i in range(N_res):
    N_0i = N_0(i, num_qubits)
    N_1i = N_1(i, num_qubits)
    H_self += q[l] * N_0i + q[l+1] * N_1i 
    l += num_rot
    if l >= num:
        break

j = 0
for i in range(N_res-1):
    N_0i = N_0(i, num_qubits)
    N_1i = N_1(i, num_qubits)
    N_0j = N_0(i+1, num_qubits)
    N_1j = N_1(i+1, num_qubits)
    H_int += v[j] * N_0i @ N_0j + v[j+1] * N_0i @ N_1j + v[j+2] * N_1i @ N_0j + v[j+3] * N_1i @ N_1j
    j += num_rot**2
    if j >= numm:
        break

H_gen = H_int + H_self

mixer_op = sum(X_op(i,num_qubits) for i in range(num_qubits))
p = 1
initial_point = np.ones(2*p)
qaoa1 = QAOA(sampler=Sampler(), optimizer=COBYLA(), reps=p, mixer=mixer_op, initial_point=initial_point)
result_gen = qaoa1.compute_minimum_eigenvalue(H_gen)
print("\nThe result of the quantum optimisation using QAOA with the number operators is: \n")
print('best measurement', result_gen.best_measurement)
print("Optimal parameters: ", qaoa1.optimal_params)
print('\nThe ground state energy with QAOA is: ', np.real(result_gen.best_measurement['value']), '\n')
print(result_gen)


IBMProvider.save_account('25a4f69c2395dfbc9990a6261b523fe99e820aa498647f92552992afb1bd6b0bbfcada97ec31a81a221c16be85104beb653845e23eeac2fe4c0cb435ec7fc6b4', overwrite=True)
provider = IBMProvider()
available_backends = provider.backends()
print([backend.name for backend in available_backends])
service = QiskitRuntimeService(channel="ibm_quantum")
backend = service.backend("ibmq_qasm_simulator")
noise_model = NoiseModel.from_backend(backend)
simulator = AerSimulator(noise_model = noise_model)
fake_backend = FakeKolkata()
noise_model = NoiseModel.from_backend(fake_backend)
options = Options()
options.simulator = {
    "noise_model": noise_model,
    "basis_gates": fake_backend.configuration().basis_gates,
    "coupling_map": fake_backend.configuration().coupling_map,
    "seed_simulator": 42
}
options.execution.shots = 1000
options.optimization_level = 0
options.resilience_level = 0

with Session(service=service, backend=backend):
    sampler = Sampler(options=options)
    qaoa2 = QAOA(sampler=sampler, optimizer=COBYLA(), reps=p, mixer=mixer_op, initial_point=initial_point)
    result2 = qaoa2.compute_minimum_eigenvalue(H_gen)

print("\n\nThe result of the noisy quantum optimisation using QAOA with number operators is: \n")
print('best measurement', result2.best_measurement)
print('Optimal parameters: ', result2.optimal_parameters)
print('The ground state energy with noisy QAOA number operators is: ', np.real(result2.best_measurement['value']) + N_res*P + k)