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#

[1]:
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.

[2]:
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.

[3]:
molecule = setting.target_molecule
molecule.basis = "6-31g"
molecule.multiplicity = 1
molecule.num_excited_states = 0
molecule.cas = qy.cas(2, 2)
[4]:
# 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.

[5]:
# 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
[6]:
# 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#

[7]:
# 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#

[8]:
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#

[9]:
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#

[10]:
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#

[11]:
# retrieve geometry optimization result
geo_result = output.geometry_optimization_result
[12]:
# nuclear coordinates after geometry optimization
opt_geom = qy.get_optimized_geometry(geo_result)
list(opt_geom.coordinates)
[12]:
[[1.0145534314938425, -0.09617455853183432, 2.7049645685439506e-16],
 [0.08162112703799226, 0.08162112703799221, 2.5630543206324393e-17],
 [-0.09617455853183349, 1.0145534314938425, 1.7891491676035803e-17]]
[13]:
# 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
[14]:
# save xyz string as an xyz file
with open("geometry_h2o_opt.xyz", "w") as f:
    f.write(xyz_str)
[15]:
# 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#

[16]:
# 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
[17]:
# 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}]
[18]:
# vibration analysis results
vib_vqe = qy.get_evaluated_property(q_result, "vibrational_analysis")
vib_casci = qy.get_evaluated_property(c_result, "vibrational_analysis")
[19]:
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).

[20]:
# 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#

[21]:
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#

[22]:
# 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

[23]:
# 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#

[24]:
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
[24]:
3
[25]:
# 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
[25]:
[1735.4405461256251, 3985.951418866509, 4142.3256928356295]
[26]:
# 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
[26]:
[[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]]
[27]:
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