Geometry optimization and vibration analysis with Qamuy

At the end of this notebook, there is a part for visualization of molecules and their vibration modes.

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/

Import Qamuy

import qamuy.chemistry as qy
from qamuy.client import Client

# You can fill in your e-mail address and password.
client = Client(email_address="YOUR_EMAIL_ADDRESS", password="YOUR_PASSWORD")

Create Input Object

Create a Qamuy input object.

setting = qy.QamuyChemistryInput()

A Qamuy input object has the following structure.

QamuyChemistryInput:

• target_molecule

• output_chemical_properties

• post_hf_methods

• quantum_device

• mapping

• solver

• ansatz

• cost_function

• optimizer

Set Molecule

Here is an example of water molecule. You can set molecular geometrical structure, basis set, active space, the number of excited state and so on.

molecule = setting.target_molecule
molecule.basis = "6-31g"
molecule.multiplicity = 1
molecule.num_excited_states = 0
molecule.cas = qy.cas(2, 2)

# Set the geometry of water molecule
atoms = ["H", "O", "H"]
coords = [[1., 0., 0.], [0., 0., 0.], [0., 1., 0.]]
init_geom = qy.molecule_geometry(atoms, coords)  # initial geometry of H2O molecule
molecule.geometry = init_geom

Alternatively, you can set the geometry from your .xyz file in this way.

# Read your pre-made xyz file
with open("geometry_h2o.xyz") as f:
    xyz_str = f.read()
print(xyz_str)

3
Initial geometry for water molecule
H   1.00000 0.00000 0.00000
O   0.00000 0.00000 0.00000
H   0.00000 1.00000 0.00000

# Set the geometry of water molecule from xyz string
init_geom = qy.molecule_geometry_from_xyz(xyz_str)  # initial geometry of H2O molecule from xyz string
molecule.geometry = init_geom

Algorithm

# Solver
setting.solver.type = "VQE"
# or "VQD", "SSVQE", "MCVQE"

# Anasatz
setting.ansatz.type = "SYMMETRY_PRESERVING"
setting.ansatz.is_state_real = True
setting.ansatz.reference_state = "RHF"
setting.ansatz.depth = 4
setting.ansatz.use_random_initial_guess = True
# or "HARDWARE_EFFICIENT", "UCCSD", ...

# Optimizer
setting.optimizer.type = "BFGS"
setting.optimizer.ftol = 1e-08
setting.optimizer.gtol = 1e-08
setting.optimizer.max_iter = 1000

# Cost function
setting.cost_function.type = "SIMPLE"

# Mapping
setting.mapping.type = "JORDAN_WIGNER"
# or "BRAVYI_KITAEV", ...

# Device
setting.quantum_device.type = "EXACT_SIMULATOR"

# Classical algorithm
setting.post_hf_methods.append(
    qy.PostHFMethod(type="CASCI")
)

Options for geometry optimization

geom_opt = setting.geometry_optimization
geom_opt.target_state = 0
geom_opt.gradientmax = 1e-6
geom_opt.gradientrms = 1e-6
geom_opt.stepmax = 2e-2
geom_opt.steprms = 1.5e-2

Properties to be evaluated

setting.output_chemical_properties.append(
    qy.output_chemical_property(
        target="vibrational_analysis", states=[0], dx=0.001,
        type="HAMILTONIAN_ANALYTICAL"
        # or
        # type="HAMILTONIAN_NUMERICAL"
    )
)

Run Qamuy

job = client.submit(setting)
results = client.wait_and_get_job_results([job])
output = results[0].output

