# !pip install qoai_qiskit
# from qoai_qiskit import load_account
# provider = load_account()
# from qiskit import
# provider = IBMQ.get_provider()
!pip install qiskit_research
from qiskit_research.protein_folding.interactions.random_interaction import (
RandomInteraction,
)
from qiskit_research.protein_folding.interactions.miyazawa_jernigan_interaction import (
MiyazawaJerniganInteraction,
)
from qiskit_research.protein_folding.peptide.peptide import Peptide
from qiskit_research.protein_folding.protein_folding_problem import (
ProteinFoldingProblem,
)
from qiskit_research.protein_folding.penalty_parameters import PenaltyParameters
from qiskit.utils import algorithm_globals, QuantumInstance
algorithm_globals.random_seed = 23
Qiskit Nature
# Before, using qiskit-nature:
from qiskit_nature.problems.sampling.protein_folding.protein_folding_problem import (
ProteinFoldingProblem,
)
Now using qiskit research
# Now, using qiskit-research:
from qiskit_research.protein_folding.protein_folding_problem import (
ProteinFoldingProblem,
)
Protein Main Chain
main_chain = "APRLRFY"
Side chain
side_chains = [""] * 7
Interaction between Aminoacids
random_interaction = RandomInteraction()
mj_interaction = MiyazawaJerniganInteraction()
Physical Constraints
penalty_back = 10
penalty_chiral = 10
penalty_1 = 10
penalty_terms = PenaltyParameters(penalty_chiral, penalty_back, penalty_1)
Peptide Definition
peptide = Peptide(main_chain, side_chains)
Protein Folding Problem
protein_folding_problem = ProteinFoldingProblem(peptide, mj_interaction, penalty_terms)
qubit_op = protein_folding_problem.qubit_op()
print(qubit_op)
Using VQE with CVaR expectation value for the solution of the problem
from qiskit.circuit.library import RealAmplitudes
from qiskit.algorithms.optimizers import COBYLA
from qiskit.algorithms import NumPyMinimumEigensolver
from qiskit.algorithms.minimum_eigensolvers import SamplingVQE
from qiskit import execute, Aer
from qiskit.primitives import Sampler
# set classical optimizer
optimizer = COBYLA(maxiter=50)
# set variational ansatz
ansatz = RealAmplitudes(reps=1)
counts = []
values = []
def store_intermediate_result(eval_count, parameters, mean, std):
counts.append(eval_count)
values.append(mean)
# initialize VQE using CVaR with alpha = 0.1
vqe = SamplingVQE(
Sampler(),
ansatz=ansatz,
optimizer=optimizer,
aggregation=0.1,
callback=store_intermediate_result,
)
raw_result = vqe.compute_minimum_eigenvalue(qubit_op)
print(raw_result)
import matplotlib.pyplot as plt
fig = plt.figure()
plt.plot(counts, values)
plt.ylabel("Conformation Energy")
plt.xlabel("VQE Iterations")
fig.add_axes([0.44, 0.51, 0.44, 0.32])
plt.plot(counts[40:], values[40:])
plt.ylabel("Conformation Energy")
plt.xlabel("VQE Iterations")
plt.show()
Visualizing the answer
result = protein_folding_problem.interpret(raw_result=raw_result)
print(
"The bitstring representing the shape of the protein during optimization is: ",
result.turn_sequence,
)
print("The expanded expression is:", result.get_result_binary_vector())
From this sequence of turns we can get the cartesian coordinates of each of the aminoacids of the protein.
print(result.protein_shape_file_gen.get_xyz_data())
And finally, we can also plot the structure of the protein in 3D. Note that when rendered with the proper backend this plot can be interactively rotated.
fig = result.get_figure(title="Protein Structure", ticks=False, grid=True)
fig.get_axes()[0].view_init(10, 70)
And here is an example with side chains.
peptide = Peptide("APRLR", ["", "", "F", "Y", ""])
protein_folding_problem = ProteinFoldingProblem(peptide, mj_interaction, penalty_terms)
qubit_op = protein_folding_problem.qubit_op()
# set classical optimizer
optimizer = COBYLA(maxiter=50)
# set variational ansatz
ansatz = RealAmplitudes(reps=1)
counts = []
values = []
def store_intermediate_result(eval_count, parameters, mean, std):
counts.append(eval_count)
values.append(mean)
# initialize VQE using CVaR with alpha = 0.1
vqe = SamplingVQE(
Sampler(),
ansatz=ansatz,
optimizer=optimizer,
aggregation=0.1,
callback=store_intermediate_result,
)
raw_result = vqe.compute_minimum_eigenvalue(qubit_op)
result_2 = protein_folding_problem.interpret(raw_result=raw_result)
fig = result_2.get_figure(title="Protein Structure", ticks=False, grid=True)
fig.get_axes()[0].view_init(10, 60)