renamed folder

code/entropy_new.py

+import numpy as np
+import sys
+import math
+from numpy.core.arrayprint import format_float_positional
+
+from numpy.core.fromnumeric import trace
+
+from joblib import delayed, Parallel
+import numpy as np
+import scipy.special
+import pickle
+
+import numba
+from numba import uint64, float32, float64
+import sys
+import os
+
+import lattice_symmetries
+from lattice_symmetries import *
+
+@numba.jit("uint64(uint64)", nogil=True, nopython=True)
+def _hamming_weight(n: int):
+    """The Hamming weight counts the non zero-values in a bit string. 
+    In the case of spin-chains this is the amount of up-spins.
+    Example: 10 -> 1010 -> hamming weight = 2
+    See https://stackoverflow.com/a/9830282
+    """
+    n = (n & 0x5555555555555555) + ((n & 0xAAAAAAAAAAAAAAAA) >> 1)
+    n = (n & 0x3333333333333333) + ((n & 0xCCCCCCCCCCCCCCCC) >> 2)
+    n = (n & 0x0F0F0F0F0F0F0F0F) + ((n & 0xF0F0F0F0F0F0F0F0) >> 4)
+    n = (n & 0x00FF00FF00FF00FF) + ((n & 0xFF00FF00FF00FF00) >> 8)
+    n = (n & 0x0000FFFF0000FFFF) + ((n & 0xFFFF0000FFFF0000) >> 16)
+    n = (n & 0x00000000FFFFFFFF) + ((n & 0xFFFFFFFF00000000) >> 32)
+    return n
+
+def sector_dm(trace_states, sub_states, amplitudes):
+    distinct_sub_states = list(set(sub_states))
+    distinct_trace_states = list(set(trace_states))
+    n = len(distinct_sub_states)
+    matrix = np.zeros((n,n))
+    for state in distinct_trace_states:
+        
+        single_trace_result = [(amplitudes[i], distinct_sub_states.index(sub_states[i])) 
+                                for i, x in enumerate(trace_states) if x == state]
+        
+        psi = np.zeros(n)
+        for x in single_trace_result:  
+            psi[x[1]] += x[0]
+        matrix += np.outer(psi,psi)
+    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):
+    """
+    :sub_dim:
+    :amplitudes:
+    :states_of_amplitudes:
+    """
+    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])
+    #distinct_states_out = list(set(states_out))
+    
+    rho = []
+    for sub_hamming in range(max(0, hamming - (dim - sub_dim)), min(hamming, sub_dim) + 1):
+        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]
+        sector_amplitudes = [amplitudes[i] for i in indices]
+        rho.append(sector_dm(sector_trace_states, sector_sub_states, sector_amplitudes))
+    return rho
+
+    
+def lambda_log_lambda(x):
+    """
+    Von Neumann Entropy
+    :x: array of eigenvalues of a sector density matrix
+    """
+    y = np.where(x > 1e-7, np.log(x), 0.0)
+    y *= x
+    return y
+
+def compute_entropy(rho):
+    """
+    :rho: list of sector density matrices
+    """
+    entropy = 0
+    for i in range(len(rho)):
+            spectrum, _ = np.linalg.eigh(rho[i])
+            entropy -= lambda_log_lambda(spectrum).sum()
+    return entropy
+
+
+def idxs_rearrangement(old_states, rearrangement):
+    """
+    Rearranges the old states according to the spin number list rearrangement
+    :old_states: index of old states
+    :rearrangement: list of new spin positions
+    :return: index of rearranged spins
+    """
+    new_states = []
+    batch_size = 1000000
+    number_spins = len(rearrangement)
+
+    # Rearranges the columns in the old spin matrix and computes the index of the rearranged spins
+    for i in range(len(old_states) // batch_size + 1):
+        print(i * 1. / (len(old_states) // batch_size + 1), "not zero if len(old_states) is large", flush=True)
+        new_states.append(
+            spin_to_index(
+                index_to_spin(
+                    old_states[i * batch_size:np.min([(i + 1) * batch_size, len(old_states)])], number_spins
+                    )[:, rearrangement].astype(bool), number_spins
+                    )
+                    )
+
+    return np.concatenate(new_states)
\ No newline at end of file

code/heisenberg.py

+import numpy as np
+import yaml
+import time
+import matplotlib.pyplot as plt
+
+from create_model import spin_model
+from entropy_new import reduced_dm, compute_entropy
+
+def test_xxx():
+    model_name = 'xxz'
+    number_spins = 4
+    periodic = True
+    param = -2.0
+    if param < -3.0:
+        hamming_weight = 0
+    else: 
+        hamming_weight=None
+    model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic, 
+                            param=param, hamming_weight=hamming_weight)
+    model.compute_ew_and_ev()
+    print('EIGENSTATE at h/J = ', param, 'is: ', model.eigenstates[:,0])
+    sub_dim = 2
+    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, hamming_weight, gs, basis_states)
+    print(rhos)
+    entropy = compute_entropy(rhos)
+    print(entropy)
+
+def test_area_law():
+    model_name = 'xxz'
+    number_spins = 12
+    periodic = True
+    hamming_weight = number_spins // 2
+    params = np.linspace(-4.0, 4.0, 41)
+    
+    entropies = {}
+    sub_dims = [i for i in range(2,number_spins-1)]
+    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)
+        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)
+            bipartite_entropies.append(compute_entropy(rhos))
+        entropies[str(param)] = bipartite_entropies
+        del model
+    
+    data = {
+        'Model': model_name,
+        'Number_spins': number_spins,
+        'Periodic': periodic,
+        'Hamming_weight': hamming_weight,
+        'Sub_dims': sub_dims,
+        'Delta_over_J': params,
+        'Entropies': entropies
+    }
+
+    timestr = time.strftime("%Y%m%d-%H%M%S")
+    filename = 'test_area_law_' + timestr
+    with open('output/' + filename + '.yaml', 'w') as outfile:
+        yaml.dump(data, outfile, default_flow_style=False)
+    print(filename)
+    plt.figure(figsize=(12,12))
+    i = 0
+    for key in entropies:
+        if i % 5 == 0:
+            plt.plot(np.array(sub_dims)/number_spins, entropies[key], label='Delta_over_J:' + key)
+        i += 1
+    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.grid(True)
+
+    plt.savefig('output/' + filename + '.jpg')
+
+
+    plt.figure(figsize=(20,20))
+    half_chain_ee = []
+    for key in entropies:
+        half_chain_ee.append(entropies[key][7])
+    plt.plot(params, half_chain_ee, label='Delta_over_J:' + key)
+    plt.legend()
+    plt.xlabel('Delta / J')
+    plt.ylabel('Bipartite Entanglement Entropy')
+    plt.title('Model:' + model_name + ' Number spins:' + str(number_spins))#, ' Periodic:' + str(periodic))
+    plt.grid(True)
+
+    plt.savefig('output/' + 'entropy' + timestr + '.jpg')

