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