Band structure calculation 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/

Import Qamuy#

[1]:
import qamuy.chemistry as qy
from qamuy.client import Client

from copy import deepcopy
import numpy as np
import matplotlib.pyplot as plt

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

input instance generation#

[2]:
input = qy.QamuyChemistryInput()
ps = input.target_periodic_system

unit cell#

[3]:
# see also https://pyscf.org/user/pbc/gto.html
BOHR = 0.529 # unit: angstrom
unitcell_x = 4. * BOHR
unitcell_yz = 1.
atoms = ["H", "H"]
coords = [[0., 0., 0.],[BOHR, 0., 0.]]
lattice_vec = [[unitcell_x, 0., 0.], [0., unitcell_yz, 0.], [0., 0., unitcell_yz]]
kpt_grid_shape = [1, 1, 1] # number of k-points in each direction.
dimension = 1

ps.geometry = qy.periodic_system_geometry(atoms, coords, lattice_vec, kpt_grid_shape)
ps.geometry.dimension = dimension
[4]:
from PIL import Image
filename = "figs/hydrogen_chain.png"
im = Image.open(filename)
im
[4]:
../_images/_usecase_04_band_calculation_10_0.png

basis, number of excited states, and multiplicity#

[5]:
ps.basis = "sto-3g"
ps.num_excited_states = 0
ps.multiplicity = 1  # required (from ver. 0.18.1)

# CAS
ps.cas = qy.cas(2, 2)

solver, ansatz, optimizer, cost function, and quantum device#

[6]:
input.solver.type = "VQE"

input.ansatz.type = "UCCD"
input.ansatz.is_state_real = True
input.ansatz.reference_state = "RHF"
input.ansatz.trotter_steps = 4
input.ansatz.initial_parameter = [0.0]
input.ansatz.init_param_random_seed = 1

input.optimizer.type = "BFGS"
input.optimizer.gtol = 1e-6

input.cost_function.type = "SIMPLE"

# exact simulator (state vector)
input.quantum_device.type = "EXACT_SIMULATOR"

properties to be evaluated#

[7]:
properties = input.output_chemical_properties
properties.append(qy.output_chemical_property(target="band_structure", algorithm="QSE"))
properties.append(qy.output_chemical_property(target="band_structure", algorithm="QEOM"))

run#

[8]:
kx_list = [round(0.05*i,2) for i in range(20)]

inputs = []
jobs = []
outputs = []

for kx in kx_list:
    temp_input = deepcopy(input)
    ps = temp_input.target_periodic_system
    ps.scaled_center = [kx, 0., 0.]
    inputs.append(temp_input)

for input_ in inputs:
    jobs.append(client.submit(input_))

results = client.wait_and_get_job_results(jobs)
outputs = [result.output for result in results]
[Parallel(n_jobs=-1)]: Using backend ThreadingBackend with 8 concurrent workers.
[Parallel(n_jobs=-1)]: Done   2 tasks      | elapsed:  6.4min
[Parallel(n_jobs=-1)]: Done   8 out of  20 | elapsed:  6.5min remaining:  9.8min
[Parallel(n_jobs=-1)]: Done  11 out of  20 | elapsed:  6.5min remaining:  5.3min
[Parallel(n_jobs=-1)]: Done  14 out of  20 | elapsed:  6.5min remaining:  2.8min
[Parallel(n_jobs=-1)]: Done  17 out of  20 | elapsed:  6.7min remaining:  1.2min
[Parallel(n_jobs=-1)]: Done  20 out of  20 | elapsed:  6.8min remaining:    0.0s
[Parallel(n_jobs=-1)]: Done  20 out of  20 | elapsed:  6.8min finished

extract band energies#

[9]:
kx_list = []
conduction_qse = []
conduction_qeom = []
conduction_hf = []
valence_qse = []
valence_qeom = []
valence_hf = []

for output in outputs:
    q_result = output.molecule_result.quantum_device_result
    gs_energy = q_result.evaluated_properties[0].energy.values[0].value
    band_energies_qse = qy.get_evaluated_properties(q_result, "band_structure")[0].values
    band_energies_qeom = qy.get_evaluated_properties(q_result, "band_structure")[1].values
    hf_result = output.molecule_result.hf_result

    kx_list.append(band_energies_qse[0].kpt[0])

    conduction = []
    valence = []
    for band_dict in band_energies_qse:
        if band_dict.band == "conduction band":
            conduction.append(band_dict.value)
        elif band_dict.band == "valence band":
            valence.append(band_dict.value)
    idx_conduction = np.abs(np.asarray(conduction) - gs_energy).argmin()
    idx_valence = np.abs(np.asarray(valence) - gs_energy).argmin()
    conduction_qse.append(conduction[idx_conduction])
    valence_qse.append(valence[idx_valence])

    conduction = []
    valence = []
    for band_dict in band_energies_qeom:
        if band_dict.band == "conduction band":
            conduction.append(band_dict.value)
        elif band_dict.band == "valence band":
            valence.append(band_dict.value)
    idx_conduction = np.abs(np.asarray(conduction) - gs_energy).argmin()
    idx_valence = np.abs(np.asarray(valence) - gs_energy).argmin()
    conduction_qeom.append(conduction[idx_conduction])
    valence_qeom.append(valence[idx_valence])

    conduction_hf.append(hf_result.mo_energy[1])
    valence_hf.append(hf_result.mo_energy[0])

plot#

[10]:
fig = plt.figure(figsize=(8, 6))
ax = fig.add_subplot(111)

ax.plot(kx_list, conduction_qse, color="black", marker="o", label="QSE")
ax.plot(kx_list, conduction_qeom, color="red", marker=".", label="qEOM")
ax.plot(kx_list, conduction_hf, ls="dashed", marker="^", color="green", label="HF")
ax.plot(kx_list, valence_qse, marker="o", color="black")
ax.plot(kx_list, valence_qeom, marker=".", color="red")
ax.plot(kx_list, valence_hf, marker="^", ls="dashed", color="green")

ax.grid(True, which='both', axis='both')
ax.set_xlabel("k")
ax.set_ylabel("Energy [Hartree]")
ax.set_title("Band structure of hydrogen chain")
ax.legend(fontsize=18)
for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] + ax.get_xticklabels() + ax.get_yticklabels()):
    item.set_fontsize(16)
fig.tight_layout()
../_images/_usecase_04_band_calculation_22_0.png