code/package_fails.py

+import numpy as np
+import yaml
+import time
+import matplotlib.pyplot as plt
+
+from create_model import spin_model
+from entropy_new import _hamming_weight, reduced_dm, compute_entropy
+
+def test_symmetries_xxx():
+    model_name = 'xxz'
+    number_spins = [8,10,12]
+    periodic = True
+    param = 1.0
+    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')
+        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')
+        model = spin_model(model_name=model_name, number_spins=length, periodic=periodic, 
+                            param=param, hamming_weight=hamming_weight, use_symmetries=False)
+        model.compute_ew_and_ev()
+        #print('EIGENSTATE at Delta/J = ', param, 'is: ', model.eigenstates[:,0])
+        del model
\ No newline at end of file

code/tfim.py

+import numpy as np
+import yaml
+import time
+import matplotlib.pyplot as plt
+
+from create_model import spin_model
+from entropy_new import reduced_dm, compute_entropy
+
+def test_tfim():
+    model_name = 'tfim'
+    number_spins = 4
+    periodic = True
+    hamming_weight = number_spins // 2
+    param = 0.0
+    model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic, 
+                            param=param, hamming_weight=None)
+    model.compute_ew_and_ev()
+    print('EIGENSTATE at h/J = ', param, 'is: ', model.eigenstates[:,0])
+    sub_dim = 2
+    model.compute_ew_and_ev()
+    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, hamming_weight, gs, basis_states)
+    print(rhos)
+    entropy = compute_entropy(rhos)
+    print(entropy)
+
+def test_area_law_tfim():
+    model_name = 'tfim'
+    number_spins = 12
+    periodic = True
+    params = np.linspace(0.0, 4.0, 21)
+    
+    entropies = {}
+    sub_dims = [i for i in range(1,number_spins)]
+    for param in params:
+        bipartite_entropies = []
+        model = spin_model(model_name=model_name, number_spins=number_spins, periodic=periodic, 
+                                param=param, hamming_weight=None)
+        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)
+            entropy = 0
+            for hamming_weight in range(1):
+                rhos = reduced_dm(sub_dim, number_spins, 6, gs, basis_states)
+                entropy += compute_entropy(rhos)
+            bipartite_entropies.append(entropy)
+        entropies[str(param)] = bipartite_entropies
+        del model
+    
+    data = {
+        'Model': model_name,
+        'Number_spins': number_spins,
+        'Periodic': periodic,
+        'Hamming_weight': hamming_weight,
+        'Sub_dims': sub_dims,
+        'Delta_over_J': params,
+        'Entropies': entropies
+    }
+
+    timestr = time.strftime("%Y%m%d-%H%M%S")
+    filename = 'test_area_law_tfim_' + timestr
+    with open('output/' + filename + '.yaml', 'w') as outfile:
+        yaml.dump(data, outfile, default_flow_style=False)
+    print(filename)
+    plt.figure(figsize=(12,12))
+    i = 0
+    for key in entropies:
+        if i % 5 == 0:
+            plt.plot(np.array(sub_dims)/number_spins, entropies[key], label='Delta_over_J:' + key)
+        i += 1
+    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.grid(True)
+
+    plt.savefig('output/' + filename + '.jpg')
+
+
+    plt.figure(figsize=(20,20))
+    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.legend()
+    plt.xlabel('Delta / J')
+    plt.ylabel('Bipartite Entanglement Entropy')
+    plt.title('Model:' + model_name + ' Number spins:' + str(number_spins))#, ' Periodic:' + str(periodic))
+    plt.grid(True)
+
+    plt.savefig('output/' + 'entropy_tfim_' + timestr + '.jpg')