[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 4 concurrent workers.
[Parallel(n_jobs=-1)]: Done   1 tasks      | elapsed:  1.2min
[Parallel(n_jobs=-1)]: Done   1 out of   1 | elapsed:  1.2min finished

Output Geometry Optimization Result

# retrieve geometry optimization result
geo_result = output.geometry_optimization_result

# nuclear coordinates after geometry optimization
opt_geom = qy.get_optimized_geometry(geo_result)
list(opt_geom.coordinates)

[[1.0145534314938425, -0.09617455853183432, 2.7049645685439506e-16],
 [0.08162112703799226, 0.08162112703799221, 2.5630543206324393e-17],
 [-0.09617455853183349, 1.0145534314938425, 1.7891491676035803e-17]]

# molecular geometry can be output as xyz string
xyz_str = qy.geometry_to_xyz(opt_geom)
print(xyz_str)

3

H           1.01455       -0.09617        0.00000
O           0.08162        0.08162        0.00000
H          -0.09617        1.01455        0.00000

# save xyz string as an xyz file
with open("geometry_h2o_opt.xyz", "w") as f:
    f.write(xyz_str)

# history of geometry optimization
energy_hist = geo_result.energy_hist
coords_hist = geo_result.coordinates_history

from pprint import pprint
for step, (eng, coords)  in enumerate(zip(energy_hist, coords_hist)):
    print("\nIteration: ", step)
    print("Energy", eng)
    print("Coordinates: ")
    pprint(list(coords))

Iteration:  0
Energy -75.97070429288176
Coordinates:
[[1.0, -6.661338147750939e-16, 2.0030704868015217e-16],
 [0.0, 0.0, 0.0],
 [6.661338147750939e-16, 1.0, 1.1371144305660282e-16]]

Iteration:  1
Energy -75.9851030503343
Coordinates:
[[1.0200722426672586, -0.07997248255539341, 2.7269202750807856e-16],
 [0.05990023988813371, 0.05990023988813357, 1.8809782984341534e-17],
 [-0.07997248255539086, 1.020072242667259, 2.2516681244334993e-17]]

Iteration:  2
Energy -75.98631072164697
Coordinates:
[[1.0103552652670076, -0.09207446711642214, 1.919113607796106e-16],
 [0.08171920184941374, 0.08171920184941354, 2.5661340510684308e-17],
 [-0.09207446711641959, 1.0103552652670078, 9.644579044646018e-17]]

Iteration:  3
Energy -75.98637118971763
Coordinates:
[[1.0147584498403261, -0.09616001126680152, 1.923287765652749e-16],
 [0.08140156142647492, 0.08140156142647476, 2.5561595544158442e-17],
 [-0.09616001126679935, 1.0147584498403264, 9.612811962732183e-17]]

Iteration:  4
Energy -75.98637148253249
Coordinates:
[[1.0145560063708294, -0.09616803644161875, 2.7049743631904296e-16],
 [0.08161203007078924, 0.0816120300707891, 2.5627686590403906e-17],
 [-0.09616803644161714, 1.0145560063708294, 1.7893368827308335e-17]]

Iteration:  5
Energy -75.98637148288499
Coordinates:
[[1.0145534314938425, -0.09617455853183432, 2.7049645685439506e-16],
 [0.08162112703799226, 0.08162112703799221, 2.5630543206324393e-17],
 [-0.09617455853183349, 1.0145534314938425, 1.7891491676035803e-17]]

Output Vibration Analysis Result

# quantum algorithm results
q_result = output.molecule_result.quantum_device_result

# classical algorithm results
c_result = output.molecule_result.post_hf_results[0]

# VQE log
vqe_log = q_result.vqe_log

# Output energy
print(f"VQE  : {qy.get_evaluated_property(q_result, 'energy').values}")
print(f"CASCI: {qy.get_evaluated_property(c_result, 'energy').values}")

VQE  : [{'value': -75.98551775662527, 'state': 0, 'sample_std': 0.0}]
CASCI: [{'value': -75.98551775662887, 'state': 0, 'sample_std': 0.0}]

# vibration analysis results
vib_vqe = qy.get_evaluated_property(q_result, "vibrational_analysis")
vib_casci = qy.get_evaluated_property(c_result, "vibrational_analysis")

eng_vqe = [mode.frequency for mode in vib_vqe.values[0].value]
eng_casci = [mode.frequency for mode in vib_casci.values[0].value]

Here, the frequency corresponds to the energy (eigen frequencies) of each vibration mode, and is stored in real and imag parts (unit: cm^-1).

# Compare VQE vs CASCI
print("# \t VQE \t\t CASCI")
for i, (e_vqe, e_casci) in enumerate(zip(eng_vqe, eng_casci)):
    print(i, f"\t{e_vqe:.2f}", f"\t{e_casci:.2f}")

#        VQE             CASCI
0       0.00+23.46j     0.00+22.56j
1       0.00+7.43j      0.00+11.40j
2       0.00+0.01j      0.00+6.76j
3       0.00+0.00j      0.00+0.70j
4       0.01+0.00j      7.81+0.00j
5       12.17+0.00j     13.40+0.00j
6       1735.44+0.00j   1735.42+0.00j
7       3985.95+0.00j   3986.05+0.00j
8       4142.33+0.00j   4142.30+0.00j

Bonus: Visualization of molecules using py3Dmol

Install libraries necessary for visualization

Run the following cell, depending on whether you are using Jupyter Notebook or Jupyter Lab.

Once you’ve installed this, you don’t have to do it next time.

# for Jupyter Notebook users
!pip install py3Dmol  # viewer library for Jupyter Notebook

# for Jupyter Lab users
!pip install py3Dmol  # viewer library for Jupyter Notebook
!jupyter labextension install jupyterlab_3dmol  # extension for Jupyter Lab

Definition of functions for visualization

import py3Dmol

def visualize_mole(geometry, vib_modes=None):
    xyz_str = qy.geometry_to_xyz(geometry, vib_modes)
    if vib_modes:
        _view3d_vib(xyz_str)
    else:
        _view3d(xyz_str)


def _view3d(xyz_str):
    view = py3Dmol.view(width=400, height=400)
    view.addModel(xyz_str, "xyz")
    view.setStyle({'stick': {}, 'sphere': {'scale': .30}}, viewer=(0,0))
    view.setBackgroundColor('#ebf4fb', viewer=(0,0))
    view.show()


def _view3d_vib(xyz_str):
    view = py3Dmol.view(width=400,height=400)
    view.addModel(xyz_str,'xyz',{'vibrate': {'frames':10,'amplitude':.5}})
    view.setStyle({'stick': {}, 'sphere': {'scale': .30}}, viewer=(0,0))
    view.setBackgroundColor('#ebf4fb', viewer=(0,0))
    view.animate({'loop': 'backAndForth'})
    view.zoomTo()
    view.show()

Visualize molecules

# Before geometry optimization
visualize_mole(init_geom)

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

# After geometry optimization
visualize_mole(opt_geom)

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

Visualize vibration modes

n_vib_modes = 3 * len(atoms) - 6  # number of vibrational freedom of non-linear molecule. (3 * len(atoms) - 5 for linear molecule)
n_vib_modes

3

# eigenenergies for each vibration modes obtained with VQE
vib_freqs = [mode.frequency.real for mode in vib_vqe.values[0].value][3*len(atoms)-n_vib_modes:]
vib_freqs

[1735.4405461256251, 3985.951418866509, 4142.3256928356295]

# displacements for each vibration modes obtained with VQE
vib_modes = [mode.normal_mode for mode in vib_vqe.values[0].value][3*len(atoms)-n_vib_modes:]
vib_modes

[[0.14333246165223532, 0.6594568668357705, 4.813447332336657e-13, -0.05057888887956966, -0.05057888888464554, -9.772242324382186e-14, 0.6594568668446047, 0.14333246164509725, 4.819297897141424e-13],
 [0.6691284649282472, -0.18273953567314644, -1.5547957902642235e-12, -0.03064441781836253, -0.030644417817726224, 7.09932929635145e-14, -0.18273953566877674, 0.6691284649098128, -1.5545532726539354e-12],
 [-0.6638342540023767, 0.1264986979872842, 4.55617045026171e-13, 0.04979408810239, -0.04979408810401285, -7.419612795371809e-20, -0.12649869799034297, 0.6638342540202272, -4.557616184912729e-13]]

for idx, (freq, mode) in enumerate(zip(vib_freqs, vib_modes)):
    print(f"\nmode {idx+1}")
    print(f"frequency: {freq: 5.2f} [cm-1]")
    visualize_mole(opt_geom, mode)

mode 1
frequency:  1735.44 [cm-1]

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

mode 2
frequency:  3985.95 [cm-1]

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

mode 3
frequency:  4142.33 [cm-1]

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