Examples#
This page contains a collection of example calculations which can be run using the SCCE engine. The examples presented here will use Python code for the readout and plotting of the data contained in the result files. However, Python is not required in general - the data can be read from the result files using any hdf5 or json compatible program or programming language, depending on the chosen output type.
Simple Examples#
Here we give two examples of one-dimensional calculations that can be run using the SCCE engine. These examples do not model any complicated physical effects, and the results could just as well be obtained using simpler numerical methods. This section is intended only as an introduction to constructing the unit cell and system yaml files, as well as using the data contained in the result files.
In this section we will always set system_size = [1, 1, 1].
Only the DMRG cluster solver will be used in the calculations presented here, and we thus skip the self-consistent SCCE loop. For small to medium sized one-dimensional systems, the DMRG solver is very fast - both faster and more accurate than using the full SCCE loop.
1D spinless Fermions#
One of the simplest possible example which can be run on the SCCE platform is a one-dimensional system of non-interacting, spinless fermions, described by the Hamiltonian
We choose hard-wall boundary conditions, with system size \(M=40\). The following python code defines, in a dictionary, the unit cell, with a hopping strength of \(t=-1\), and the system, which is given to the LoopSCCE class as an argument. The loop function then launches the computation.
import scce.loop
from scce import save_result,LoopSCCE
# Create a unit-cell
unitcell = {"atoms": [['0', 'up', [0.0, 0, 0], 0],
],
"bonds": [['0', '0', [1, 0, 0], -1, 0],
],
"lattice_vectors": [[1, 0, 0]]
}
# Create a system
system = {"cluster_size": [40, 1, 1], # measured in unit cells
"system_size": [1, 1, 1], # measured in clusters
"system_boundary_condition": "hard-wall",
"N":1,
"site_type": 'spinless_fermions',
}
indict = {'unitcell': unitcell, 'system': system}
loop = LoopSCCE(indict,
DMRG_dirResults="./results/DMRG_runs",
DMRG_dirScratch="./results/SCCE")
loop.NumStates = 4
loop(maxIter=20)
save_result(loop, filename='1D_fermions', filedir='./results', overwrite=True)
Here we've chosen to have just one fermion in the system, \(N=1\), and to calculate four eigenstates. The combination \(N=1\) with hard-wall boundary conditions allows us to easily verify that the result files contain the expected sinusoidal solutions for the fermion density of the eigenstates. For more interesting density distributions we recommend changing the particle number to \(N=5\).
The fermion densities of the eigenstates in the 1D_fermions created by the SCCE engine can be plotted with the following Python code:
import h5py
import numpy as np
from matplotlib import pyplot as plt
RESULT_FILE = 'results/1D_fermions.h5'
with h5py.File(RESULT_FILE,'r') as result:
# ~ create figuure and axes
fig_density, ax_density = plt.subplots(1, 1, figsize=(9.5,5.9))
# ~ get number of unit cells in system
xx=np.arange(0,result["system_cluster_size"][()])
# ~ get diagonal of reduced density matrix
for M in range(4):
yy=np.diagonal(result["result_state_{}_1rdm_real".format(M)][()])
# ~ plot density vs number of sites
ax_density.plot(xx,yy,label='state {}'.format(M), lw=2)
# ~ create legend
ax_density.legend(fontsize=12)
ax_density.set_title('Fermion density N=1', fontsize=12)
ax_density.set_xlabel('position', fontsize=12)
ax_density.set_ylabel('Fermion density', fontsize=12)
# ~ show figure
plt.show()
The Python code will produce the following plot of the fermion density:
Plot of the 1D fermion density for a 40-site system with one particle and hard-wall boundary conditions#
1D Hubbard Model#
To show an example with spin, let us use the one-dimensional Hubbard model with hard-wall boundary conditions:
For this example we choose \(M=20\), six electrons in the system, \(N=6\), and a strongly attractive interaction, \(U=-20\).
This example can be run with the following code.
import scce.loop
from scce import save_result,LoopSCCE
# Create a unit-cell
# use a hopping of 1 for spin up fermions, and 0.9 for down spin fermions to avoid degeneracies
unitcell = {"atoms": [['0', 'up', [0.0, 0, 0], -20.],
],
"bonds": [['0', '0', [1, 0, 0], [-1, -0.9] , 0],
],
"lattice_vectors": [[1, 0, 0]]
}
# Create a system
system = {"cluster_size": [20, 1, 1], # measured in unit cells
"system_size": [1, 1, 1], # measured in clusters
"system_boundary_condition": "hard-wall",
"N":6,
"site_type": 'spinful_fermions',
}
indict = {'unitcell': unitcell, 'system': system}
loop = LoopSCCE(indict,
DMRG_dirResults="./results/DMRG_runs",
DMRG_dirScratch="./results/SCCE")
loop.NumStates = 4
loop(maxIter=20)
save_result(loop, filename='hubbard_model', filedir='./results', overwrite=True)
We plot the spin resolved density and the double occupancy on all cluster sites with Python:
import h5py
import numpy as np
from matplotlib import pyplot as plt
# ~ h5 file name
RESULT_FILE="./results/hubbard_model.h5"
with h5py.File(RESULT_FILE,'r') as result:
# ~ create figuure and axes
fig_density, ax_density = plt.subplots(1, 1, figsize=(9.5,5.9))
# ~ get number of unit cells in system
xx=np.arange(0,result["system_cluster_size"][()])
# ~ get diagonal of reduced density matrix
for M in range(4):
yy=np.diagonal(result["result_state_{}_1rdm_up_up_real".format(M)][()])
# ~ plot density vs number of sites
ax_density.plot(xx,yy,label='state {}'.format(M), lw=9 - 2 * M)
# ~ create legend
ax_density.legend(fontsize=12)
ax_density.set_title('Spin-resolved density: Density of up electrons, U=-20', fontsize=12)
ax_density.set_xlabel('position', fontsize=12)
ax_density.set_ylabel('density', fontsize=12)
# ~ show figure
plt.show()
# ~ create figuure and axes
fig_density, ax_density = plt.subplots(1, 1, figsize=(9.5,5.9))
# ~ get number of unit cells in system
xx=np.arange(0,result["system_cluster_size"][()])
# ~ get double occupancy correlator
for M in range(4):
yy=np.diagonal(result["result_state_{}_double_occupancy_correlator".format(M)][()])
# ~ plot density vs number of sites
ax_density.plot(xx,yy,label='state {}'.format(M), lw=9 - 2 * M)
# ~ create legend
ax_density.legend(fontsize=12)
ax_density.set_title('Double occupancy, U=-20', fontsize=12)
ax_density.set_xlabel('position', fontsize=12)
ax_density.set_ylabel('density', fontsize=12)
# ~ show figure
plt.show()
Comparing the two plots, we see that the density of spin up electrons is approximately equal to the double occupancy. In the strongly attractive system, all electrons are paired up. To explore the SCCE platform further, we recommend experimenting with different values for U. For example, choosing repulsive interactions with \(U=2\) will almost completely suppress the double occupancy.
Plot of the spin-up density in the 20-site 1D Hubbard model with \(U=-20\)#
Plot of the double occupancy in the 20-site 1D Hubbard model with \(U=-20\)#
Spinless fermions with nearest-neighbor interaction#
Here we give an example of a use case of the effective non-interacting model that is generated by SCCE. For this we turn to a model of spinless fermions with nearest-neighbor interaction, which can be mapped on a quantum spin 1/2 chain, the XXZ model, via a Jordan-Wigner transformation. Its Hamiltonian reads:
For this example we choose a system of \(M=42\) lattice sites occupied by \(N=21\) spinless fermions. The nearest-neighbor interaction is varied as \(U=-2.0, \dots, 2.0\).
We then run a SCCE calculation for each value of \(U\). In this example, we use steps of \(\Delta U = 0.05\).
import scce.loop
from scce import save_result,LoopSCCE
import numpy as np
for nn_interaction in np.arange(-2.0, 2.0, 0.05):
# Create a unit-cell
unitcell = {"atoms": [['0', 'X', [0.0, 0, 0], 0],
],
"bonds": [['0', '0', [1, 0, 0], -1, float(nn_interaction)],
],
"lattice_vectors": [[1, 0, 0]]
}
# Create a system
system = {"cluster_size": [42, 1, 1], # measured in unit cells
"system_size": [1, 1, 1], # measured in clusters
"system_boundary_condition": ["periodic"],
"N":21,
"site_type": 'spinless_fermions',
}
indict = {'unitcell': unitcell, 'system': system}
loop = LoopSCCE(indict,
DMRG_dirResults="./results/DMRG_runs",
DMRG_dirScratch="./results/SCCE")
loop.NumStates = 1
loop(maxIter=20)
save_result(loop, filename="3_1_RESULT_FILE_{val:.3f}.h5".format(val=nn_interaction), filedir='./results', overwrite=True)
Then, we plot the results, using the following code.
import h5py
import numpy as np
import matplotlib.pyplot as plt
import os
import re
RESULT_DIR="./results/"
def extract_U(s):
return float(re.search(r"E_(.*)\.h5",s)[1])
list_of_vF = []
filenames = [x for x in os.listdir(RESULT_DIR) if x.startswith("3_1_RESULT_FILE")]
filenames.sort(key=extract_U)
for filename in filenames:
with h5py.File(RESULT_DIR +"/"+filename, 'r') as f:
U = extract_U(filename)
M = int(f["system_size"][()])
effective_hamiltonian = f['result_cluster_mean_field_hamiltonian_real'][()]
vF_num = -(effective_hamiltonian[M // 2, M //2 - 1] + effective_hamiltonian[M // 2 - 1, M // 2 - 2])
eta = 0.5 * np.arccos(U / -2.0)
vF_ana = np.pi * np.sin(2 * eta) / ( np.pi - 2 * eta)
list_of_vF.append([U, vF_num, vF_ana])
vF_data = np.array(list_of_vF)
fig, ax = plt.subplots(figsize=(9.5, 5.9))
ax.plot(vF_data[:, 0], vF_data[:, 1], 'o', label='inverse Hartree-Fock')
ax.plot(vF_data[:, 0], vF_data[:, 2], '-', label='Bethe Ansatz')
legend = ax.legend(loc='upper left')
ax.set_title("Renormalized Fermi velocity $v_F$")
ax.set_xlabel("Interaction strength $U/t$")
ax.set_ylabel("Fermi velocity $v_F$")
plt.show()
Plot of the renormalized Fermi velocity as a function of the interaction strength \(U\) with comparison to the analytical result from Bethe ansatz#
1D Heisenberg model (frozen magnon calculation)#
In this example we will study the 1D isotropic Heisenberg model using SCCE's feature to simulate spin systems. The Heisenberg model is a prototypical spin model and we will investigate its low energy excitations, the so-called magnons.
The Hamiltonian of the isotropic Heisenberg model is given by
where we consider periodic boundary conditions, i.e., site \(j=M+1\) is identical to site \(j=1\). In the "frozen magnon approach" we apply an external magnetic field \(\vec{B}_j\) in order to force them into a spin-wave configuration. The external magnetic field is given by
This constraining magnetic field induces as static spin-wave configuration
characterized by the wave vector \(q\), identical to the wave vector of the magnetic field, and an opening angle \(\phi\), which, in general, is different from the opening angle of the magnetic field.
We start by setting up the system input file:
system:
cluster_size: [1, 1, 1]
system_size: [1, 1, 1]
boundary_condition: 'periodic'
site_type: 'spins'
mod_Sz: 2
where we use the input option site_type to specify that we want to run a simulation of a
spin system and the input option mod_Sz = 2 to tell the solver that \(S_z\) is no longer
a good quantum number due to the presence of the constraining magnetic field. Note that
\(S_z \bmod 2\) is still a good quantum number. The above yaml input must be stored in a
file named magnon_system.yaml for the following Python scripts to work.
In order to describe the frozen-magnon state we define \(M=16\) spins in our unit cell, each with its constraining magnetic field. Instead of using the web-interface to define unit cells for each wave vector \(q\) commensurate with our M-spin unit cell (super cell), we employ the following python script to generate a yaml input files for the unit cells which than can be loaded in the calculation step.
# Generic imports
import io
import numpy as np
import yaml
# set system size
nSpins = 16
# set number of Brioullin zones to sample magnon wave vector
nBZ = 1
# set isotropic nearest neighbor FM coupling
J = 1.
# set strength of the constraining magnetic field
B = 1.
# set opening angle theta for contraining magnetic field
# this is the *small parameter* in the classical spin-wave theory
# ... feel free to explore differences to classical spin-wave theory
# when choosing the openening angle theta bigger for the frozen magnon calculation
theta = 1.e-2 * np.pi
# set isotropic Heisenberg interaction between nearest neighbor spins
bonds = [
{
"id_from": i,
"id_to": (i + 1) % nSpins,
"translation": [(i + 1) // nSpins, 0, 0],
"J": J
}
for i in range(nSpins)
]
# set longitudinal and transverse amplitudes of magnetic field
ct = np.cos(theta)
st = np.sin(theta)
# generate set of commensurate magnon wave vectors
Q = nBZ * 2 * (np.linspace(0., 1., num=(nBZ * nSpins + 1), endpoint=True) - .5)
for q in Q:
# construct unit cell containing all spins (super cell) and constraining magnetic field
atoms = [
{
"id": i,
"name": f"spin {i+1}",
"position": [i, 0, 0],
"Bx": float(B * st * np.cos(q * i * np.pi)),
"By": float(B * st * np.sin(q * i * np.pi)),
"Bz": float(B * ct)
}
for i in range(nSpins)
]
# set unitcell
unitcell = {'atoms': atoms,
'lattice_vectors': [[1, 0, 0], [0, 0, 0], [0, 0, 0]],
'bonds': bonds,
}
# Write unit cell YAML file for frozen magnon calculation
with io.open(f"heisenberg1D_FM_q{q:+7.4f}.yaml", 'w', encoding='utf8') as file:
yaml.dump({'unitcell': unitcell}, file, default_flow_style=False, indent=2)
After we have created the input files for the magnon wave vectors, we can start the calculation
by combining each stored unit-cell input generated in the previous Python script with the common
system input directly created as yaml file. This is accomplished by the following Python script
# Generic imports
import io
import yaml
from glob import glob
from pprint import pprint
# HQS imports
from lattice_validator import Validator as LV
from lattice_builder import Builder as LB
import scce.loop
from scce import save_result, LoopSCCE
# set 'results', 'output' and 'scratch' directories for SCCE / DMRG solver
DMRG_results_dir = "./results/DMRG_runs"
DMRG_scratch_dir = "/tmp/DMRG_scratch"
SCCE_output_dir = "./results"
# loading system configuration (... same for all wave vectors)
with io.open("magnon_system.yaml", "r", encoding='utf8') as file:
config_system = yaml.load(file, Loader=yaml.SafeLoader)
# getting list of all unit cells for various wave vectors
magnon_unitcells = glob("heisenberg1D_FM_q*.*.yaml")
for input_file in magnon_unitcells:
# extract magnon wave vector from file name
Q = float(input_file.split("heisenberg1D_FM_q")[1].split(".yaml")[0])
# open input file for unit cell at current magnon wave vector
with io.open(f"{input_file}", 'r', encoding='utf8') as file:
config_uc = yaml.load(file, Loader=yaml.SafeLoader)
# merging unit cell and system into full configuration
config = config_uc | config_system
# validate input
lv = LV()
valid, normalized_configuration = lv.validateConfiguration(config)
# exit in case of error
if not valid:
# pretty print error message
pprint(normalized_configuration)
exit(1)
# initialize lattice builder
lb = LB(normalized_configuration)
print(f'lb.is_complex: {lb.is_complex}')
print(f'lb.site_type: {lb.site_type}')
print(f'lb.is_spin_matrix: {lb.is_spin_matrix}')
# initialize SCCE loop
loop = LoopSCCE(
normalized_configuration,
DMRG_dirResults=DMRG_results_dir,
DMRG_dirScratch=DMRG_scratch_dir,
)
# run SCCE loop
loop.clear()
loop(maxIter=20, skipCalc=False)
# store results from SCCE run
save_result(loop, filename=f"1D_magnon_Q{Q:+7.4f}", filedir=SCCE_output_dir, overwrite=True)
We can extract the magnon dispersion from the energies computed at each wave vector \(q\). Note that the energy \(E(q)\), provided by SCCE, contains the energy contribution due to the coupling to the constraining magnetic field. Since SCCE provides the expectation value of the spin operator, \(\vec{S}_j = \langle \hat{\vec{S}}_j \rangle\), we can extract the intrinsic energy
Furthermore, we can extract the opening angle \(\phi\) from \(\vec{S}_j\), which is needed to extract the magnon dispersion, defined by
where \(s = 1/2\) is the spin quantum number. The results are shown in the following two plots:
Plot of the magnon dispersion
2) Plot comparing the opening angles of the spin magnetization, \(\phi\), to the opening angle of the constraining magnetic field, \(\theta\).
Below we provide a python script performing the steps outlined above:
# Generic imports
import io
import h5py
from glob import glob
import yaml
import numpy as np
from matplotlib import pylab as plt
# set spin quantum number ... important for comparison to classical spin-wave theory
s = 1 / 2
# getting list of all wave vectors
magnon_unitcells = glob("heisenberg1D_FM_q*.*.yaml")
q_list: list[float] = []
for input_file in magnon_unitcells:
# extract magnon wave vector from file name
q_list.append(float(input_file.split("heisenberg1D_FM_q")[1].split(".yaml")[0]))
# sort magnon wave vectors in ascending order
q_list.sort()
q_mag = np.zeros(len(q_list), dtype=float)
# initialize array for energies
E = np.zeros_like(q_mag)
# initialize array for opening angles of spin magnetization
phi = np.zeros_like(q_mag)
# initialize array for opening angles of constraining magnetic field
theta = np.zeros_like(q_mag)
for q_idx, q in enumerate(q_list):
q_mag[q_idx] = q
with h5py.File(f"./results/1D_magnon_Q{q:+7.4f}.h5", "r") as result_file:
# extract number of spins from result file
nSpins = int(result_file["unitcell_number_atoms"][()])
# initialize array holding spin magnetization
S = np.zeros((nSpins, 3), dtype=float)
# extract spin magnetization from output file
S[:, 0] = result_file["result_state_0_Sx"]
S[:, 1] = result_file["result_state_0_Sy"]
S[:, 2] = result_file["result_state_0_Sz"]
# determine opening angle of spin magnetization
phi[q_idx] = np.arctan(np.sqrt(S[:, 0]**2 + S[:, 1]**2) / S[:, 2]).mean()
# initialize array for constraining magnetic field
B = np.zeros((nSpins, 3), dtype=float)
# extract constraining magnetic field from input file
input_file = f"./heisenberg1D_FM_q{q:+7.4f}.yaml"
with io.open(f"{input_file}", 'r', encoding='utf8') as input_file:
config_uc = yaml.load(input_file, Loader=yaml.SafeLoader)
for i, atom in enumerate(config_uc["unitcell"]["atoms"]):
for j, B_comp in enumerate(["Bx", "By", "Bz"]):
B[i, j] = atom[B_comp]
# determine opening angle of constraining magnetic field
theta[q_idx] = np.arctan(np.sqrt(B[:, 0]**2 + B[:, 1]**2) / B[:, 2]).mean()
# extract energy and remove energy due to constraining field
e_total = float(result_file["system_cluster_energy"][0])
E[q_idx] = e_total + (S * B).sum()
if abs(q) < 1.e-3:
# store reference energy at zero wave vector
E0 = E[q_idx]
# set results for wave vectors
w_mag = (E - E0) / (np.sin(phi)**2 * s * s * nSpins)
# generate classical spin-wave result for magnon dispersion
x = np.linspace(-1., 1., num=101, endpoint=True)
wx = 1. - np.cos(x * np.pi)
# plot magnon dispersion
plt.figure()
plt.title(f"1D Heisenberg model ({nSpins} spins) - magnon dispersion")
plt.plot(q_mag, w_mag, lw=3, ls="-", color="red", label=r"$\omega(q)$ SCCE simulation")
plt.plot(x, wx, lw=2, ls="--", color="black",
label=r"$\omega(q) = 2 J [1 - \cos(q)]$ (classical spin-wave theory)")
plt.grid()
plt.xlim((-1., 1.))
plt.ylim((-0.1, 2.3))
plt.xlabel(r"wave vector $q$ [$\pi / a$]")
plt.ylabel(r"magnon dispersion [$J$]")
plt.legend(loc=0)
plt.savefig("1DHeisenbergModelDispersion.svg")
# plot opening angles fo constraining magnetic field and spin magnetization
plt.figure()
plt.title(f"1D Heisenberg model ({nSpins} spins) - opening angles")
plt.plot(q_mag, theta / np.pi, lw=3, ls="-", label=r"$\theta(q)$")
plt.plot(q_mag, phi / np.pi, lw=3, ls="-", label=r"$\phi(q)$")
plt.grid()
plt.xlim((-1., 1.))
plt.ylim((-.1e-2, 1.1e-2))
plt.xlabel(r"wave vector $q$ [$\pi / a$]")
plt.ylabel(r"opening angle [$\pi$]")
plt.legend(loc=0)
plt.savefig("1DHeisenbergModelAngles.svg")
Running the script above produces two svg plots showing the magnon dispersion (including a comparison to classical spin-wave theory) and a comparison of the opening angles of the magnetization compared to the opening angles of the constraining magnetic field.
Finally, we provide the result for 32 spins
