diff --git a/local_hw_copy.py b/local_hw_copy.py
new file mode 100644
index 00000000..3d40f371
--- /dev/null
+++ b/local_hw_copy.py
@@ -0,0 +1,384 @@
+# Point 2 of constraint studies for paper, Ising model with local penalties
+
+# Script to optimise the Hamiltonian, starting directly from the Ising Hamiltonian
+# %%
+import numpy as np
+import pandas as pd
+import time
+from copy import deepcopy
+
+num_rot = 2
+file_path = "RESULTS/localpenalty-QAOA/15res-2rot"
+
+########################### Configure the hamiltonian from the values calculated classically with pyrosetta ############################
+df1 = pd.read_csv("energy_files/one_body_terms.csv")
+q = df1['E_ii'].values
+num = len(q)
+N = int(num/num_rot)
+num_qubits = num
+
+print('Qii values: \n', q)
+
+df = pd.read_csv("energy_files/two_body_terms.csv")
+value = df['E_ij'].values
+Q = np.zeros((num,num))
+n = 0
+
+for i in range(0, num-2):
+    if i%2 == 0:
+        Q[i][i+2] = deepcopy(value[n])
+        Q[i+2][i] = deepcopy(value[n])
+        Q[i][i+3] = deepcopy(value[n+1])
+        Q[i+3][i] = deepcopy(value[n+1])
+        n += 2
+    elif i%2 != 0:
+        Q[i][i+1] = deepcopy(value[n])
+        Q[i+1][i] = deepcopy(value[n])
+        Q[i][i+2] = deepcopy(value[n+1])
+        Q[i+2][i] = deepcopy(value[n+1])
+        n += 2
+
+print('\nQij values: \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))
+
+print('\nH: \n', H)
+
+with open(file_path, "w") as file:
+    file.write(f"H : {H} \n")
+
+# 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
+
+P = 6
+
+def generate_pairs(N):
+    pairs = [(i, i+1) for i in range(0, 2*N, 2)]
+    return pairs
+
+pairs = generate_pairs(N)
+
+M = deepcopy(H)
+M = add_penalty_term(M, P, pairs)
+
+# %% ################################################ Classical optimisation ###########################################################
+from scipy.sparse.linalg import eigsh
+
+Z_matrix = np.array([[1, 0], [0, -1]])
+identity = np.eye(2)
+
+def construct_operator(qubit_indices, num_qubits):
+    operator = np.eye(1)
+    for qubit in range(num_qubits):
+        if qubit in qubit_indices:
+            operator = np.kron(operator, Z_matrix)
+        else:
+            operator = np.kron(operator, identity)
+    return operator
+
+C = np.zeros((2**num_qubits, 2**num_qubits))
+
+for i in range(num_qubits):
+    operator = construct_operator([i], num_qubits)
+    C += H[i][i] * operator
+
+for i in range(num_qubits):
+    for j in range(i+1, num_qubits):
+        operator = construct_operator([i, j], num_qubits)
+        C += H[i][j] * operator
+
+print('C :\n', C)
+
+def create_hamiltonian(pairs, P, num_qubits):
+    H_pen = np.zeros((2**num_qubits, 2**num_qubits))
+    def tensor_term(term_indices):
+        term = [Z_matrix if i in term_indices else identity for i in range(num_qubits)]
+        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_qubits)
+H_tot = C + H_penalty
+
+# Extract the ground state energy and wavefunction
+# using sparse representation so as to be able to generalise to larger systems
+eigenvalues, eigenvectors = eigsh(H_tot, k=num, which='SA')
+print("\n\nClassical optimisation results. \n")
+print("Ground energy eigsh: ", eigenvalues[0])
+print("ground state wavefuncion eigsh: ", eigenvectors[:,0])
+print('\n\n')
+
+with open(file_path, "a") as file:
+    file.write("\n\nClassical optimisation results.\n")
+    file.write(f"Ground energy eigsh: {eigenvalues[0]}\n")
+    file.write(f"Ground state wavefunction eigsh: {eigenvectors[:,0]}\n")
+
+
+# %% ############################################ Quantum optimisation ########################################################################
+from qiskit_algorithms.minimum_eigensolvers import QAOA
+from qiskit.quantum_info.operators import Pauli, SparsePauliOp
+from qiskit_algorithms.optimizers import COBYLA
+from qiskit.primitives import Sampler
+
+def X_op(i, num_qubits):
+    """Return an X Pauli operator on the specified qubit in a num-qubit system."""
+    op_list = ['I'] * num_qubits
+    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_qubits), coeffs=[0])
+
+for i in range(num_qubits):
+    for j in range(i+1, num_qubits):
+        if M[i][j] != 0:
+            pauli = generate_pauli_zij(num_qubits, i, j)
+            op = SparsePauliOp(pauli, coeffs=[M[i][j]])
+            q_hamiltonian += op
+
+for i in range(num_qubits):
+    pauli = generate_pauli_zij(num_qubits, 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))
+
+#the mixer in QAOA should be a quantum operator representing transitions between configurations
+mixer_op = sum(X_op(i,num_qubits) for i in range(num_qubits))
+p = 1  # Number of QAOA layers
+initial_point = np.ones(2 * p)
+
+# %%
+start_time = time.time()
+qaoa = QAOA(sampler=Sampler(), optimizer=COBYLA(), reps=p, mixer=mixer_op, initial_point=initial_point)
+result = qaoa.compute_minimum_eigenvalue(q_hamiltonian)
+end_time = time.time()
+
+print("\n\nThe result of the quantum optimisation using QAOA is: \n")
+print('best measurement', result.best_measurement)
+elapsed_time = end_time - start_time
+print(f"Local Simulation run time: {elapsed_time} seconds")
+print('\n\n')
+
+with open(file_path, "a") as file:
+    file.write("\n\nThe result of the quantum optimisation using QAOA is: \n")
+    file.write(f"'best measurement' {result.best_measurement}\n")
+    file.write(f"Local Simulation run time: {elapsed_time} seconds\n")
+
+
+# %% ############################################ Simulators ##########################################################################
+from qiskit_aer import Aer
+from qiskit_ibm_provider import IBMProvider
+from qiskit_aer.noise import NoiseModel
+from qiskit_aer.primitives import Sampler
+from qiskit.primitives import Sampler, BackendSampler
+from qiskit.transpiler import PassManager
+
+simulator = Aer.get_backend('qasm_simulator')
+provider = IBMProvider()
+available_backends = provider.backends()
+print("Available Backends:", available_backends)
+device_backend = provider.get_backend('ibm_torino')
+noise_model = NoiseModel.from_backend(device_backend)
+
+options= {
+    "noise_model": noise_model,
+    "basis_gates": simulator.configuration().basis_gates,
+    "coupling_map": simulator.configuration().coupling_map,
+    "seed_simulator": 42,
+    "shots": 1000,
+    "optimization_level": 3,
+    "resilience_level": 0
+}
+
+noisy_sampler = BackendSampler(backend=simulator, options=options, bound_pass_manager=PassManager())
+
+start_time1 = time.time()
+qaoa1 = QAOA(sampler=noisy_sampler, optimizer=COBYLA(), reps=p, mixer=mixer_op, initial_point=initial_point)
+result1 = qaoa1.compute_minimum_eigenvalue(q_hamiltonian)
+end_time1 = time.time()
+
+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']))
+elapsed_time1 = end_time1 - start_time1
+print(f"Aer Simulator run time: {elapsed_time1} seconds")
+print('\n\n')
+
+with open(file_path, "a") as file:
+    file.write("\n\nThe result of the noisy quantum optimisation using QAOA is: \n")
+    file.write(f"'best measurement' {result1.best_measurement}")
+    file.write(f"Optimal parameters: {result1.optimal_parameters}")
+    file.write(f"'The ground state energy with noisy QAOA is: ' {np.real(result1.best_measurement['value'])}")
+    file.write(f"Aer Simulator run time: {elapsed_time1} seconds")
+
+# %% ############################################# Hardware with QAOAAnastz ##################################################################
+from qiskit.circuit.library import QAOAAnsatz
+from qiskit_algorithms import SamplingVQE
+from qiskit_ibm_runtime import QiskitRuntimeService, Session, Sampler
+from qiskit import transpile, QuantumCircuit, QuantumRegister
+from qiskit.transpiler import CouplingMap, Layout
+
+service = QiskitRuntimeService()
+backend = service.backend("ibm_nazca")
+print('Coupling Map of hardware: ', backend.configuration().coupling_map)
+
+ansatz = QAOAAnsatz(q_hamiltonian, mixer_operator=mixer_op, reps=p)
+print('\n\nQAOAAnsatz: ', ansatz)
+
+target = backend.target
+# %%
+# real_coupling_map = backend.configuration().coupling_map
+# coupling_map = CouplingMap(couplinglist=real_coupling_map)
+
+def generate_linear_coupling_map(num_qubits):
+
+    coupling_list = [[i, i + 1] for i in range(num_qubits - 1)]
+    
+    return CouplingMap(couplinglist=coupling_list)
+
+# linear_coupling_map = generate_linear_coupling_map(num_qubits)
+# coupling_map = CouplingMap(couplinglist=[[0, 1],[0, 15], [1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7], [7, 8], [8, 9], [9, 10], [10, 11], [11, 12], [12, 13], [13, 14]])
+coupling_map = CouplingMap(couplinglist=[[0, 1], [0, 14], [1, 2], [3, 2], [3, 4], [4, 15], [5, 4], [6, 5], [6, 7], [7, 8], [8, 16], [9, 8], [9, 10], [10, 11], [12, 11], [13, 12], [15, 22], [17, 12], [18, 14], [18, 19], [20, 19], [20, 21], [20, 33], [21, 22], [23, 22], [24, 23], [25, 24], [25, 26], [26, 16], [27, 26], [27, 28], [28, 29]])
+qr = QuantumRegister(num_qubits, 'q')
+circuit = QuantumCircuit(qr)
+trivial_layout = Layout({qr[i]: i for i in range(num_qubits)})
+ansatz_isa = transpile(ansatz, backend=backend, initial_layout=trivial_layout, coupling_map=coupling_map,
+                       optimization_level= 3, layout_method='dense', routing_method='stochastic')
+print("\n\nAnsatz layout after explicit transpilation:", ansatz_isa._layout)
+
+hamiltonian_isa = q_hamiltonian.apply_layout(ansatz_isa.layout)
+print("\n\nAnsatz layout after transpilation:", hamiltonian_isa)
+
+# %%
+ansatz_isa.decompose().draw('mpl')
+
+op_counts = ansatz_isa.count_ops()
+total_gates = sum(op_counts.values())
+depth = ansatz_isa.depth()
+print("Operation counts:", op_counts)
+print("Total number of gates:", total_gates)
+print("Depth of the circuit: ", depth)
+
+# %%
+session = Session(backend=backend)
+print('\nhere 1')
+sampler = Sampler(backend=backend, session=session)
+print('here 2')
+qaoa2 = SamplingVQE(sampler=sampler, ansatz=ansatz_isa, optimizer=COBYLA(), initial_point=initial_point)
+print('here 3')
+result2 = qaoa2.compute_minimum_eigenvalue(hamiltonian_isa)
+
+print("\n\nThe result of the noisy quantum optimisation using QAOAAnsatz is: \n")
+print('best measurement', result2.best_measurement)
+print('Optimal parameters: ', result2.optimal_parameters)
+print('The ground state energy with noisy QAOA is: ', np.real(result2.best_measurement['value']))
+
+# %%
+jobs = service.jobs(session_id='crsnrga7jqmg008ze0eg')
+
+for job in jobs:
+    if job.status().name == 'DONE':
+        results = job.result()
+        print("Job completed successfully")
+else:
+    print("Job status:", job.status())
+
+# %%
+def get_best_measurement_from_sampler_result(sampler_result):
+    if not hasattr(sampler_result, 'quasi_dists') or not isinstance(sampler_result.quasi_dists, list):
+        raise ValueError("SamplerResult does not contain 'quasi_dists' as a list")
+
+    best_bitstring = None
+    highest_probability = -1
+
+    for quasi_distribution in sampler_result.quasi_dists:
+        for bitstring, probability in quasi_distribution.items():
+            if probability > highest_probability:
+                highest_probability = probability
+                best_bitstring = bitstring
+
+    return best_bitstring, highest_probability
+
+best_measurement, probability = get_best_measurement_from_sampler_result(results)
+print(f"Best measurement: {best_measurement} with probability {probability}")
+
+def index_to_bitstring(index, num_qubits):
+    bitstring = format(index, '0{}b'.format(num_qubits))
+    return bitstring
+
+bitstring = index_to_bitstring(best_measurement, num_qubits)
+print("The bitstring representation of the highest probability state is:", bitstring)
+
+total_usage_time = 0
+for job in jobs:
+    job_result = job.usage_estimation['quantum_seconds']
+    total_usage_time += job_result
+
+print(f"Total Usage Time Hardware: {total_usage_time} seconds")
+print('\n\n')
+
+with open(file_path, "a") as file:
+    file.write("\n\nThe result of the noisy quantum optimisation using QAOAAnsatz is: \n")
+    file.write(f"'best measurement' {result2.best_measurement}")
+    file.write(f"Optimal parameters: {result2.optimal_parameters}")
+    file.write(f"'The ground state energy with noisy QAOA is: ' {np.real(result2.best_measurement['value'])}")
+    file.write(f"Total Usage Time Hardware: {total_usage_time} seconds")
+    file.write(f"Total number of gates: {total_gates}\n")   
+    file.write(f"Depth of circuit: {depth}\n")
+
+# %%
+index = ansatz_isa.layout.final_index_layout() # Maps logical qubit index to its position in bitstring
+
+# measured_bitstring = result2.best_measurement['bitstring']
+measured_bitstring = bitstring
+original_bitstring = ['']*num_qubits
+
+for i, logical in enumerate(index):
+        original_bitstring[i] = measured_bitstring[logical]
+
+original_bitstring = ''.join(original_bitstring)
+print("Original bitstring:", original_bitstring)
+# %%