This tutorial can be downloaded link.
Advanced Tutorial 2: Detailed Analysis of QDET Many-Body States¶
In this tutorial, we show how the full diagonalization of effective Hamiltonian in QDET enables a detailed analysis of the excited states of the defect.
As we focus on the analysis in this tutorial, we download the input and output files from the QDET calculation. For more details about how to perform the underlying QDET calculations, see Basic Tutorial 5.
[1]:
%%bash
mkdir -p west.wfreq.save
wget -N -q https://west-code.org/doc/training/nv_diamond_63/wfreq.in
wget -N -q https://west-code.org/doc/training/nv_diamond_63/wfreq.json -O west.wfreq.save/wfreq.json
Solution of QDET for \(m_S=0\)¶
We diagonalize the QDET Hamiltonian as usual:
[2]:
from westpy.qdet import QDETResult
# construct object for effective Hamiltonian
effective_hamiltonian = QDETResult(filename='west.wfreq.save/wfreq.json')
# diagonalize Hamiltonian
solution = effective_hamiltonian.solve()
_ _ _____ _____ _____
| | | | ___/ ___|_ _|
| | | | |__ \ `--. | |_ __ _ _
| |/\| | __| `--. \ | | '_ \| | | |
\ /\ / |___/\__/ / | | |_) | |_| |
\/ \/\____/\____/ \_/ .__/ \__, |
| | __/ |
|_| |___/
WEST version : 5.5.0
Today : 2024-02-02 14:09:28.333370
-----------------------------------------------------
===============================================================
Building effective Hamiltonian...
nspin: 1
occupations: [[2. 2. 2. 2. 1. 1.]]
===============================================================
| diag[1RDM - 1RDM(GS)] | ||||||||
|---|---|---|---|---|---|---|---|---|
| E [eV] | char | 87 | 122 | 123 | 126 | 127 | 128 | |
| 0 | 0.000 | 3- | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
| 1 | 0.436 | 1- | -0.001 | -0.009 | -0.018 | -0.067 | 0.004 | 0.091 |
| 2 | 0.436 | 1- | -0.001 | -0.009 | -0.018 | -0.067 | 0.092 | 0.002 |
| 3 | 1.251 | 1- | -0.002 | -0.019 | -0.023 | -0.067 | 0.054 | 0.057 |
| 4 | 1.939 | 3- | -0.003 | -0.010 | -0.127 | -0.860 | 1.000 | 0.000 |
| 5 | 1.940 | 3- | -0.003 | -0.010 | -0.127 | -0.860 | 0.000 | 1.000 |
| 6 | 2.935 | 1- | -0.000 | -0.032 | -0.043 | -0.855 | 0.929 | 0.002 |
| 7 | 2.936 | 1- | -0.000 | -0.032 | -0.043 | -0.855 | 0.002 | 0.929 |
| 8 | 4.661 | 1- | -0.006 | -0.054 | -0.188 | -1.672 | 0.960 | 0.960 |
| 9 | 5.080 | 3- | -0.014 | -0.698 | -0.213 | -0.075 | 1.000 | 0.000 |
-----------------------------------------------------
The solution was obtained for 10 electrons (5 in each spin channel) and 6 orbitals:
[3]:
print(solution['nelec'])
print(solution['norb'])
(5, 5)
6
These numbers can be compared with the basis and occupation stored in the effective Hamiltonian.
[4]:
print(effective_hamiltonian.basis)
print(len(effective_hamiltonian.basis))
[ 87 122 123 126 127 128]
6
[5]:
import numpy as np
print(effective_hamiltonian.occupation)
print(np.sum(effective_hamiltonian.occupation))
[[2. 2. 2. 2. 1. 1.]]
10.0
This is the case of solutions with \(m_S=0\). We will later discuss the case with \(m_S=1\).
To check for the value of \(S^2\) and spin multiplicity, one can easily:
[6]:
import pyscf
fcisolver = pyscf.fci.direct_spin1.FCISolver()
for i, state in enumerate(solution['evcs']):
en = solution['evs'][i]
ss, mult = fcisolver.spin_square(state, norb=solution['norb'], nelec=solution['nelec'])
print(f"i={i}, energy [eV]={en:.3f}, <S^2>={ss:.3f}, mult={mult:.3f}")
i=0, energy [eV]=0.000, <S^2>=2.000, mult=3.000
i=1, energy [eV]=0.436, <S^2>=0.000, mult=1.000
i=2, energy [eV]=0.436, <S^2>=0.000, mult=1.000
i=3, energy [eV]=1.251, <S^2>=0.000, mult=1.000
i=4, energy [eV]=1.939, <S^2>=2.000, mult=3.000
i=5, energy [eV]=1.940, <S^2>=2.000, mult=3.000
i=6, energy [eV]=2.935, <S^2>=0.000, mult=1.000
i=7, energy [eV]=2.936, <S^2>=0.000, mult=1.000
i=8, energy [eV]=4.661, <S^2>=0.000, mult=1.000
i=9, energy [eV]=5.080, <S^2>=2.000, mult=3.000
The eigenstate of each many-body state is a matrix of size (number of Slater determinants in the spin-alpha channel, number of Slater determinants in the spin-beta channel), where each entry represents the contribution of that product of Slater determinants to the many-body states. We can determine the number of Slater determinants in each spin-channel as follows:
[7]:
from pyscf.fci import cistring
# print the number of possible Slater determinants in each spin channel
for ispin, spinlabel in enumerate(['alpha','beta']):
print(f'No. Slater determinants in spin-channel {spinlabel} = {len(cistring.make_strings(range(solution["norb"]),solution["nelec"][ispin]))}')
No. Slater determinants in spin-channel alpha = 6
No. Slater determinants in spin-channel beta = 6
We check that the shape of the eigenvector matrix is 6x6
[8]:
print(solution['evcs'][0].shape)
(6, 6)
We can now inspect the contribution of the various Slater determinants to a given excited many-body state.
[9]:
from westpy.qdet import visualize_correlated_state
for i, state in enumerate(solution['evcs']):
en = solution['evs'][i]
print(f'i= {i}, en [eV]={en}, ',visualize_correlated_state(state, solution['norb'], solution['nelec'], cutoff=10 **(-3)))
print( )
i= 0, en [eV]=0.0, -0.707|011111>|101111>+0.707|101111>|011111>
i= 1, en [eV]=0.4359576505738596, +0.313|011111>|011111>+0.596|011111>|101111>-0.035|011111>|110111>-0.018|011111>|111011>+0.012|011111>|111101>+0.004|011111>|111110>+0.596|101111>|011111>-0.312|101111>|101111>+0.180|101111>|110111>+0.093|101111>|111011>-0.064|101111>|111101>-0.021|101111>|111110>-0.035|110111>|011111>+0.180|110111>|101111>-0.018|111011>|011111>+0.093|111011>|101111>+0.012|111101>|011111>-0.064|111101>|101111>+0.004|111110>|011111>-0.021|111110>|101111>
i= 2, en [eV]=0.43629670382916214, +0.596|011111>|011111>-0.312|011111>|101111>-0.180|011111>|110111>-0.093|011111>|111011>+0.064|011111>|111101>+0.021|011111>|111110>-0.312|101111>|011111>-0.595|101111>|101111>-0.035|101111>|110111>-0.018|101111>|111011>+0.012|101111>|111101>+0.004|101111>|111110>-0.180|110111>|011111>-0.035|110111>|101111>-0.093|111011>|011111>-0.018|111011>|101111>+0.064|111101>|011111>+0.012|111101>|101111>+0.021|111110>|011111>+0.004|111110>|101111>
i= 3, en [eV]=1.2505228566304396, +0.687|011111>|011111>+0.688|101111>|101111>-0.145|110111>|110111>-0.082|110111>|111011>+0.074|110111>|111101>+0.022|110111>|111110>-0.082|111011>|110111>-0.051|111011>|111011>+0.046|111011>|111101>+0.014|111011>|111110>+0.074|111101>|110111>+0.046|111101>|111011>-0.040|111101>|111101>-0.012|111101>|111110>+0.022|111110>|110111>+0.014|111110>|111011>-0.012|111110>|111101>-0.007|111110>|111110>
i= 4, en [eV]=1.9393654779292357, +0.656|011111>|110111>+0.252|011111>|111011>-0.070|011111>|111101>-0.037|011111>|111110>+0.013|101111>|110111>+0.005|101111>|111011>-0.001|101111>|111101>-0.656|110111>|011111>-0.013|110111>|101111>-0.252|111011>|011111>-0.005|111011>|101111>+0.070|111101>|011111>+0.001|111101>|101111>+0.037|111110>|011111>
i= 5, en [eV]=1.940122558236312, +0.013|011111>|110111>+0.005|011111>|111011>-0.001|011111>|111101>-0.656|101111>|110111>-0.252|101111>|111011>+0.070|101111>|111101>+0.037|101111>|111110>-0.013|110111>|011111>+0.656|110111>|101111>-0.005|111011>|011111>+0.252|111011>|101111>+0.001|111101>|011111>-0.070|111101>|101111>-0.037|111110>|101111>
i= 6, en [eV]=2.9352449224310595, +0.176|011111>|011111>-0.061|011111>|101111>+0.653|011111>|110111>+0.146|011111>|111011>+0.127|011111>|111101>+0.004|011111>|111110>-0.061|101111>|011111>-0.176|101111>|101111>+0.027|101111>|110111>+0.006|101111>|111011>+0.005|101111>|111101>+0.653|110111>|011111>+0.027|110111>|101111>+0.146|111011>|011111>+0.006|111011>|101111>+0.127|111101>|011111>+0.005|111101>|101111>+0.004|111110>|011111>
i= 7, en [eV]=2.9363046696777113, -0.061|011111>|011111>-0.176|011111>|101111>-0.027|011111>|110111>-0.006|011111>|111011>-0.005|011111>|111101>-0.176|101111>|011111>+0.061|101111>|101111>+0.653|101111>|110111>+0.146|101111>|111011>+0.127|101111>|111101>+0.004|101111>|111110>-0.027|110111>|011111>+0.653|110111>|101111>-0.006|111011>|011111>+0.146|111011>|101111>-0.005|111101>|011111>+0.127|111101>|101111>+0.004|111110>|101111>
i= 8, en [eV]=4.660579248302518, -0.141|011111>|011111>-0.141|101111>|101111>-0.856|110111>|110111>-0.294|110111>|111011>+0.125|110111>|111101>+0.049|110111>|111110>-0.294|111011>|110111>-0.084|111011>|111011>+0.018|111011>|111101>+0.009|111011>|111110>+0.125|111101>|110111>+0.018|111101>|111011>+0.106|111101>|111101>+0.016|111101>|111110>+0.049|111110>|110111>+0.009|111110>|111011>+0.016|111110>|111101>+0.006|111110>|111110>
i= 9, en [eV]=5.079583342421291, +0.194|011111>|110111>-0.326|011111>|111011>+0.591|011111>|111101>+0.084|011111>|111110>+0.004|101111>|110111>-0.006|101111>|111011>+0.011|101111>|111101>+0.002|101111>|111110>-0.194|110111>|011111>-0.004|110111>|101111>+0.326|111011>|011111>+0.006|111011>|101111>-0.591|111101>|011111>-0.011|111101>|101111>-0.084|111110>|011111>-0.002|111110>|101111>
The ground state is a linear combination of products of two Slater determinants, namely \(|011111\rangle| 101111\rangle\), where the first one (\(|011111\rangle\)) is the Slater determinant in the alpha-channel, and the second one (\(|101111\rangle\)) is the determinant in the beta-channel.
Slater determinants are in the Fock-representation, where the orbitals with the lowest energy is on the right. For example, \(|011111\rangle\) means that the 5 orbitals with the lowest energy are occupied and the orbital with the highest energy is unoccupied.
Finally, the cutoff in visualize_correlated_state is the cutoff for Slater determinants. Only Slater determinants with prefactors larger (in absolute value) than the cutoff are printed.
We can visualize the contribution of the various Slater determinants to a given excited many-body state.
[10]:
import matplotlib.pyplot as plt
index1 = 0
index2 = 1
fig = plt.figure()
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.set_title('Ground state')
ax2.set_title('1st Excited state')
ax1.set_xlabel('Slater determinant index')
ax1.set_ylabel('Contribution of Slater determinant')
ax2.set_xlabel('Slater determinant index')
ax1.bar([i for i in range(solution['evcs'][index1].flatten().shape[0])],
abs(solution['evcs'][index1].flatten()))
ax2.bar([i for i in range(solution['evcs'][index2].flatten().shape[0])],
abs(solution['evcs'][index2].flatten()))
plt.show()
Solution of QDET for \(m_S=1\)¶
We demonstate how to diagonalize the effective Hamiltonian for \(m_S=1\). In this case the alpha spin cannel will have more electrons than the beta spin-channel.
We re-download the dataset obtained with WEST.
[11]:
%%bash
mkdir -p west.wfreq.save
wget -N -q https://west-code.org/doc/training/nv_diamond_63/wfreq.in
wget -N -q https://west-code.org/doc/training/nv_diamond_63/wfreq.json -O west.wfreq.save/wfreq.json
We now diagonalize the effective Hamiltonian choosing manually the number of electrons.
[12]:
from westpy.qdet import QDETResult
# construct object for effective Hamiltonian
effective_hamiltonian = QDETResult(filename='west.wfreq.save/wfreq.json')
nelec=(int(np.sum(effective_hamiltonian.occupation)//2)+1,int(np.sum(effective_hamiltonian.occupation)//2)-1)
print(nelec)
# diagonalize Hamiltonian
solution = effective_hamiltonian.solve(nelec=nelec)
(6, 4)
-----------------------------------------------------
===============================================================
Building effective Hamiltonian...
nspin: 1
occupations: [[2. 2. 2. 2. 1. 1.]]
===============================================================
| diag[1RDM - 1RDM(GS)] | ||||||||
|---|---|---|---|---|---|---|---|---|
| E [eV] | char | 87 | 122 | 123 | 126 | 127 | 128 | |
| 0 | 0.000 | 3- | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
| 1 | 1.939 | 3- | -0.003 | -0.010 | -0.127 | -0.860 | 1.000 | 0.000 |
| 2 | 1.940 | 3- | -0.003 | -0.010 | -0.127 | -0.860 | 0.000 | 1.000 |
| 3 | 5.080 | 3- | -0.014 | -0.698 | -0.213 | -0.075 | 1.000 | 0.000 |
| 4 | 5.080 | 3- | -0.014 | -0.698 | -0.213 | -0.075 | 0.000 | 1.000 |
| 5 | 6.163 | 3- | -0.001 | -0.278 | -0.657 | -0.064 | 0.999 | 0.001 |
| 6 | 6.163 | 3- | -0.001 | -0.278 | -0.657 | -0.064 | 0.001 | 0.999 |
| 7 | 7.030 | 3- | -0.017 | -0.691 | -0.353 | -0.939 | 1.000 | 1.000 |
| 8 | 8.016 | 3- | -0.003 | -0.303 | -0.761 | -0.933 | 1.000 | 1.000 |
| 9 | 11.138 | 3- | -0.019 | -0.976 | -0.877 | -0.127 | 1.000 | 1.000 |
-----------------------------------------------------
We see that we have found only triplet states. In order to find singlet and triplet states, it is necessary to have \(m_S=0\), which is the assumption made by default in the solve method.