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)
../_images/_usecase_03_molecular_dynamics_16_0.png
[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