"""
This module provides various types of initial states. The initial states for
the qIMPS routine can be defined via
* :py:func:`State.GenerateqInitialState` using the class
:py:class:`qProductState`.
Non-symmetry conserving MPS states may be obtained through
* :py:class:`State.ProductState`
* :py:class:`State.PartlyProduct`
"""
import math as ma
from math import sqrt, floor, ceil
import numpy as np
#import numpy.linalg as la
from copy import deepcopy
from sys import version_info
if(version_info[0] < 3):
from itertools import izip as myzip
else:
myzip = zip
[docs]def nint(x):
"""
Return the integer nearest x (only defined for :math:`x \ge 0`).
"""
y = round(x) - .5
return int(y) + (y > 0)
[docs]def tupleadd(xs, ys, nqns, ncs):
"""
Addition of tuples with intermediate conversion to numpy array.
"""
tmp = np.array(xs, dtype=np.int16) + np.array(ys, dtype=np.int16)
for ii in range(nqns, ncs): tmp[ii] = tmp[ii]%2
return tuple(tmp)
[docs]def tupledist(xl, xr, mf, nqns):
"""
Distance of tuples calculated with conversion to numpy array.
"""
mf_ = np.array(mf)
return np.sqrt(np.sum((np.array(xl)[:nqns] - mf_[:nqns])**2) +
np.sum((np.array(xr)[:nqns] - mf_[:nqns])**2))
[docs]def GenerateqInitialState(fileID, qnums, pnums, L, qOperators):
"""
Generates an initial state obeying the symmetry according to
:py:class:`State.qProductState`.
**Arguments**
fileID : str
Basis for the complete file. The complete file name is
fileID + ``InitialState.dat``
qnums : list
List containing the quantum numbers of the U(1) symmetry.
pnums : list
List containint the quantum numbers for the Z2 symmetry.
L : int
number of sites in the system
qOperators : instance of :py:class:`Ops.OperatorList`
The symmetry conserving operators have to be generated in
:py:class:`Ops.OperatorList` before creating a qProductState.
"""
qstate = qProductState(qnums, pnums, L, qOperators)
return qstate.write(fileID)
[docs]def GenerateParityState(fileID,L,N0,N1):
"""
Generate a product state on L sites with N0 particles in state
:math:`|0 \\rangle` and N1 particles in state :math:`|1 \\rangle`.
N0 + N1 cannot be greater than the system size, but there can
be empty site :math:`| E \\rangle`. The particles are placed in
the middle:
.. math::
| E \ldots E 0 \ldots 0 1 \ldots 1 E \ldots E \\rangle
**Arguments**
fileID : str
Basis for the complete file. The complete file name is
fileID + ``InitialState.dat``
L : int
system size.
N0 : int
Number of particles in the state 0.
N1 : int
Number of particles in the state 1.
"""
N=N0+N1
#Place them all in the middle
empty=L-N
le=int(floor(0.5*empty))
re=int(ceil(0.5*empty))
qInitialFileName=fileID+'InitialState.dat'
qFile=open(qInitialFileName,'w')
site=1
print('le re N0 N1 sum',le, re,N0, N1,le+re+N0+N1)
for i in range(le):
print('empty site',site)
qFile.write('%16i'%(0)+'%16i'%(1)+'\n')#empty state, dimension 1
qFile.write('%30.15E'%(1.0)+'\n')#empty state, dimension 1
site+=1
for i in range(N0):
print('N_0 site',site)
qFile.write('%16i'%(1)+'%16i'%(2)+'\n')#N=1 state, dimension 2
qFile.write('%30.15E'%(1.0)+'%30.15E'%(0.0)+'\n')#Block format [c_0, c_1]
site+=1
for i in range(N1):
print('N_1 site',site)
qFile.write('%16i'%(1)+'%16i'%(2)+'\n')#N=1 state, dimension 2
qFile.write('%30.15E'%(0.0)+'%30.15E'%(1.0)+'\n')#Block format [c_0, c_1]
site+=1
for i in range(re):
print('empty site',site)
qFile.write('%16i'%(0)+'%16i'%(1)+'\n')#empty state, dimension 1
qFile.write('%30.15E'%(1.0)+'\n')#empty state, dimension 1
site+=1
qFile.close()
return qInitialFileName
[docs]class qProductState:
"""
Generate an initial state on L sites with the quantum numbers qnums and
the Hilbert space in qOperators. The routine attempts to make the state
as uniform as possible by defining a mean filling and filling as many
sites as possible with that filling. The "exceptional" sites are then
filled with the excess quantum numbers.
This routine is a work in progress for multiple quantum numbers-MLW 090111.
**Arguments**
qnums : list
List containing the quantum numbers of the U(1) symmetry.
pnums : list
List containint the quantum numbers for the Z2 symmetry.
L : int
number of sites in the system
qOperators : instance of :py:class:`Ops.OperatorList`
The symmetry conserving operators have to be generated in
:py:class:`Ops.OperatorList` before creating a qProductState.
**Details**
1) In the first step generate four lists containing:
* meanfilling: rounded quantum number in system divided by
system size L. This is done for every quantum number.
* excess: Difference between the mean filling just calculated and
the targeted quantum number. Negative when mean filling too big.
Never used again.
* nexceptionall (!!! two l's at the end !!!): number of exceptional
sites.
* nrmex: sign/signum function of excess.
2) Specify particles injected on each iteration. ninject contains the
two times the mean filling (except single site) and substracts/adds
exceptionall on the last injections
3) Injections are carried out on two-site blocks. Therefore the quantum
numbers on a two-site block is derived and stored in a dictionary.
4) Select from those two-site blocks, for each injection, the injection
which is closest to the mean-filling. E.g. for one quantum number
with U(1) symmetry, choose (1, 1) over (2, 0) or (0, 2) for mean-filling
of two particles. Z2 symmetries have always a distance of 0.
5) Write out states is covered through
:py:func:`State.qProductState.write`.
"""
def __init__(self, qnums, pnums, L, qOperators):
"""
"""
# number of iterations of two-site iMPS
niter=int(ceil(0.5*L))
# Switch for single site at the end (occurs for odd #s of sites)
singlesite=int(niter-floor(0.5*L))
#print('total quantum numbers',qnums, L, niter, singlesite)
# Number of symetries and total number of symmetrie
nqns = len(qnums)
npns = len(pnums)
ncs = nqns + npns
allnums = deepcopy(qnums) + deepcopy(pnums)
# Find the meanfilling per site, the number of "exceptional" sites which
# do not have this filling, and the (signed) number of excess particles
# which must go into the exceptional sites
meanfilling = np.zeros((ncs), dtype=np.int16)
nexceptionall = np.zeros((ncs), dtype=np.int16)
nrmex = np.zeros((ncs), dtype=np.int16)
# Iterate through abelian and discrete symmetries
for ii in range(ncs):
meanfilling[ii] = nint(float(allnums[ii]) / L)
ex = allnums[ii] - meanfilling[ii] * L
nrmex[ii] = np.sign(ex)
nexceptionall[ii] = int(ceil(abs(ex)))
nexceptional = np.max(nexceptionall)
#print('nex niter-nex',nexceptional,niter-nexceptional)
#print('excess',excess,nrmex)
# Calculate the numbers of particles injected into each iteration of iMPS.
# Injection is mean-field plus correction if applicable
ninject = np.zeros((niter, ncs), dtype=np.int16)
for ii in range(ncs):
ninject[:(niter - nexceptionall[ii]), ii] = 2 * meanfilling[ii]
ninject[(niter - nexceptionall[ii]):, ii] = 2 * meanfilling[ii] \
+ nrmex[ii]
# Correct last one only if single-sites
if(singlesite):
for ii in range(ncs):
ninject[-1, ii] = meanfilling[ii] + nrmex[ii]
# find all possible total quantum numbers of a two-site block
TwoSiteKeys={}
count=0
for i in qOperators['q_d']:
for j in qOperators['q_d']:
comb = tupleadd(i, j, nqns, ncs)
if comb in TwoSiteKeys:
TwoSiteKeys[comb].append([i, j])
else:
TwoSiteKeys[comb] = [[i, j]]
# Find the set of single-site tuples such that the total quantum numbers
# and equal to the injected quantum numbers and the distance for U(1)
# |q_i-mean|^2+|q_j-mean|^2 is minimized. Note that there may be many such
# tuples, and we include them all to avoid breaking e.g. reflection
# symmetries. Then check that each tuple set represents a valid set of
# quantum numbers. (Z2 symmetry does not contribute to distance).
# Each element will contain info for one injection
blocks = []
tups = []
for ii in range(niter - singlesite):
try:
filling = tuple(ninject[ii, :])
halffilling = tuple(0.5 * ninject[ii, :])
tup = []
maxdist = np.inf
for tuples in TwoSiteKeys[filling]:
totdist = tupledist(tuples[0], tuples[1], halffilling, nqns)
if(totdist == maxdist):
# Equal distance, append to list
tup.append(tuples)
elif(totdist < maxdist):
# Smaller distance, reset
maxdist = totdist
tup = [tuples]
tups.append(tup)
blocks.append(len(tup))
except:
# Should be made more robust in the case where multiple quantum
# numbers exist
print("possible two site keys", TwoSiteKeys)
print("needed two site key", ninject[ii, :])
raise Exception("Cannot build a state with the given quantum " +
"numbers!")
# fill the remaining site, and be sure that the tuple for this site
# is legitimate
if singlesite:
filling = tuple(ninject[niter - 1, :])
if(filling in qOperators['q_d']):
blocks.append(1)
tups.append([[filling]])
else:
raise Exception("Cannot build a state with the given "
+ "quantum numbers!")
self.L = L
self.tups = tups
self.singlesite = singlesite
self.blocks = blocks
self.niter = niter
self.qOps = qOperators
[docs] def write(self, fileID):
"""
Write out the state.
**Arguments**
fileID : str
Basis for the complete file. The complete file name is
fileID + ``InitialState.dat``
"""
qInitialFileName = fileID + 'InitialState.dat'
qFile = open(qInitialFileName, 'w')
for i in range(self.niter - self.singlesite):
# number of blocks
qFile.write('%16i\n'%(self.blocks[i]))
for j in range(self.blocks[i]):
itupstr = ''
for q in self.tups[i][j][0]:
itupstr += '%16i'%(q)
jtupstr = ''
for q in self.tups[i][j][1]:
jtupstr += '%16i'%(q)
# tuples of the blocks
qFile.write(itupstr + jtupstr + '\n')
# write deg dimensions
qFile.write('%16i'%(self.qOps['q_d'][self.tups[i][j][0]]) +
'%16i'%(self.qOps['q_d'][self.tups[i][j][1]]) + '\n')
if(self.singlesite):
qFile.write('%16i\n'%(1))
itupstr = ''
for q in self.tups[self.niter - 1][0][0]:
itupstr += '%16i'%(q)
# tuples of the blocks
qFile.write(itupstr + '\n')
# write deg dimensions
tmp = self.tups[self.niter - 1][0][0]
qFile.write('%16i'%(self.qOps['q_d'][tmp])+'\n')
qFile.close()
return qInitialFileName
[docs]class QProductState():
"""
Generate an initial state defined by the user. The state is defined
in terms of the quantum numbers and must be a product state.
**Arguments**
qOperators : instance of :py:class:`Ops.OperatorList`
The symmetry conserving operators have to be generated in
:py:class:`Ops.OperatorList` before creating a qProductState.
qnums_per_site : 2d numpy array (floats), optional
The sites run over the columns, the different quantum numbers
run over the rows. At system with 10 sites and 2 abelian
quantum numbers has a (2 x 10) shape.
Default to None.
pnums_per_site : 2d numpy array (floats), optional
The sites run over the columns, the different quantum numbers
run over the rows. At system with 10 sites and 2 discrete
quantum numbers has a (2 x 10) shape.
Default to None.
func_subspace : callable, optional
Can be used to pass state for degenerate subspace. If not
given, an equal superposition is assumed. The function
must take the following arguments in exactly this order:
(1) dimension of the subspace (2) site index (3) the
q-quantum numbers for U1-type symmetries, and (4) the
p-quantum numbers for Z2-type symmetries. It must return
a vector of the dimension of the subspace.
**Details**
For degenerate subspaces in one symmetry sector, equal weight is
assigned to each state. The functions takes care of the shifting
and scaling the quantum numbers to integer-quantum numbers starting
at zero.
The system size is deducted from the second dimension of the array
qnums_per_site and pnums_per_site. Either qnums_per_site or
pnums_per_site must be given.
"""
def __init__(self, qOperators, qnums_per_site=None, pnums_per_site=None,
func_subspace=None):
# Get total number of symmetries, copy quantum numbers for shift & scale
if((qnums_per_site is None) and (pnums_per_site is None)):
raise Exception('No quantum numbers given.')
elif(qnums_per_site is None):
ll = pnums_per_site.shape[1]
psite = deepcopy(pnums_per_site)
qsite = np.zeros((0, ll))
elif(pnums_per_site is None):
ll = qnums_per_site.shape[1]
qsite = deepcopy(qnums_per_site)
psite = np.zeros((0, ll))
else:
ll = qnums_per_site.shape[1]
qsite = deepcopy(qnums_per_site)
psite = deepcopy(pnums_per_site)
totnums = qsite.shape[0] + psite.shape[0]
# Shift and scale
kk = 0
for ii in range(qsite.shape[0]):
qsite[ii, :] -= qOperators['q_shifts'][kk]
qsite[ii, :] /= qOperators['q_scales'][kk]
kk += 1
for ii in range(psite.shape[0]):
psite[ii, :] -= qOperators['q_shifts'][kk]
psite[ii, :] /= qOperators['q_scales'][kk]
kk += 1
# Now the quantum numbers can be converted to integers
qsite = np.array(qsite, dtype=np.int32)
psite = np.array(psite, dtype=np.int32)
for ii in range(ll):
filling = tuple(list(qsite[:, ii])
+ list(psite[:, ii]))
if(filling not in qOperators['q_d']):
raise Exception("The given quantum numbers are not available "
+ "in the choosen basis.")
self.L = ll
self.qOps = qOperators
self.qnums = qsite.shape[0]
self.pnums = psite.shape[0]
self.totnums = totnums
self.qsite = qsite
self.psite = psite
self.func = func_subspace
[docs] def write(self, filename):
"""
Write out the state.
**Arguments**
filename : str
The complete filename.
"""
fh = open(filename, 'w')
# Real state, q-numbers
fh.write(' T\n')
fh.write(' T\n')
# Number of sites / orthogonality center
fh.write('%12d\n'%(self.L))
fh.write('%12d\n'%(1))
# Write the canonization and the haslambda together
fh.write(' c F\n')
for ii in range(1, self.L):
fh.write(' r F\n')
fh.write(' F\n')
# Get the quantum numbers
Qleft = np.zeros((self.totnums), dtype=np.int32)
Qsite = np.zeros((self.totnums), dtype=np.int32)
Qrigh = np.zeros((self.totnums), dtype=np.int32)
# Write the tensors
for ii in range(self.L):
# Number of blocks (always one; product state)
fh.write('%12d %12d %12d\n'%(1, self.qnums, self.pnums))
# Write the hash and the size for the array with the qnumbers
fh.write('%15.12f%12d\n'%(0.0, 3 * self.totnums))
# Write the quantum numbers
Qsite[:self.qnums] = self.qsite[:, ii]
Qsite[self.qnums:] = self.psite[:, ii]
Qrigh[:self.qnums] = Qleft[:self.qnums] + self.qsite[:, ii]
Qrigh[self.qnums:] = (Qleft[self.qnums:] + self.psite[:, ii])%2
qstr = ''
for qq in Qleft:
qstr += '%12d'%(qq)
for qq in Qsite:
qstr += '%12d'%(qq)
for qq in Qrigh:
qstr += '%12d'%(qq)
fh.write(qstr + '\n')
# Write the tensor - rank
fh.write('%12d\n'%(3))
# Tensor: dimensions, entry by entry
locdim = self.qOps['q_d'][tuple(Qsite)]
fh.write('%12d%12d%12d'%(1, locdim, 1))
fh.write('%12d%12d%12d'%(1, locdim, 1))
fh.write('%12d%12d%12d\n'%(1, 2, 3))
tstr = ''
if(locdim == 1):
# No degeneracy, weight must be one.
tstr += '%15.12f\n'%(1.0)
elif(self.func is None):
# No function given, choose equal superposition
weight = 1.0 / np.sqrt(locdim)
for jj in range(locdim):
tstr += '%15.12f\n'%(weight)
else:
# Weight is given by user-specified function
vec = self.func(locdim, ii, self.qsite[:, ii],
self.psite[:, ii])
# Checks
deltanorm = np.abs(1 - np.sum(np.dot(vec, vec.conj())))
imagpart = np.sum(np.abs(np.imag(vec)))
if(len(vec) != locdim): raise Exception('Mismatch locdim.')
if(deltanorm > 1e-14): raise Exception('Norm violation')
if(imagpart != 0.0): raise Exception('Subspace has imag part.')
for jj in range(locdim):
tstr += '%15.12f\n'%(vec[jj])
fh.write(tstr)
# New left incoming must be right outgoing
Qleft[:] = Qrigh[:]
# The Lambdas would go here, for a product state they
# are all equal to 1, but we do not write them.
fh.close()
return
[docs]class ProductState():
"""
Create an unnormalized Product (or Fock) state. State will be normalized
upon writing. This class can currently not take care of symmetries.
**Arguments**
state_matrix : numpy 2d array
The number of sites in the system is defined through the number of
columns and the local dimension (equal for each site) through the
number of rows. In conclusion each column defines a local state
of the global Product/Fock state.
"""
def __init__(self, state_matrix):
self.state_matrix = state_matrix
[docs] def write(self, filename, enforce_cmplx=False):
"""
Write state to file
**Arguments**
filename : string
Filename for saving the state to file.
enforce_cmplx : Bool, optional
Enforce writting a complex MPS even if all entries in the state
are real which is necessary when complex Hamiltonians are passed.
Default to `False`.
"""
# Check that the ending is not `.bin`.
if(filename[-4:] == '.bin'):
raise ValueError("fock: extension `.bin` reserved for binary files!")
if(np.any(np.iscomplex(self.state_matrix)) or enforce_cmplx):
self.write_complex(filename)
else:
self.write_real(filename)
return
[docs] def write_real(self, flnm):
"""
Write the unnormalized Fock state with completely real entries.
**Arguments**
flnm : string
Filename for saving the state to file.
"""
# Open the file
fh = open(flnm, 'w+')
# Get the local dimension and the number of sites and string for dimensions
ld, ll = self.state_matrix.shape
dstr = "3 \n"
dstr += "1 " + str(ld) + " 1 "
dstr += "1 " + str(ld) + " 1 "
dstr += "1 2 3 \n"
# Write real/complex MPS/qMPS, system size and orthogonality center
fh.write("T \n")
fh.write("F \n")
fh.write(str(ll) + " \n")
fh.write("1 \n")
# Write the canonization toghether with haslambda
fh.write("c F \n")
for ii in range(ll - 1):
fh.write("r F \n")
fh.write("F \n")
# Write all the tensors
for ii in range(ll):
# Normalize the (local) state
self.state_matrix[:, ii] /= np.sqrt(np.sum(self.state_matrix[:, ii]**2))
# Enable check if seems to be necessary
#if(np.abs(np.sum(self.state_matrix[:, ii]
# * np.conj(self.state_matrix[:, ii])) - 1.0) > 1e-13):
# print(np.sum(self.state_matrix[:, ii]
# * np.conj(self.state_matrix[:, ii])))
# raise ValueError("Unnormalized state in Fock state!")
#else:
# print(np.sum(self.state_matrix[:, ii]
# * np.conj(self.state_matrix[:, ii])))
# Write dimensions and entries
fh.write(dstr)
for jj in range(ld):
fh.write(str(self.state_matrix[jj, ii]) + " \n")
fh.write("\n")
fh.close()
return
[docs] def write_complex(self, flnm):
"""
Write an unnormalized Fock state with complex entries.
**Arguments**
flnm : string
Filename for saving the state to file.
"""
# Open the file
fh = open(flnm, 'w+')
# Get the local dimension and the number of sites and string for dimensions
ld, ll = self.state_matrix.shape
dstr = "3 \n"
dstr += "1 " + str(ld) + " 1 "
dstr += "1 " + str(ld) + " 1 "
dstr += "1 2 3 \n"
# Write real/complex MPS/qMPS, system size and orthogonality center
fh.write("F \n")
fh.write("F \n")
fh.write(str(ll) + " \n")
fh.write("1 \n")
# Write the canonization toghether with haslambda
fh.write("c F \n")
for ii in range(ll - 1):
fh.write("r F \n")
fh.write("F \n")
# Write all the tensors
for ii in range(ll):
# Normalize the (local) state
self.state_matrix[:, ii] /= np.sqrt(np.sum(np.real(self.state_matrix[:, ii] * np.conj(self.state_matrix[:, ii]))))
# Enable check if it seems to be necessary
#if(np.abs(np.sum(self.state_matrix[:, ii]
# * np.conj(self.state_matrix[:, ii])) - 1.0) > 1e-13):
# print(np.sum(self.state_matrix[:, ii]
# * np.conj(self.state_matrix[:, ii])))
# raise ValueError("Unnormalized state in Fock state!")
#else:
# print(np.sum(self.state_matrix[:, ii]
# * np.conj(self.state_matrix[:, ii])))
# Write dimensions and entries
fh.write(dstr)
for jj in range(ld):
fh.write("(" + str(np.real(self.state_matrix[jj, ii])) + "," +
str(np.imag(self.state_matrix[jj, ii])) + ") \n")
fh.write("\n")
fh.close()
return
[docs]class PartlyProduct():
"""
Define an input state build of pure states possibly defined on several
lattice sites. Use the iadd (`+=`) to add vectors and
:py:func:`mps_state.write` to write your state to a file. This class can
currently not take care of symmetries.
**Arguments**
ld : int
local dimension of the Hilbert space
**Variables**
ld : int
local dimension of the Hilbert space
purestates : list of numpy 1d array
the list contains state vectors possibly spanning multiple sites.
The order in which they are added corresponds the first (left) site
to the last (right) site.
cutoff : float
cut singular values smaller than `cutoff`.
default to 1e-14
cutchi : int
cut bond dimension above `cutchi`
default to 500
**Details**
Where the Fock state (:math:`| \psi_1 \\rangle \otimes
| \psi_2 \\rangle \otimes \\ldots \otimes | psi_L \\rangle`) this kind of
state can handle entangled sites, e.g. in the case of qubits up to
20 neighboring sites entangled with each other:
:math:`| \psi_{1,2,3} \\rangle \otimes | \psi_{4} \\rangle \otimes \\ldots
| \psi_{L-1, L} \\rangle`
"""
def __init__(self, ld):
self.ld = ld
self.purestates = []
self.cutoff = 1e-14
self.cutchi = 500
return
[docs] def get_entry_real(self, entry):
"""
Return the (real-valued) entry of a tensor as a string (internal use).
**Arguments**
entry : float
entry of a tensor, real-valued
"""
return str(entry) + " \n"
[docs] def get_entry_complex(self, entry):
"""
Return the entry of a tensor as a string of a complex number (internal
use).
**Arguments**
entry : float, complex float
entry of a tensor
"""
return "(" + str(np.real(entry)) + ", " + str(np.imag(entry)) + ") \n"
[docs] def __iadd__(self, purestate):
"""
Add a state list to the list.
**Arguments**
purestate : numpy 1d array
state spanning possible multiple sites. If :math:`n` is the
number of sites spanned, the shape with :math:`l_d` the local
dimension should be :math:`(l_d)^n`.
"""
self.purestates.append(purestate)
return self
[docs] def write(self, flnm, enforce_complex=False):
"""
Write the initial mps state to a file.
**Arguments**
flnm : string
filename for the initial mps state.
enforce_cmplx : Bool, optional
Enforce writting a complex MPS even if all entries in the state are
real which is necessary when complex Hamiltonians are passed.
Default to `False`.
"""
# Check that the ending is not `.bin`.
if(flnm[-4:] == '.bin'):
raise ValueError("fock: extension `.bin` reserved for binary files!")
# Check for complex entries
is_complex = self.has_complex_state() or enforce_complex
if(is_complex):
self.get_entry = self.get_entry_complex
self.get_complex = lambda : "F \n"
else:
self.get_entry = self.get_entry_real
self.get_complex = lambda : "T \n"
# Decompose states
tensors = self.decompose()
#self.check_isometry(tensors)
#raise ValueError("Manual error")
# write actual file
self.write_file(tensors, self.ld, flnm)
return
[docs] def has_complex_state(self):
"""
Check for complex states in the list. Return ``True`` if complex numbers
occur, else ``False`` (internal use).
"""
for state in self.purestates:
if(np.any(np.iscomplex(state))): return True
return False
[docs] def check_isometry(self, tensors):
"""
Check isometry and normalization condition
**Arguments**
tensors : list of numpy rank-3 arrays
containing the tensors for the initial MPS with orthogonality
center on the first site.
"""
first = True
for tensor in tensors:
d1, d2, d3 = tensor.shape
if(first):
tmp = np.dot(np.reshape(tensor, [d1 * d2, d3]),
np.conj(np.transpose(np.reshape(tensor, [d1 * d2, d3]))))
print(np.trace(tmp))
first = False
else:
tmp = np.dot(np.reshape(tensor, [d1, d2 * d3]),
np.conj(np.transpose(np.reshape(tensor, [d1, d2 * d3]))))
print(np.sum(np.abs(tmp - np.eye(tmp.shape[1]))))
[docs] def decompose(self):
"""
Decompose pure states spanning multiple sites into tensors. Returned
is a list with the local tensors (internal use).
"""
# shortcut
ld = self.ld
# Save tensors here
mps_tensors = []
for state in self.purestates:
# Check if the state is only defined on one site
if(state.shape[0] == ld):
norm = np.sum(np.conj(state) * state)
mps_tensors.append(np.reshape(state / np.sqrt(norm),
[1, ld, 1]))
continue
# Find out how many sites the state spans
nsites = int(np.round(ma.log(state.shape[0], ld)))
if(ld**nsites != state.shape[0]): raise ValueError
# Initial dimension for first SVD
d1 = ld**(nsites - 1)
d2 = ld
# position to insert the tensors
pos = len(mps_tensors)
# initialize u
u = deepcopy(state)
for ii in range(nsites - 1):
[u, s, v] = la.svd(np.reshape(u, [d1, d2]),
full_matrices=False)
# cutback
s = s / np.sqrt(np.sum(s**2))
idx = np.argmax(s < self.cutoff)
idx = s.shape[0] if(idx == 0) else idx
if(idx < self.cutchi):
u = u[:, :idx]
s = s[:idx]
v = v[:idx, :]
else:
cut = np.sum(s[self.cutchi:]**2)
u = u[:, :self.cutchi]
s = s[:self.cutchi]
v = v[:self.cutchi, :]
# insert new tensor
mps_tensors.insert(pos, np.reshape(v, [s.shape[0], ld, d2 // ld]))
# normalize and contract singular values
s = s / np.sqrt(np.sum(s**2))
u = np.dot(u, np.diag(s))
# set dimensions
d2 = s.shape[0] * ld
d1 = d1 // ld
mps_tensors.insert(pos, np.reshape(u, [1, ld, s.shape[0]]))
return mps_tensors
[docs] def write_file(self, ts, ld, flnm):
"""
Open file handle and write the state (internal use).
**Arguments**
ts : list contain 3d numpy arrays
representing the normalized ond orthogonalized tensors
for the initial state.
ld : integer
local dimension of the sites
filename : string
Filename for saving the state to file.
"""
# Open the file
fh = open(flnm, 'w+')
# Get the number of sites
ll = len(ts)
# Write real/complex MPS/qMPS, system size and orthogonality center
fh.write(self.get_complex())
fh.write("F \n")
fh.write(str(ll) + " \n")
fh.write("1 \n")
# Write the canonization toghether with haslambda
fh.write("c F \n")
for ii in range(ll - 1):
fh.write("r F \n")
fh.write("F \n")
# Write all the tensors
for local in ts:
fh.write("3 \n")
fh.write(str(local.shape[0]) + " " + str(local.shape[1]) + " " +
str(local.shape[2]) + " ")
fh.write(str(local.shape[0]) + " " + str(local.shape[1]) + " " +
str(local.shape[2]) + " ")
fh.write("1 2 3 \n")
for i3 in range(local.shape[2]):
for i2 in range(local.shape[1]):
for i1 in range(local.shape[0]):
fh.write(self.get_entry(local[i1, i2, i3]))
fh.write("\n")
fh.close()
return
[docs]class State():
"""
This class represents a State object on L sites. (dj: Apparantly it is
not used in the library and it is unclear to me how to use it.)
**Arguments**
L : int
Number of sites for the state.
**Variables**
data : dictionary
Dictionary contains keys for system size ``L`` and ``states``
and ``weights``. The latter two are represented as lists.
"""
def __init__(self,L):
self.data={}
self.data['L']=L
self.data['states']=[]
self.data['weights']=[]
[docs] def addProductState(self,state,weight=1.0):
"""
Add a product state to the State object. If weight is supplied, the
product state is weighted by this factor.
"""
if isinstance(state,np.ndarray):
state=map(float,state)
if len(state)!=self.data['L']:
raise Exception("Length of product state added to State "
+ "disagrees with L used to initialize the state!")
self.data['states'].append(state)
self.data['weights'].append(weight)
[docs] def normalizeWeights(self):
"""
Normalize the weights of a State, in effect normalizing the state.
"""
norm=0.0
for weight in self.data['weights']:
norm=weight*weight
for weight in self.data['weights']:
weight=weight/sqrt(norm)
[docs] def writeqState(self,fileStub,qOperators):
"""
Translate the states in the State to the qState representation
using the qOperators structure, and then write to a file.
"""
Statefilename=fileStub+'qState.dat'
Statefile=open(Statefilename,'w')
self.normalizeWeights()
Statefile.write('%16i'%(len(self.data['weights']))+'\n')
count=0
for state in self.data['states']:
for i in range(self.data['L']):
Statefile.write('%30.15E'%(self.data['weights'][count])+'\n') #write weight
# state[i]=np.dot(qOperators['transformation'],state[i]) how to transform state? state[i]->
runningj=0
for jindex in qOperators['q_numbers']:
dj=qOperators['q_d'][jindex]
print('rj dj si', runningj, dj, state[i])
if runningj<=state[i] and state[i]<=runningj+dj:
tstr=''
for j in jindex:
tstr+='%16i'%(j)
Statefile.write(tstr+'\n')#write quantum number tuple
Statefile.write('%16i'%(dj)+'\n') #write deg dim of current irrep
#dum=np.zeros(dj)
Statefile.write('%16i'%(state[i]-runningj+1)+'\n')#write position within degeneracy dimension
runningj=runningj+dj
#dum[state[i]-runningj+1]=1.0
#tstr=''
#for j in dum:
# tstr+='%30.15E'%(j)
# qn[count].append(jindex)
# block[count].append(np.zeros(dj))
# block[count][i][state[i]-runningj+1]=self.data['weights'][count][i]
count+=1
Statefile.close()
return Statefilename
[docs] def writeState(self,fileStub):
"""
Write the state represented by the State object to a file.
"""
Statefilename=fileStub+'State.dat'
Statefile=open(Statefilename,'w')
self.normalizeWeights()
Statefile.write('%16i'%(len(self.data['weights']))+'\n')
count=0
for state in self.data['states']:
Statefile.write('%30.15E'%(self.data['weights'][count])+'\n')
for i in range(self.data['L']):
Statefile.write('%16i'%(state[i])+'\n')#write position and weight
count+=1
Statefile.close()
return Statefilename