Commit 799e42bf authored by mgassner's avatar mgassner
Browse files

some changes

parent 629cc12a
......@@ -34,7 +34,7 @@ def generate_symmetries(model_name, number_spins, periodic):
sites = np.arange(number_spins)
# Momentum in x direction with eigenvalue π
if periodic:
T = ((sites + 1) % number_spins).tolist()
T = (sites + 1) % number_spins
symmetries.append(ls.Symmetry(T, sector=number_spins // 2))
# Parity with eigenvalue π
P = sites[::-1]
......@@ -49,20 +49,20 @@ def create_edges(number_spins, periodic):
class spin_model:
#constructor
def __init__(self, model_name, number_spins, periodic, param, hamming_weight, use_symmetries=False):
def __init__(self, model_name, number_spins, periodic, param, hamming_weight, use_symmetries=False, spin_inversion=None):
self.__number_spins = number_spins
self.__periodic = periodic
self.symmetries = []
self.__spin_inversion = None
self.__symmetries = []
self.__spin_inversion = spin_inversion
if use_symmetries:
self.symmetries = generate_symmetries(model_name, number_spins, periodic)
self.__spin_inversion = -1
self.__symmetries = generate_symmetries(model_name, number_spins, periodic)
self.__name = model_name
self.__matrix = create_2spin_matrix(model_name, param)
self.edges = create_edges(number_spins, periodic)
self.__hamming_weight = hamming_weight
self.basis = None
self.create_basis()
#creating Hamiltonian
self.__matrix = create_2spin_matrix(model_name, param)
self.edges = create_edges(number_spins, periodic)
self.hamiltonian = None
self.create_hamiltonian()
self.eigenvalues = None
......@@ -76,7 +76,7 @@ class spin_model:
def create_basis(self):
# Constructing the basis
symmetry_group = ls.Group(self.symmetries)
symmetry_group = ls.Group(self.__symmetries)
self.basis = ls.SpinBasis(symmetry_group, number_spins=self.__number_spins,
hamming_weight=self.__hamming_weight, spin_inversion=self.__spin_inversion)
self.basis.build()
......
......@@ -33,10 +33,11 @@ def _hamming_weight(n: int):
n = (n & 0x00000000FFFFFFFF) + ((n & 0xFFFFFFFF00000000) >> 32)
return n
def sector_dm(trace_states, sub_states, amplitudes):
def sector_dm(trace_states, sub_states, amplitudes, spin_inversion):
distinct_sub_states = list(set(sub_states))
distinct_trace_states = list(set(trace_states))
n = len(distinct_sub_states)
print(trace_states, sub_states)
matrix = np.zeros((n,n))
for state in distinct_trace_states:
......@@ -46,17 +47,14 @@ def sector_dm(trace_states, sub_states, amplitudes):
psi = np.zeros(n)
for x in single_trace_result:
psi[x[1]] += x[0]
matrix += np.outer(psi,psi)
matrix_add = np.outer(psi,psi)
if spin_inversion is not None:
matrix += 0.5 * (matrix_add + matrix_add[::-1,::-1])
else:
matrix += matrix_add
return matrix
"""for sub_state in distinct_states_sub:
for out_state in distinct_states_out:
amplitudes_sub = [(amplitudes[i], states_sub[i]) for i, x in enumerate(states_out) if x == state]
amplitudes_sub = [sum(x[0] for x in amplitudes_sub if x[1]==i) for i in distinct_states_sub]
rho.append(np.outer(amplitudes_sub, amplitudes_sub))
print(rho)"""
#Goal:
def reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes):
def reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes, spin_inversion=False):
"""
:sub_dim: amount of spins in the subsystem (uint)
:dim: amount of spins of total system (uint)
......@@ -65,21 +63,24 @@ def reduced_dm(sub_dim, dim, hamming, amplitudes, states_of_amplitudes):
:states_of_amplitudes: state which belongs to amplitude (1d-array of uints)
:return: List of sector density matrices
"""
if spin_inversion is not None:
dim -= 1
sub_states = states_of_amplitudes >> (dim - sub_dim)
trace_states = states_of_amplitudes & (2**(dim-sub_dim)-1)
hamming_of_sub_states = np.array([_hamming_weight(n) for n in sub_states])
#print(sub_states, trace_states, hamming_of_sub_states, amplitudes)
rho = []
if hamming is None:
sub_hammings = np.arange(0, sub_dim+1)
rho.append(sector_dm(trace_states, sub_states, amplitudes, spin_inversion))
else:
sub_hammings = np.arange(max(0, hamming - (dim - sub_dim)), min(hamming, sub_dim) + 1)
for sub_hamming in sub_hammings:
for sub_hamming in range(max(0, hamming - (dim - sub_dim)), min(hamming, sub_dim) + 1):
print('Sub hamming is', sub_hamming)
indices = np.where(hamming_of_sub_states == sub_hamming)[0].tolist()
sector_sub_states = [sub_states[i] for i in indices]
sector_trace_states = [trace_states[i] for i in indices]
#print(sector_trace_states)
sector_amplitudes = [amplitudes[i] for i in indices]
rho.append(sector_dm(sector_trace_states, sector_sub_states, sector_amplitudes))
rho.append(sector_dm(sector_trace_states, sector_sub_states, sector_amplitudes, spin_inversion))
return rho
......@@ -89,7 +90,7 @@ def lambda_log_lambda(x):
:x: array of eigenvalues of a sector density matrix
:return: negative Bipartite Entanglement entropy of a sector density_matrix
"""
y = np.where(x > 1e-7, np.log(x) / np.log(2), 0.0)
y = np.where(x > 1e-7, np.log(x) / np.log(2), 1e-17)
y *= x
return y
......
from create_model import generate_symmetries
import numpy as np
import yaml
import time
......@@ -6,50 +7,59 @@ import matplotlib.pyplot as plt
from create_model import spin_model
from entropy_new import reduced_dm, compute_entropy
def test_xxx():
def test_xxz():
model_name = 'xxz'
number_spins = 4
number_spins = 10
periodic = True
param = -2.0
if param < -3.0:
hamming_weight = 0
spin_inversion = -1
param = 1.0
if param < -1.0:
hamming_weight = None
else:
hamming_weight=None
hamming_weight = number_spins // 2
model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic,
param=param, hamming_weight=hamming_weight)
param=param, hamming_weight=hamming_weight, use_symmetries=True,spin_inversion=spin_inversion)
model.compute_ew_and_ev()
print('EIGENSTATE at h/J = ', param, 'is: ', model.eigenstates[:,0])
sub_dim = 2
sub_dims = np.arange(1, number_spins)
number_spins = model.number_spins()
basis_states = model.basis.states
gs = model.eigenstates[:,0]
print('Number Spins: ', model.basis.number_spins)
print('States', model.basis.states)
rhos = reduced_dm(sub_dim, number_spins, 2, gs, basis_states)
print(rhos)
"""for sub_dim in sub_dims:
print('Sub Dimension is ', sub_dim)
rhos = reduced_dm(sub_dim, number_spins, hamming_weight, gs, basis_states, spin_inversion)
#print(rhos)
entropy = compute_entropy(rhos)
print(entropy)
print('-------------------------------------------')
"""
def test_area_law():
model_name = 'xxz'
number_spins = 12
number_spins = 14
periodic = True
hamming_weight = number_spins // 2
params = np.linspace(-4.0, 4.0, 41)
params = np.linspace(-2.0, 2.0, 41)
spin_inversion=True
entropies = {}
sub_dims = [i for i in range(2,number_spins-1)]
sub_dims = [i for i in range(1,number_spins)]
for param in params:
if param < -1.0:
hamming = None
else:
hamming = hamming_weight
bipartite_entropies = []
model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic,
param=param, hamming_weight=hamming_weight)
param=param, hamming_weight=hamming, spin_inversion=spin_inversion)
model.compute_ew_and_ev()
gs = model.eigenstates[:,0]
basis_states = model.basis.states
for sub_dim in sub_dims:
print('----------------------------------')
print('Sub Dimension is ', sub_dim)
rhos = reduced_dm(sub_dim, number_spins, hamming_weight, gs, basis_states)
rhos = reduced_dm(sub_dim, number_spins, hamming_weight, gs, basis_states, spin_inversion)
bipartite_entropies.append(compute_entropy(rhos))
entropies[str(param)] = bipartite_entropies
del model
......@@ -78,7 +88,7 @@ def test_area_law():
plt.legend()
plt.xlabel('x/L')
plt.ylabel('Bipartite Entanglement Entropy')
plt.title('Model:' + model_name + ' Number spins:' + str(number_spins), ' Periodic:' + str(periodic))
plt.title('Model:' + model_name + ' Number spins:' + str(number_spins))#, ' Periodic:' + str(periodic))
plt.grid(True)
plt.savefig('output/' + filename + '.jpg')
......
import numpy as np
import scipy.sparse.linalg
import lattice_symmetries as ls
def github_test(number_spins=10, spin_inversion=-1, use_symmetries=True):
number_spins = number_spins # System size
hamming_weight = number_spins // 2 # Hamming weight (i.e. number of spin ups)
# Constructing symmetries
symmetries = []
if use_symmetries:
sites = np.arange(number_spins)
# Momentum in x direction with eigenvalue π
T = (sites + 1) % number_spins
symmetries.append(ls.Symmetry(T, sector=number_spins // 2))
# Parity with eigenvalue π
P = sites[::-1]
symmetries.append(ls.Symmetry(P, sector=1))
# Constructing the group
symmetry_group = ls.Group(symmetries)
#print("Symmetry group contains {} elements".format(len(symmetry_group)))
# Constructing the basis
basis = ls.SpinBasis(
symmetry_group, number_spins=number_spins, hamming_weight=hamming_weight, spin_inversion=spin_inversion
)
basis.build() # Build the list of representatives, we need it since we're doing ED
print("Hilbert space dimension is {}".format(basis.number_states))
# Heisenberg Hamiltonian
# fmt: off
σ_x = np.array([ [0, 1]
, [1, 0] ])
σ_y = np.array([ [0 , -1j]
, [1j, 0] ])
σ_z = np.array([ [1, 0]
, [0, -1] ])
# fmt: on
σ_p = σ_x + 1j * σ_y
σ_m = σ_x - 1j * σ_y
matrix = 0.5 * (np.kron(σ_p, σ_m) + np.kron(σ_m, σ_p)) + np.kron(σ_z, σ_z)
edges = [(i, (i + 1) % number_spins) for i in range(number_spins)]
hamiltonian = ls.Operator(basis, [ls.Interaction(matrix, edges)])
# Diagonalize the Hamiltonian using ARPACK
eigenvalues, eigenstates = ls.diagonalize(hamiltonian, k=1)
print("Ground state energy is {:.10f}".format(eigenvalues[0]))
check = {'4': None, '6': None, '8': None, '10': -18.06178542}
#assert np.isclose(eigenvalues[0], check[str(number_spins)])
\ No newline at end of file
from heisenberg import test_xxx, test_area_law
from heisenberg import test_xxz, test_area_law
from tfim import test_tfim, test_area_law_tfim
from package_fails import test_symmetries_xxx
from ls_example_from_github import github_test
if __name__ == "__main__":
#test_xxx()
#test_xxz()
#test_area_law()
test_area_law_tfim()
#test_area_law_tfim()
#test_tfim()
#test_symmetries_xxx()
\ No newline at end of file
test_symmetries_xxx()
"""print('------------------------------------------------')
print('---With symmetries and spin inversion---')
github_test()
print('---Without symmetries, but spin inversion---')
github_test(use_symmetries=False)
print('---With symmetries, but no spin inversion---')
github_test(spin_inversion=None)
print('---Without symmetries and spin inversion---')
github_test(spin_inversion=None, use_symmetries=False)
print('------------------------------------------------')
print('---With symmetries and spin inversion---')
github_test()
print('---Without symmetries, but spin inversion---')
github_test(use_symmetries=False)
print('---With symmetries, but no spin inversion---')
github_test(spin_inversion=None)
print('---Without symmetries and spin inversion---')
github_test(spin_inversion=None, use_symmetries=False)"""
\ No newline at end of file
......@@ -14,16 +14,28 @@ def test_symmetries_xxx():
print('Test if gs-energy stays the same if we use symmetries for the ls package:')
for length in number_spins:
hamming_weight = length // 2
print('------- Number spins is' + str(length) + '-----------' )
print('With symmetries')
print('\n' , '------- Number spins is: ' + str(length) + '-----------' )
print('---With symmetries and spin inversion---')
model = spin_model(model_name=model_name, number_spins=length, periodic=periodic,
param=param, hamming_weight=hamming_weight, use_symmetries=True, spin_inversion=-1)
model.compute_ew_and_ev()
#print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
del model
print('---Without symmetries, but spin inversion---')
model = spin_model(model_name=model_name, number_spins=length, periodic=periodic,
param=param, hamming_weight=hamming_weight, spin_inversion=-1)
model.compute_ew_and_ev()
#print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
del model
print('---With symmetries, but without spin inversion---')
model = spin_model(model_name=model_name, number_spins=length, periodic=periodic,
param=param, hamming_weight=hamming_weight, use_symmetries=True)
model.compute_ew_and_ev()
#print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
del model
print('Without symmetries')
print('---Without symmetries and spin inversion---')
model = spin_model(model_name=model_name, number_spins=length, periodic=periodic,
param=param, hamming_weight=hamming_weight, use_symmetries=False)
param=param, hamming_weight=hamming_weight)
model.compute_ew_and_ev()
#print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
del model
\ No newline at end of file
......@@ -8,29 +8,33 @@ from entropy_new import reduced_dm, compute_entropy
def test_tfim():
model_name = 'tfim'
number_spins = 4
number_spins = 6
periodic = True
hamming_weight = number_spins // 2
hamming_weight = None
param = 1.0
spin_inversion= None
model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic,
param=param, hamming_weight=None)#hamming_weight)
param=param, hamming_weight=None, spin_inversion=spin_inversion)
model.compute_ew_and_ev()
print('EIGENSTATE at h/J = ', param, 'is: ', model.eigenstates[:,0])
sub_dim = 2
model.compute_ew_and_ev()
sub_dims = np.arange(1, number_spins)
bipartite_entropies = []
gs = model.eigenstates[:,0]
basis_states = model.basis.states
print('Number Spins: ', model.basis.number_spins)
print('States', model.basis.states)
rhos = reduced_dm(sub_dim, number_spins, 0, gs, basis_states)
print(rhos)
entropy = compute_entropy(rhos)
print(entropy)
for sub_dim in sub_dims:
print('----------------------------------')
print('Sub Dimension is ', sub_dim)
rhos = reduced_dm(sub_dim, number_spins, hamming_weight, gs, basis_states, spin_inversion)
bipartite_entropies.append(compute_entropy(rhos))
del model
print(sub_dims, bipartite_entropies)
def test_area_law_tfim():
model_name = 'tfim'
number_spins = 14
periodic = True
spin_inversion = True
params = np.linspace(0.0, 1.5, 21)
hamming_weight = None # TODO: ?????
entropies = {}
......@@ -38,7 +42,7 @@ def test_area_law_tfim():
for param in params:
bipartite_entropies = []
model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic,
param=param, hamming_weight=hamming_weight)
param=param, hamming_weight=hamming_weight, spin_inversion=spin_inversion)
model.compute_ew_and_ev()
gs = model.eigenstates[:,0]
basis_states = model.basis.states
......@@ -67,9 +71,12 @@ def test_area_law_tfim():
print(filename)
plt.figure(figsize=(12,12))
i = 0
x = np.array(sub_dims)/number_spins
if spin_inversion:
x -= 1/number_spins
for key in entropies:
if i % 5 == 0:
plt.plot(np.array(sub_dims)/number_spins, entropies[key], label='Delta_over_J:' + key)
plt.plot(np.array(sub_dims)/number_spins, entropies[key], label='h over J:' + key)
i += 1
plt.legend()
plt.xlabel('x/L')
......@@ -84,7 +91,7 @@ def test_area_law_tfim():
half_chain_ee = []
for key in entropies:
half_chain_ee.append(entropies[key][5])
plt.plot(params, half_chain_ee, label='Delta_over_J:' + key)
plt.plot(params, half_chain_ee, label='h over J:' + key)
plt.legend()
plt.xlabel('Delta / J')
plt.ylabel('Bipartite Entanglement Entropy')
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment