Molecular dynamics simulation with Qamuy#
Install Qamuy Client SDK#
If you are running this notebook on Google Colaboratory, then install Qamuy Client SDK and login to Qamuy by running the following 2 cells. Otherwise, run commands in the following 2 cells on a terminal since they require input from standard input, which cannot be handled on Jupyter Notebook.
[ ]:
!python -m pip install qamuy-client --extra-index-url https://download.qamuy.qunasys.com/simple/
Define Methods for MD#
[1]:
import os
from functools import reduce
import numpy as np
import argparse
import yaml
import matplotlib.pyplot as plt
from qamuy.client import Client
from qamuy.utils.file_io import save, load_input
# You can fill in your e-mail address and password.
client = Client(email_address="YOUR_EMAIL_ADDRESS", password="YOUR_PASSWORD")
ao = 5.2917721067e-11 # bohr
NA = 6.022140857e23 # mol-1
me = 9.10938356e-31 # kg / a. u.
hbarOEh = 2.418884326509e-17 # s / a. u.
kt = hbarOEh * 1.0e15 # a.u. to fs
kr = ao * 1.0e10 # a.u. to A
mass_list = {
"H": 1.00782504,
"C": 12.00000000,
"N": 14.00307401,
"O": 15.99491464,
"F": 18.9984032,
"Cl": 35.453,
}
mass_list_au = {key: val * (1.0e-3 / NA / me) for key, val in mass_list.items()}
folder = "MD/"
temp_folder = folder + "TEMP/"
output_folder = folder + "OUTPUT/"
qamuy_setting_file = folder + "qamuy_setting.yaml"
qamuy_input_file = temp_folder + "qamuy_input.yaml"
qamuy_output_file = temp_folder + "qamuy_output.yaml"
output_file = output_folder + "output.log"
output_energy_file = output_folder + "energy.log"
output_xyz_file = output_folder + "geometry.xyz"
output_vel_file = output_folder + "velocity.log"
init_xyz_file = folder + "init.xyz"
init_vel_file = folder + "init.vel"
class Setting:
def __init__(self, time_step, time_limit, state, classical=False):
self.time_step = time_step / kt
self.time_limit = time_limit / kt
self.qamuy_setting = self.read_qamuy_setting()
self.classical = classical
self.client = client
self.state = state
def read_qamuy_setting(self) -> dict:
with open(qamuy_setting_file, "r") as f:
return yaml.safe_load(f)
def write_qamuy_input(self):
with open(qamuy_input_file, "w") as f:
yaml.dump(self.qamuy_setting, f, indent=2)
def read_qamuy_output(self):
with open(qamuy_output_file, "r") as f:
output = yaml.safe_load(f)
if self.classical:
results = output["molecule_result"]["post_hf_results"][0]
opt_params = []
else:
results = output["molecule_result"]["quantum_device_result"]
opt_params = results["vqe_log"]["opt_params"]
properies = results["evaluated_properties"]
ene = properies[0]["energy"]["values"][self.state]["value"]
grad = properies[4]["gradient"]["values"][0]["value"]
qamuy_log = Qamuy_log(
opt_params=opt_params, potential=ene, force=-np.array(grad)
)
return qamuy_log
def run_qamuy(self, step):
atoms = step.atoms
rr = step.position
geometry = self.qamuy_setting["target_molecule"]["geometry"]
geometry["atoms"] = atoms
geometry["coordinates"] = [(rr[3 * i : 3 * i + 3] * kr).tolist() for i in range(len(atoms))]
self.qamuy_setting["ansatz"]["initial_parameter"] = (
step.qamuy_log.opt_params if step.qamuy_log else []
)
self.write_qamuy_input()
input_pb = load_input(qamuy_input_file)
job = self.client.submit(input_pb)
output_pb = self.client.wait_and_get_job_results([job])[0]
save(output_pb.output, filename=qamuy_output_file)
return self.read_qamuy_output()
class Qamuy_log:
def __init__(self, opt_params, potential, force):
self.opt_params = opt_params
self.potential = potential
self.force = force
class Step:
def __init__(
self, atoms, velocity, position, time=0.0, filename="output", qamuy_log=None
):
self.filename = filename
self.atoms = atoms
self.time = time
self.velocity = velocity
self.position = position
self.qamuy_log = qamuy_log
self.mass = np.array(
reduce(lambda ss, x: ss + [mass_list_au[x]] * 3, atoms, [])
)
def cal_PE(self) -> float:
return self.qamuy_log.potential
def cal_KE(self) -> float:
v = self.velocity
return 0.5 * sum(step.mass * v * v)
def cal_totE(self) -> float:
return self.cal_PE() + self.cal_KE()
def velocity_verlet(setting: Setting, step: Step) -> Step:
# all terms in a.u.
dt = setting.time_step
r = step.position
v = step.velocity
F = step.qamuy_log.force
m = step.mass
step.time = step.time + dt
step.position = r + v * dt + 0.5 * dt * dt / m * F
q_log = setting.run_qamuy(step)
if q_log is not None:
new_F = q_log.force
step.velocity = v + 0.5 * dt / step.mass * (F + new_F)
step.qamuy_log = q_log
return step
def run_MD(setting: Setting, step: Step):
os.makedirs(temp_folder, exist_ok=True)
os.makedirs(output_folder, exist_ok=True)
step.qamuy_log = setting.run_qamuy(step)
dt = setting.time_step
N = int((setting.time_limit - step.time) / dt)
write_outputs(step, "w")
for i in range(N):
step = velocity_verlet(setting, step)
write_outputs(step, "a")
def write_outputs(step, mode):
atoms = step.atoms
t = step.time
length = len(atoms)
# write xyz
header = "%2d\n% 3.15e %5.6f\n" % (length, step.qamuy_log.potential, t * kt)
dd = reduce(
lambda ss, i: ss
+ "%2s % 3.15e % 3.15e % 3.15e\n"
% ((atoms[i],) + tuple(step.position[3 * i : 3 * i + 3] * kr)),
range(length),
header,
)
with open(output_xyz_file, mode) as f:
f.write(dd)
# write velocity
dd = reduce(
lambda ss, i: ss
+ "%2s % 3.15e % 3.15e % 3.15e\n"
% ((atoms[i],) + tuple(step.velocity[3 * i : 3 * i + 3])),
range(length),
header,
)
with open(output_vel_file, mode) as f:
f.write(dd)
# write energy
if mode == "w":
with open(output_energy_file, mode) as f:
f.write(
"# %7s %19s %19s %19s\n"
% ("Time", "Potential Eenergy", "Kinetic Eenergy", "Total Eenergy")
)
with open(output_energy_file, "a") as f:
f.write(
"% 5.6f % 3.15e % 3.15e % 3.15e\n"
% (t * kt, step.cal_PE(), step.cal_KE(), step.cal_totE())
)
def read_initial_file():
with open(init_xyz_file, "r") as f:
data = f.read()
lines = data.split("\n")
n = int(lines[0])
atoms, coords = zip(
*[(l.split()[0], [float(x) for x in l.split()[1:4]]) for l in lines[2 : n + 2]]
)
with open(init_vel_file, "r") as f:
data = f.read()
vels = np.array([float(x) for x in data.split()])
return list(atoms), (np.array(coords).ravel() / kr), vels
Setting for VQE#
[2]:
os.makedirs(folder, exist_ok=True)
setting = (
"""target_molecule:
geometry:
atoms: []
coordinates: []
basis: 6-31g
cas:
active_orb: 2
active_ele: 2
charge: 0
multiplicity: 1
num_excited_states: 0
# Set Chemical Properties to Be Calculated
output_chemical_properties:
- gradient:
type: ANALYTICAL_PYSCF
parameters:
- state: 0
# Set VQE Solver Type
solver:
type: VQE
# Set Qubit Mapping
mapping:
type: JORDAN_WIGNER
# Set Cost Function
cost_function:
type: SIMPLE
particle_number_weight: 10.0
s2_number_weight: 10.0
sz_number_weight: 10.0
# Set Ansatz
ansatz:
type: SYMMETRY_PRESERVING
initial_parameter: []
use_random_initial_guess: true
is_state_real: true
reference_state: RHF
depth: 5
# Set Optimizer
optimizer:
type: BFGS
ftol: 1.0e-06
gtol: 1.0e-06
max_fev: 100000
max_iter: 100000
max_run: 100000
"""
)
with open(qamuy_setting_file, "w") as f:
f.write(setting)
Generate initial nulcear coodinates and velocities#
[3]:
# coordinates in Å
xyz = (
"""3
H -1.5 -1.5 0.0
O 0.0 0.0 0.0
H 1.5 -1.5 0.0
"""
)
with open(init_xyz_file, "w") as f:
f.write(xyz)
# velocities in atomic unit
vel = (
"""0.005 0.0 0.0
0.0 0.0 0.0
-0.005 0.0 0.0
"""
)
with open(init_vel_file, "w") as f:
f.write(vel)
[4]:
import py3Dmol
view = py3Dmol.view(width=400, height=300)
view.addModel(xyz, "xyz")
view.setStyle({"sphere": {}})
view.show()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol
Setting for MD#
[5]:
time_step = 0.5 # fs
time_limit = 5.0 # fs
init_state = 0 # ground state
setting = Setting(
time_step=time_step,
time_limit=time_limit,
state=init_state,
)
Run MD#
[6]:
atoms, coords, vels = read_initial_file()
step = Step(atoms=atoms, position=coords, velocity=vels)
run_MD(setting, step)
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 27.6s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 27.6s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 21.7s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 21.7s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 16.3s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 16.3s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 21.3s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 21.3s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 21.3s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 21.3s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 10.9s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 11.0s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 11.0s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 11.0s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 11.1s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 11.1s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 16.3s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 16.3s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 21.3s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 21.3s finished
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done 1 tasks | elapsed: 26.8s
[Parallel(n_jobs=-1)]: Done 1 out of 1 | elapsed: 26.8s finished
Plot Results#
[7]:
def plot_energy(ff):
es = np.loadtxt(ff, comments="#").T
plt.plot(es[0], es[1], label="Potential energy")
plt.plot(es[0], es[3], label="Total energy")
plt.xlabel("Time / fs")
plt.ylabel("Energy / Hartree")
plt.legend()
plt.show()
plot_energy(output_energy_file)
[8]:
with open(output_xyz_file) as f:
output_xyz = f.read()
view = py3Dmol.view(width=400, height=300)
view.addModelsAsFrames(output_xyz, "xyz")
view.setStyle({"sphere": {}})
view.animate({'loop': "forward"})
view.show()
You appear to be running in JupyterLab (or JavaScript failed to load for some other reason). You need to install the 3dmol extension:
jupyter labextension install jupyterlab_3dmol