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

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

input = qy.QamuyChemistryInput()
ps = input.target_periodic_system

unit cell

# 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

from PIL import Image
filename = "figs/hydrogen_chain.png"
im = Image.open(filename)
im

basis, number of excited states, and multiplicity

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

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

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

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

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

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()