"""exactdiag.py
Module defines a pure state and a density matrix based on a general abstract
quantum system. It defines time evolution methods for those systems.
"""
import sys
import numpy as np
#import numpy.linalg as nla
#import scipy.linalg as sla
#import scipy.sparse as sp
#import scipy.sparse.linalg as spla
from copy import deepcopy
from math import floor
import subprocess
import multiprocessing as mup
from itertools import product as Cartesian
from itertools import combinations
from time import sleep
if(sys.version_info > (3,6)):
from time import perf_counter as clock
else:
from time import clock
from os import remove as remove_file
#import MPSPyLib as mps
#import MPSPyLib.QILib as qi
if(int(np.version.version.split('.')[1]) > 8):
linspace = np.linspace
else:
def linspace(*args, **kwargs):
"""
Provide linspace function which allows dtype arguments for
older python versions.
"""
if('dtype' in kwargs):
dtype = kwargs['dtype']
del kwargs['dtype']
else:
dtype = np.float64
return np.array(np.linspace(*args, **kwargs), dtype=dtype)
get_spectrum = False
[docs]def set_spectrum(val):
"""
Set the value for the calculation of the spectrum.
**Arguments**
val : boolean
For calculating the spectrum need ``True``, otherwise ``False``.
"""
global get_spectrum
get_spectrum = val
return
# To-do list
# ----------
# _QuantumSystem.dot has bad implementation refering to self.rho
# Commented part inside trotter_layer
# Comment function line 1730 _get_bond_entropy_symm_sub
# get_bond_entropies for density matrices.
# fidelity of density matrices
[docs]def runED(params, Operators, Hmpos, sparse=True, maxdim=None):
"""
Run all simulation in the list built of parameter dictionaries.
**Arguments**
params : list
list of dictionaries where each dictionaries represents a MPS
simulation.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
Hmpos : instance of :py:class:`MPO.MPO` or list of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
sparse : bool, optional
Flag on the use of sparse matrices.
Default to True (use sparse matrices).
maxdim : None or int, optional
Maximal dimension for the complete Hilbert space in order to
prevent to large systems for ED.
Default to None (no limitation).
"""
if(hasattr(Hmpos, '__len__')):
get_mpo = lambda ii : Hmpos[ii]
else:
get_mpo = lambda ii : Hmpos
for ii in range(len(params)):
qtid = 0 if('QT' not in params[ii]) else params[ii]['QT']
if(qtid == 0):
if(params[ii]['verbose'] > 0):
print('Executing now ' + params[ii]['Write_Directory'] +
params[ii]['job_ID'] + params[ii]['unique_ID'])
ed_simulation(get_mpo(ii), Operators, params[ii], sparse, maxdim, 0)
else:
for jj in range(qtid):
if(params[ii]['verbose'] > 0):
print('Executing now ' + params[ii]['Write_Directory']
+ params[ii]['job_ID'] + params[ii]['unique_ID']
+ ', QT=' + str(jj + 1))
ed_simulation(get_mpo(ii), Operators, params[ii], sparse,
maxdim, jj + 1)
return
[docs]def runED_threading(params, Operators, Hmpos, threads, sparse=True,
maxdim=None, tsleep=30.0):
"""
Allows to run ED simulations with threading on a single node.
**Arguments**
params : list
list of dictionaries where each dictionaries represents a MPS
simulation.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
Hmpos : instance of :py:class:`MPO.MPO` or list of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
sparse : bool, optional
Flag on the use of sparse matrices.
Default to True (use sparse matrices).
maxdim : None or int, optional
Maximal dimension for the complete Hilbert space in order to
prevent to large systems for ED.
Default to None (no limitation).
tsleep : float, optional
Time interval to check for finished jobs.
**Details**
The multiprocessing module should be taking care to obtain thread-safe
processes, but complications might arise from there.
"""
if(hasattr(Hmpos, '__len__')):
get_mpo = lambda Hmpos, ii : Hmpos[ii]
else:
get_mpo = lambda Hmpos, ii : Hmpos
jobs = []
started = [0] * len(params)
ended = [False] * len(params)
# Set initial indices/data
idx = 0
qidx = 0
qinner = False
qtid = 0 if('QT' not in params[idx]) else params[idx]['QT']
# Start initial set of jobs
for ii in range(min(len(params), threads)):
Ham = get_mpo(Hmpos, idx)
print('Executing now ' + params[idx]['Write_Directory']
+ params[idx]['job_ID'] + params[idx]['unique_ID'] + ' '
+ str(qidx))
jobs.append(mup.Process(target=ed_simulation,
args=(Ham, Operators, params[idx],
sparse, maxdim, qidx)))
jobs[-1].start()
if(qinner and (qidx + 1 == qtid)):
# Last trajectory, go to next simulation
qinner = False
qidx = 0
idx += 1
if(idx < len(params)):
qtid = 0 if('QT' not in params[idx]) else params[idx]['QT']
elif(qinner):
# Go to the next quantum trajectory
qidx += 1
elif((qtid == 0) or (qtid == 1)):
# No quantum trajectories (or exactly one)
idx += 1
if(idx < len(params)):
qtid = 0 if('QT' not in params[idx]) else params[idx]['QT']
else:
# Go from the first trajectory to the next
qinner = True
qidx += 1
# Wait for jobs to finish before starting new jobs
ii = 0
while(idx != len(params)):
"""status = ''
for jj in range(len(jobs)):
status += str(jobs[jj].isAlive())
if(ii == 0): print('status', status)"""
if(not jobs[ii].is_alive()):
Ham = get_mpo(Hmpos, idx)
print('Executing now ' + params[idx]['Write_Directory']
+ params[idx]['job_ID'] + params[idx]['unique_ID'] + ' '
+ str(qidx))
jobs[ii] = mup.Process(target=ed_simulation,
args=(Ham, Operators,
params[idx], sparse,
maxdim, qidx))
jobs[ii].start()
if(qinner and (qidx + 1 == qtid)):
# Last trajectory, go to next simulation
qinner = False
qidx = 0
idx += 1
if(idx < len(params)):
qtid = 0 if('QT' not in params[idx]) else params[idx]['QT']
elif(qinner):
# Go to the next quantum trajectory
qidx += 1
elif((qtid == 0) or (qtid == 1)):
# No quantum trajectories (or exactly one)
idx += 1
if(idx < len(params)):
qtid = 0 if('QT' not in params[idx]) else params[idx]['QT']
else:
# Go from the first trajectory to the next
qinner = True
qidx += 1
ii = (ii + 1) % len(jobs)
if(ii == 0): sleep(tsleep)
# Wait for all jobs to finish
for ii in range(len(jobs)):
jobs[ii].join()
return
[docs]def ed_simulation(Hmpo, Operators, param, sparse, maxdim, qtid):
"""
Carry out a single static and dynamics simulations for one single set of
parameters.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
dictionary representing a MPS simulation.
sparse : bool
Flag on the use of sparse matrices.
maxdim : int
Maximal dimension for the complete Hilbert space in order to
prevent to large systems for ED.
qtid : integer
index for multiple realizations of the same simulations, e.g. with
quantum trajectories.
"""
tic = clock()
# Setup for symmetries
if('SymmSec' in param):
if(isinstance(param['SymmSec'], SymmetrySector)):
SymmSec = param['SymmSec']
else:
# Assume string leading to filename
SymmSec = SymmetrySector(None, None, sectorbase=param['SymmSec'])
else:
SymmSec = SymmetrySector(Operators, param)
if(SymmSec.sectordic is None): SymmSec = None
# Statics
if(len(param['timeevo_mps_initial']) > 0):
Psi = np.load(param['timeevo_mps_initial'] + '.npy')
if(len(Psi.shape) == 1):
Psi = PureState(Psi, Operators['I'].shape[0], param['L'], SymmSec)
elif(len(Psi.shape) == 2):
Psi = DensityMatrix(Psi, Operators['I'].shape[0], param['L'],
SymmSec)
else:
raise Exception('Rank-3 or higher for quantum state')
else:
Psi = statics_simulation(Hmpo, Operators, param, SymmSec, sparse,
maxdim)
Psi0 = deepcopy(Psi)
# Dynamics
if('Quenches' in param):
nlind = len(Hmpo.data['lind1']) \
+ len(Hmpo.data['MBSLind']) \
+ len(Hmpo.data['lindspec']) \
+ len(Hmpo.data['lindlspec'])
qt = (nlind > 0) and (qtid != 0)
if(qt):
np.random.seed(qtid)
elif((nlind > 0) and isinstance(Psi, PureState)):
# Convert to density matrix for open system evolution
Psi = Psi.densitymatrix()
dynamic_simulation(Psi, Psi0, Hmpo, Operators, param, sparse,
maxdim, qtid)
# Write calculation time to log-file
toc = clock()
fh = open(param['Output_Directory'] + param['job_ID'] +
param['unique_ID'] + '.log', 'w+')
fh.write('time taken:%30.15e seconds!\n'%(toc - tic))
fh.close()
# Delete temporary files
cmd = ['rm', '-rf', param['Write_Directory'] + param['job_ID']
+ param['unique_ID'] + '*']
cmd = subprocess.list2cmdline(cmd)
cmd = subprocess.call(cmd, shell=True)
return
# ==============================================================================
# METHODS FOR DYNAMICS
# ==============================================================================
[docs]def dynamic_simulation(Psi, Psi0, Hmpo, Operators, param, sparse, maxdim, qtid):
"""
Run the dynamics for a single parameter set.
**Arguments**
Psi : Pure state or density matrix.
Quantum state for time evolution. It can either be represented as
instance of :py:class:`exactdiag.PureState` or as instance of
:py:class:`exactdiag.DensityMatrix`
Psi0 : instance of :py:class:`exactdiag.PureState`
Quantum state to calculate overlap with the state at every measurement.
Meant to be used with the state at time t=0 for the Loschmidt echo.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
dictionary representing a MPS simulation.
sparse : bool
Flag on the use of sparse matrices.
maxdim : int
Maximal dimension for the complete Hilbert space in order to
prevent to large systems for ED.
qtid : integer
If > 0, Psi as an instance of :py:class:`exactdiag.PureState` is
evolved with quantum trajectories (open system). If 0, instance of
:py:class:`exactdiag.PureState` is evolved under the Schrodinger
equation. If Psi is an instance of :py:class:`exactdiag.DensityMatrix`,
it is evolved in the Liouville space with the Lindblad master equation
and the flag qtid has no effect.
"""
if(param['Quenches'].classname() == 'QuantumCircuits'):
dynamic_circuits(Psi, Psi0, Operators, param)
return
qt = (qtid != 0)
ext = '_qt%09d.dat'%(qtid) if(qt) else '.dat'
for ii in range(len(param['Quenches'].data)):
dt = param['Quenches'].data[ii]['deltat']
Time = param['Quenches'].data[ii]['time']
nsteps = int(floor(Time / dt))
extradt = Time - nsteps * dt
if(extradt <= 0.000001): extradt = None
sfo = param['Quenches'].data[ii]['stepsforoutput']
time = 0.0
for jj in range(nsteps):
Psi.propagate(dt, Hmpo, Operators, param, sparse, maxdim,
(ii, time + 0.5 * dt), qt)
time = time + dt
if((jj + 1)%sfo == 0):
# Small measurements (energy, Loschmidt echo)
en = np.real(Psi.meas_mpo(Hmpo, Operators, param,
(ii, time), maxdim))
le = Psi.fidelity(Psi0)
ler = np.real(le)
lec = np.imag(le)
pur = Psi.purity()
Obsfile = param['Output_Directory'] + param['job_ID'] \
+ param['unique_ID'] + 'ObsOutDynamics/Dynamics_' \
+ str(ii + 1) + 'step_%06d'%(jj + 1) + ext
fh = open(Obsfile, 'w+')
# Bond dimension and error set to zero
fh.write('%30.15e%30.15e%16d%16d%30.15e%30.15e%30.15e\n'%(
time, en, 0, 0, ler, lec, pur))
Psi.meas_observables(Operators, param, fh,
'DynamicsObservables')
fh.close()
if(not extradt is None):
Psi.propagate(dt, Hmpo, Operators, param, sparse, maxdim,
(ii, time + 0.5 * extradt), qt)
time = time + extradt
# Small measurements (energy, Loschmidt echo)
en = np.real(Psi.meas_mpo(Hmpo, Operators, param,
(ii, time), maxdim))
le = Psi.fidelity(Psi0)
ler = np.real(le)
lec = np.imag(le)
pur = Psi.purity()
Obsfile = param['Output_Directory'] + param['job_ID'] \
+ param['unique_ID'] + 'ObsOutDynamics/Dynamics_' \
+ str(ii + 1) + 'step_extra' + ext
fh = open(Obsfile, 'w+')
# Bond dimension and error set to zero
fh.write('%30.15e%30.15e%16d%16d%30.15e%30.15e%30.15e\n'%(
time, en, 0, 0, ler, lec, pur))
Psi.meas_observables(Operators, param, fh, 'DynamicsObservables')
fh.close()
# Overwrite last propagator in any case
Psi.lastprop = None
# Overwrite subblocks of same symmetry in any case
Psi.bidx = None
return
[docs]def dynamic_circuits(Psi, Psi0, Operators, param, QCin=None):
"""
Evolve the quantum system with circuits.
**Arguments**
Psi : Pure state or density matrix.
Quantum state for time evolution. It can either be represented as
instance of :py:class:`exactdiag.PureState` or as instance of
:py:class:`exactdiag.DensityMatrix`
Psi0 : instance of :py:class:`exactdiag.PureState`
Quantum state to calculate overlap with the state at every measurement.
Meant to be used with the state at time t=0 for the Loschmidt echo.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
dictionary representing a MPS simulation.
caller : instance of QuantumCircuits or None, optional
If QCin is not ``None``, measurements are skipped and the final state
is returned.
Default to ``None``.
"""
if(QCin is None):
QC = param['Quenches']
else:
QC = QCin
# Extension (never quantum trajectories)
ext = '.dat'
time = 0.0
jj = 0
for QCsub in QC:
# Convert to rank-L tensor (or 2L for density matrix)
Psi._trotter_initialize(False)
for gate, sites, timeii in QCsub:
if(not Psi.SymmSec is None):
# Check if bidx indices are correct (number of sites)
if(not Psi.bidx is None):
for key in Psi.bidx.keys():
nn = Psi.bidx[key][0]
nn = len(nn)
break
# Have to delete mapping - wrong number of sites
if(nn != len(sites)): Psi.bidx = None
if(Psi.bidx is None):
Psi.bidx = idx_symmetry_blocks(Operators, param,
ll=len(sites))
Psi.apply_op(gate, sites)
time += timeii
# Convert back to state or density matrix
Psi._trotter_finalize(False, None, None, None)
jj += 1
if(QCin is None):
# Small measurements (energy, Loschmidt echo)
en = 0.0#np.real(Psi.meas_mpo(Hmpo, Operators, param,
# (ii, time), maxdim))
le = Psi.fidelity(Psi0)
ler = np.real(le)
lec = np.imag(le)
pur = Psi.purity()
Obsfile = param['Output_Directory'] + param['job_ID'] \
+ param['unique_ID'] + 'ObsOutDynamics/Dynamics_' \
+ str(1) + 'step_%06d'%(jj) + ext
fh = open(Obsfile, 'w+')
# Bond dimension and error set to zero
fh.write('%30.15e%30.15e%16d%16d%30.15e%30.15e%30.15e\n'%(
time, en, 0, 0, ler, lec, pur))
Psi.meas_observables(Operators, param, fh,
'DynamicsObservables')
fh.close()
if(QCin is None):
return
else:
return Psi
[docs]def expm_symm_tridiag(diag, sdiag, sc):
"""
Calculate exponential of symmetric tridiagonal matrix.
**Arguments**
diag : 1d numpy array
diagonal entries (dimension of the matrix)
sdiag : 1d numpy array
supdiagonal entries
sc : float, complex
scaling inside of the exp(.)
"""
dtype = np.float64 if(np.complex(sc) == 0.0) else np.complex128
if(diag.shape[0] == 1):
# One entry
prop = np.zeros((1, 1), dtype=dtype)
prop[0, 0] = np.exp(sc * diag[0])
return prop
arr = np.zeros((2, diag.shape[0]), dtype=dtype)
arr[0, 1:] = sdiag
arr[1, :] = diag
vals, vecs = sla.eig_banded(arr)
prop = np.dot(vecs, np.dot(np.diag(np.exp(sc * vals)),
np.conj(np.transpose(vecs))))
return prop
# ==============================================================================
# METHODS FOR SYMMETRIES
# ==============================================================================
[docs]class SymmetrySector():
"""
Taking care of the symmtry. Generates all basis states in the Hilbert space
and returns the indices of the states which are conform with the quantum
numbers defined for the simulations. The result returned is a list of
integers, unless there are no symmetries defined (None).
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
Contains all operators of the simulation, especially the diagonal
generators of the symmetry.
param : dictionary
Contains the simulation setup.
sectorbase : 2d numpy array, optional
If basis states are calculated manually, a 2d numpy array with
a basis state in each row can be passed to initialize the symmetry
sector. Number of rows is the number of basis states, number of columns
corresponds to the number of sites. For ``Operators`` and ``param`` are
not referenced and ``None`` can be passed to them.
Default to None.
empty : bool, optional
For internal use, no variables are set. Allows for setting them manually.
Default to False
**Variables**
len : int
number of basis states in the symmetry sector.
sectorbase : 2d numpy array, str
local basis for each basis state of the symmetry sector.
If str, filename of the 2d numpy array to be loaded.
sectordic : dict
For a tuple with the local basis state the corresponding
index is returned. Does not handle key errors.
empty : bool, optional
for internal use only.
Default to False
"""
def __init__(self, Operators, param, sectorbase=None, empty=False):
# Manual setup of SymmetrySector
if(empty): return
# Use given sectorbase for setup
if(not (sectorbase is None)):
if(isinstance(sectorbase, str)):
sectorbase = np.load(sectorbase + '.npy')
self.len = sectorbase.shape[0]
self.sectorbase = sectorbase
self.sectordic = {}
for ii in range(self.len):
self.sectordic[tuple(sectorbase[ii, :])] = ii
return
# Count number of symmetries and check for fast return
try:
U1Gen = param['Abelian_generators']
U1sector = param['Abelian_quantum_numbers']
nU1 = len(U1Gen)
except KeyError:
nU1 = 0
if("Abelian_generators" in param):
U1Gen = param['Abelian_generators']
U1sector = param['Abelian_quantum_numbers']
nU1 = len(U1Gen)
else:
nU1 = 0
try:
Z2Gen = param['Discrete_generators']
Z2sector = param['Discrete_quantum_numbers']
nZ2 = len(Z2Gen)
except KeyError:
nZ2 = 0
nq = nU1 + nZ2
if(nq == 0):
self.sectorbase = None
self.sectordic = None
return
ll = param['L']
ld = Operators['I'].shape[0]
sectorbase_ = []
#self.sectoridx = []
for ii in range(ld**ll):
base = np.base_repr(ii, base=ld)
base = list(map(int, list('0' * (ll - len(base)) + base)))
insector = True
for jj in range(nU1):
q = 0
Op = np.diag(Operators[U1Gen[jj]])
for kk in range(ll):
q += Op[base[kk]]
if(abs(q - U1sector[jj]) > 1e-13):
insector = False
break
if(insector):
for jj in range(nZ2):
q = 0
Op = np.diag(Operators[Z2Gen[jj]])
for kk in range(ll):
q += Op[base[kk]]
if(abs(q%2 - Z2sector[jj]) > 1e-13):
insector = False
break
if(insector):
sectorbase_.append(np.array(base, dtype=np.int8))
#self.sectoridx.append(ii)
self.sectordic = {}
self.sectorbase = np.zeros((len(sectorbase_), sectorbase_[0].shape[0]),
dtype=np.int8)
self.len = len(sectorbase_)
for ii in range(self.len):
self.sectordic[tuple(sectorbase_[ii])] = ii
self.sectorbase[ii, :] = sectorbase_[ii]
return
[docs] def __getitem__(self, key):
"""
Getitem is defined to return the index of a basis state.
**Arguments**
key : list, array (1d)
contains the local basis for each site
"""
return self.sectordic[tuple(key)]
[docs] def __len__(self):
"""
Returns length of SymmetrySector as number of basis states.
"""
return self.len
[docs] def permute(self, perm):
"""
Permute sites according to permutation.
**Arguments**
perm : list of ints
Permutation vector.
"""
return self.sectorbase[:, perm]
[docs] def __call__(self, ii):
"""
Return the indices of the local basis for the ii-th state.
**Arguments**
ii : int
Consider to ii-th state in the symmetry-sector based basis.
"""
return self.sectorbase[ii, :]
[docs] def get_subspace(self, idxs):
"""
Return the symmetry adapted basis of a subspace.
**Arguments**
idxs : list of integers
these are the sites to be considered for the new space.
"""
hashtab = {}
mapping = {}
jj = 0
for ii in range(self.sectorbase.shape[0]):
if(tuple(self.sectorbase[ii, idxs]) not in hashtab):
hashtab[tuple(self.sectorbase[ii, idxs])] = jj
mapping[jj] = ii
jj += 1
SubSec = SymmetrySector(None, None, empty=True)
SubSec.len = jj
SubSec.sectordic = hashtab
SubSec.sectorbase = np.zeros((jj, len(idxs)), dtype=np.int8)
for ii in range(jj):
SubSec.sectorbase[ii, :] = self.sectorbase[mapping[ii], idxs]
return SubSec
[docs] def write(self, filename):
"""
Save the sectorbase of the Symmetry class to hard disk.
This is all information necessary to recreate it.
**Arguments**
filename : str
Basis string for the file name extended with ``.npy``
for the sectorbase matrix.
"""
np.save(filename + '.npy', self.sectorbase)
return
[docs]def walls_particles(ll, fill):
"""
The distribution of fill particles on ll sites can be described
via fill particles and (ll - 1) walls. This iterator yields
all possible combinations where 0s represent a wall and 1s
represent the a particle.
**Arguments**
ll : int
number of sites in the system.
fill : int
number of particles distributed on the system.
"""
for elem in combinations(range(ll - 1 + fill), fill):
state = [0] * (ll - 1 + fill)
for subelem in elem:
state[subelem] = 1
yield state
[docs]def SingleU1Symmetry(Operators, param):
"""
Create symmetry sector for Bose-Hubbard type models with minimum
number of particles zero up to some maximum filling.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
contains all operators including the generator of the symmetry.
It is checked that the symmetry generator is the number operator.
param : dictionary
Contains the simulation setup. It is checked that only one
U(1) symmetry is present.
"""
try:
U1Gen = param['Abelian_generators']
U1sector = param['Abelian_quantum_numbers']
except KeyError:
raise Exception('U(1) symmetry has to be specified')
try:
Z2Gen = param['Discrete_generators']
if(len(Z2Gen) > 0): raise Exception('Z2 symmetry not allowed.')
except KeyError:
# Perfectly fine, we do not expect this key
pass
if(len(U1Gen) != 1): raise Exception('Not exactly one U(1) symmetry.')
# Check that generator is number operator
ld = Operators['I'].shape[0]
ref = np.cumsum(np.ones(ld)) - 1.0
if(np.sum(np.abs(np.diag(Operators[U1Gen[0]]) - ref)) > 1e-14):
raise Exception('Generator of U(1) not the number operator.')
ll = param['L']
fill = U1sector[0]
nmax = ld - 1
states = []
for elem in walls_particles(ll, fill):
fillii = 0
state = []
for ii in elem:
if(ii == 0):
state.append(fillii)
fillii = 0
else:
fillii += 1
state.append(fillii)
if(max(state) > nmax): continue
states.append(state)
dim = len(states)
sectorbase = np.zeros((dim, ll), dtype=np.int8)
for ii in range(dim):
sectorbase[ii, :] = np.array(states[ii], dtype=np.int8)
del states
return SymmetrySector(None, None, sectorbase=sectorbase)
# ==============================================================================
# METHODS FOR STATIC SIMULATIONS
# ==============================================================================
[docs]def statics_simulation(Hmpo, Operators, param, SymmSec, sparse, maxdim):
"""
Complete a static simulation including finding the ground and excited states
and measuring those. Returns the ground state if needed for dynamics.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian written in terms of MPO rules.
Operators : instance of the operators class :py:class:`Ops.OperatorList`
containing the operators to generate Hamiltonian.
param : dictionary with simulation parameters
settings of the whole simulation.
SymmSec : list or None
list of indices containing the states in the given basis which
are in the symmetry sector considered.
sparse : bool
flag if sparse matrix should be used for Hamiltonian
maxdim : int
limits the dimension of the complete Hilbert space.
"""
# Shortcut to iterator over states
statics_state = statics_state_sparse if(sparse) else statics_state_dense
# Get states
nn = param['n_excited_states']
Psis = []
ens = []
ll = param['L']
ld = Operators['I'].shape[0]
for elem in statics_state(Hmpo, param, nn, SymmSec, maxdim):
Psis.append(PureState(elem[0], ld, ll, SymmSec))
ens.append(elem[1])
# Carry out the measurements
# --------------------------
Obsfilebase = param['Output_Directory'] + param['job_ID'] + \
param['unique_ID'] + 'ObsOut' + str(0) + '_1.dat'
fh = open(Obsfilebase, 'w+')
fh.write('%30.15e%30.15e T%16d%30.15e\n'%(ens[0], 0.0, 0, -9999999.0))
Psis[0].meas_observables(Operators, param, fh, 'MPSObservables')
fh.close()
return Psis[0]
[docs]def statics_state_sparse(Hmpo, param, ii, SymmSec, maxdim):
"""
Find a specific state of the Hamiltonian. Returns state and energy (in this
order). Function for dense Hamiltonians.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian written in terms of MPO rules.
param : dictionary with simulation parameters
settings of the whole simulation.
ii : int, optional
index of the state to be returned. If 0, ground state is returned,
if 1 first excited state and so on. No specific handling for
degenerate states.
Default to 0 (ground state)
SymmSec : list or None
list of indices containing the states in the given basis which
are in the symmetry sector considered.
maxdim : int, optional
limits the dimension of the complete Hilbert space.
"""
Ham = Hmpo.build_hamiltonian(param, SymmSec, True, maxdim)
nn = min(4 * (ii + 1), Ham.shape[0] - 1)
vals, vecs = spla.eigsh(Ham, k=nn, which='SA')
for jj in range(ii + 1):
idx = vals.argsort()[jj]
yield (vecs[:, idx], vals[idx])
[docs]def statics_state_dense(Hmpo, param, ii, SymmSec, maxdim):
"""
Find a specific state of the Hamiltonian. Returns state and energy (in this
order). Function for dense Hamiltonians.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian written in terms of MPO rules.
param : dictionary with simulation parameters
settings of the whole simulation.
ii : int, optional
index of the state to be returned. If 0, ground state is returned,
if 1 first excited state and so on. No specific handling for
degenerate states.
Default to 0 (ground state)
SymmSec : list or None
list of indices containing the states in the given basis which
are in the symmetry sector considered.
maxdim : int, optional
limits the dimension of the complete Hilbert space.
"""
Ham = Hmpo.build_hamiltonian(param, SymmSec, False, maxdim)
vals, vecs = nla.eigh(Ham)
for jj in range(ii + 1):
idx = vals.argsort()[jj]
yield (vecs[:, idx], vals[idx])
# ==============================================================================
# METHODS FOR GENERAL QUANTUM SYSTEMS
# ==============================================================================
[docs]class _QuantumSystem():
"""
Abstract type for a quantum system, either pure state or density
matrix. Implements general methods equal for both.
"""
[docs] def traceout(self, trce, keep=None, enforce_dense=False):
"""
Carry out the partial trace on the system returning the density
matrix for the remaining system.
**Arguments**
trce : list of int
Site indices to trace out. Must be between 0 < ll - 1 and unique.
keep : list of int, optional
Site indices to be kept in reduced density matrix. Must be between
0 < ll - 1 and unique. Consistency with trce is not checked.
Default to None (indices are generated from trce)
enforce_dense : bool, optional
Return reduced density matrix on complete Hilbert space without
using symmetries (True). Otherwise symmetries might be used if
considered convenient.
Default to False.
"""
if(self.SymmSec is None):
# No symmetry - return rho on complete space
return self._traceout_nosymm(trce, keep)
if(enforce_dense):
# Complete space enforced
return self._traceout_symm2nosymm(trce, keep)
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
SubSec = self.SymmSec.get_subspace(keep)
if(isinstance(self.ld, int)):
if(self.ld**(self.ll - len(trce)) > len(SubSec)):
return self._traceout_symm(trce, keep, SubSec=SubSec)
else:
return self._traceout_symm2nosymm(trce, keep, SubSec=SubSec)
else:
if(np.prod(self.ld) / np.prod(self.ld[trce]) > len(SubSec)):
return self._traceout_symm(trce, keep, SubSec=SubSec)
else:
return self._traceout_symm2nosymm(trce, keep, SubSec=SubSec)
return
[docs] def reduced(self, keep, enforce_dense=False):
"""
Build the reduced density matrix of the system for a specified set
of sites.
**Arguments**
keep : list of int, optional
Site indices to be kept in reduced density matrix. Must be between
0 < ll - 1 and unique. Consistency with trce is not checked.
Default to None (indices are generated from trce)
enforce_dense : bool, optional
Return reduced density matrix on complete Hilbert space without
using symmetries (True). Otherwise symmetries might be used if
considered convenient.
Default to False.
"""
mask = np.ones((self.ll), dtype=np.bool)
mask[keep] = False
trce = list(np.array(list(range(self.ll)))[mask])
return self.traceout(trce, keep=keep, enforce_dense=enforce_dense)
[docs] def delete_rhos(self):
"""
Empty dictionary with reduced density matrices.
"""
self.rhos = {}
return
[docs] def getset_rho(self, key):
"""
Return the reduced density matrix from the complete wave function.
**Arguments**
key : tuple
containing one or two integers specifying the site indices.
"""
try:
return self.rhos[key]
except KeyError:
pass
self.rhos[key] = self.reduced(list(key), enforce_dense=True)
return self.rhos[key]
[docs] def dot(self, *args):
"""
General dot product split into different cases, which are detected
by the data type.
"""
if(isinstance(args[0], np.ndarray)):
if(len(args[0].shape) == 2):
# Assume operator
return np.dot(self.rho, args[0])
else:
raise Exception('dot: Shape not implemented')
else:
raise Exception('dot: unknown argument')
[docs] def meas_observables(self, Operators, param, fh, obstype):
"""
Measure the observables.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
All operators defined for the simulation.
param : dict
containing simulation setup, necessary for MPO measurement.
fh : filehandle
Open file. Results are written to this file handle.
obstype : string
Key to get Observables class in param.
"""
Obs = param[obstype]
ld = Operators['I'].shape[0]
for ii in range(len(Obs['site'])):
# Local measurements rho
Op = Operators[Obs['site']['Op'][ii]] * \
Obs.data['site']['weight'][ii]
local = ''
for jj in range(self.ll):
this = self.getset_rho((jj,))
local += '%30.15e'%(np.real(np.trace(this.dot(Op))))
fh.write(local + '\n')
if(Obs['SiteEntropy'].has):
# Set site entropy:
sitee = ''
for jj in range(self.ll):
Mat = self.getset_rho((jj,))
vals = sla.eigh(Mat.rho, eigvals_only=True)
vals = vals[vals > 1e-16]
sitee += '%30.15e'%(- np.sum(vals * np.log(vals)))
fh.write(sitee + '\n')
if(Obs['BondEntropy'].has):
# Set bond entropy
bonde = ''
for be in self.get_bond_entropies():
bonde += '%30.15e'%(be)
fh.write(bonde + '\n')
for ii in range(len(Obs['corr'])):
# Correlation measurements
Op1 = Operators[Obs.data['corr']['Op'][ii][0]] * \
Obs.data['corr']['weight'][ii]
Op2 = Operators[Obs.data['corr']['Op'][ii][1]]
Op = np.kron(Op1, Op2)
Ops = np.dot(Op1, Op2)
is_cmplx = (self.getset_rho((0,)).rho.dtype == np.complex64) or \
(Op.dtype == np.complex64)
is_cmplx = (is_cmplx and (not Obs['corr'].isherm[ii]))
if(is_cmplx):
fh.write('C\n')
else:
fh.write('R\n')
for jj in range(self.ll):
corrr = ''
corrc = ''
for kk in range(self.ll):
if(jj == kk):
Rho1 = self.getset_rho((jj,))
res = np.trace(Rho1.dot(Ops))
else:
Rho2 = self.getset_rho((jj, kk))
res = np.trace(Rho2.dot(Op))
corrr += '%30.15e'%(np.real(res))
corrc += '%30.15e'%(np.imag(res))
fh.write(corrr + '\n')
if(is_cmplx): fh.write(corrc + '\n')
for ii in range(len(Obs['fermicorr'])):
raise NotImplementedError('Fermi-correlation measures not' +
' implemented.')
for ii in range(len(Obs['string'])):
raise NotImplementedError('String measures not implemented.')
for ii in range(len(Obs['MPO'])):
maxdim = None
res = self.meas_mpo(Obs['MPO']['Op'][ii], Operators, param,
None, maxdim)
fh.write('%30.15e\n'%(np.real(res)))
if(Obs['DensityMatrix_i'].has_rhoi):
for ii in Obs['DensityMatrix_i'].ij:
# Single site density matrices
rho = self.getset_rho((ii - 1,))
strr = ''
strc = ''
for jj in range(ld):
for kk in range(ld):
# Comply with column-major format in fortran
strr += '%30.15e'%(np.real(rho.rho[kk, jj]))
strc += '%30.15e'%(np.imag(rho.rho[kk, jj]))
fh.write(strr + '\n')
if(rho.rho.dtype == np.complex128): fh.write(strc + '\n')
if(Obs['DensityMatrix_ij'].has_rhoij):
for i1, i2 in Obs['DensityMatrix_ij'].ij:
# Two site density matrices
rho = self.getset_rho((i1 - 1, i2 - 1))
# Permute for column-major order of local Hilbert spaces
tmp = np.transpose(np.reshape(rho.rho, [rho.ld] * 4),
[1, 0, 3, 2])
tmp = np.reshape(tmp, [rho.ld**2] * 2)
strr = ''
strc = ''
for jj in range(ld**2):
for kk in range(ld**2):
# Comply with column-major format in fortran
strr += '%30.15e'%(np.real(tmp[kk, jj]))
strc += '%30.15e'%(np.imag(tmp[kk, jj]))
fh.write(strr + '\n')
if(rho.rho.dtype == np.complex128): fh.write(strc + '\n')
if(len(Obs['DensityMatrix_red']) > 0):
for elem in Obs['DensityMatrix_red'].inds:
rho = self.getset_rho(tuple(elem))
# Permute for column-major order of local Hilbert spaces
# in sub-index of matrix
nelem = len(elem)
inda = np.linspace(0, nelem - 1, nelem, dtype=np.int32)[::-1]
indb = inda + nelem
ind = list(inda) + list(indb)
tmp = np.transpose(np.reshape(rho.rho,
[rho.ld] * (2 * nelem)), ind)
tmp = np.reshape(tmp, [rho.ld**nelem] * 2)
strr = ''
strc = ''
for jj in range(ld**nelem):
for kk in range(ld**nelem):
# Comply with column-major format in fortran
strr += '%30.15e'%(np.real(tmp[kk, jj]))
strc += '%30.15e'%(np.imag(tmp[kk, jj]))
fh.write(strr + '\n')
if(rho.rho.dtype == np.complex128): fh.write(strc + '\n')
if(Obs['MI'].has):
# Mutual information matrix
mim = np.zeros((self.ll, self.ll))
# entropy natural basis / base d
entrn = lambda ve : - np.sum(ve[ve > 1e-15] * np.log(ve[ve > 1e-15]))
for jj in range(self.ll):
# Diagonal element
tmp = self.getset_rho((jj,)).rho
mim[jj, jj] = entrn(nla.eigh(tmp)[0]) / np.log(ld)
for kk in range(jj + 1, self.ll):
tmp = self.getset_rho((jj, kk)).rho
mim[jj, kk] = entrn(nla.eigh(tmp)[0]) / np.log(ld)
for jj in range(1, self.ll):
for kk in range(jj + 1, self.ll):
mim[jj, kk] = 0.5 * (mim[jj, jj] + mim[kk, kk]
- mim[jj, kk])
mim[kk, jj] = mim[jj, kk]
for jj in range(self.ll):
strr = ''
for kk in range(self.ll):
strr += '%30.15e'%(mim[jj, kk])
fh.write(strr + '\n')
if(Obs.data['Lambda'].has):
# Singular values bipartition
raise NotImplementedError('Singular values of bipartitions.')
for ii in range(len(Obs['dist_psi'])):
# Distance measures to pure states (ED does not care if pure
# or mixed).
State = np.load(Obs['dist_psi']['term'][ii] + '.npy')
disttype = Obs['dist_psi']['dist'][ii]
fh.write('%30.15e\n'%(self.distance(State, disttype)))
del State
for ii in range(len(Obs['corr2nn'])):
raise NotImplementedError('Corr2NN not implemented.')
for ii in range(len(Obs['Gates'])):
# Gate measurements
Phi = deepcopy(self)
Phi = dynamic_circuits(Phi, Phi, Operators, param,
QCin=Obs['Gates'].Gates[ii])
over = self.fidelity(Phi) * Obs.data['Gates'].weight[ii]
fh.write('%30.15e%30.15e\n'%(np.real(over), np.imag(over)))
for ii in range(len(Obs['dist_rho'])):
# Distance measures mixed states / density matrices
State = np.load(Obs['dist_rho'].term[ii] + '.npy')
disttype = Obs['dist_rho'].dist[ii]
fh.write('%30.15e\n'%(self.distance(State, disttype)))
del State
if(Obs['State'].has):
self.save(fh.name[:-4] + '_State.npy')
self.delete_rhos()
return
[docs] def save(self, flnm):
"""
Save a state as numpy array. Must be implemented by classes
themselves.
"""
raise NotImplementedError('Method not implemented for this type')
[docs] def propagate(self, dt, Hmpo, Operators, param, sparse, maxdim,
quenchtime, qt):
"""
Time evolution of a quantum system for one time step with either of
the following three methods: matrix exponential, Krylov approximation
to the matrix exponential, or Trotter decomposition. Method is
determined via the convergence parameters.
**Arguments**
dt : float
time step (discretization of time).
Hmpo : instance of :py:class:`MPO.MPO`
Defines the Hamiltonnian of the system.
Operators : instance of :py:class:`Ops.OperatorList`
Contains the operators used to build the Hamiltonian.
param : dictionary
Setup for the simulation.
sparse : bool
Flag for the use of sparse methods (only referenced for matrix
exponential method and (possibly) Krylov).
maxdim : int
maximal dimension of the Hilbert space. Only referenced in ...
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
qt : boolean
If ``True``, this as an instance of :py:class:`exactdiag.PureState`
is evolved with quantum trajectories (open system). If ``False``,
instance of :py:class:`exactdiag.PureState` is evolved under the
Schrodinger equation. If Psi is an instance of a
:py:class:`exactdiag.DensityMatrix`, it is evolved in the Liouville
space with the Lindblad master equation and the flag qt has no
effect.
"""
ConvParam = param['Quenches'].data[0]['ConvergenceParameters']
if(isinstance(ConvParam, mps.KrylovConvParam)):
self.propagate_krylov(dt, Hmpo, Operators, param, sparse,
maxdim, quenchtime, qt)
elif(isinstance(ConvParam, mps.TEBDConvParam)):
self.propagate_trotter(dt, Hmpo, Operators, param, quenchtime, qt)
elif(isinstance(ConvParam, mps.ExpmConvParam)):
self.propagate_expm(dt, Hmpo, param, sparse, maxdim, quenchtime, qt)
else:
raise ValueError('Unknown time evolution method')
return
[docs] def propagate_trotter(self, dt, Hmpo, Operators, param, quenchtime, qt):
"""
Propagate the quantum system with a Trotter decomposition. This
function calls depending on the order of the decomposition the
:py:func:`exactdiag.PureState.trotter_layer` or
:py:func:`exactdiag.DensityMatrix.trotter_layer` for the odd/even
part.
**Arguments**
dt : float
size of the time step/discretization.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian as rule set for the Trotter decomposition.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
Contains the simulation setup.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
qt : bool
If ``True`` quantum states are evolved under quantum trajectories,
else (``False``) quantum states are evolved under the Schrodinger
equation. For DensityMatrix this should be always ``False``.
**Details**
The default order for the Trotter decomposition is the 2nd order.
The default is overwritten with the order in
:py:class:`convergence.TEBDConvParam` if this class is given. Possible
values are therefore 2 and 4.
"""
self._trotter_initialize(qt)
if(self.lastprop is None): self.lastprop = {}
# Set order of Trotter decomposition
if(isinstance(param['Quenches'].data[0]['ConvergenceParameters'],
mps.TEBDConvParam)):
order = param['Quenches'].data[quenchtime[0]]\
['ConvergenceParameters'].data[0]['order']
else:
order = 2
if(order == 2):
self.trotter_layer(0.5 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
self.trotter_layer(dt, Hmpo, Operators, param,
quenchtime, False, qt)
self.trotter_layer(0.5 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
elif(order == 4):
c3 = 1.0 / (2 * (4 - 4**(1.0 / 3)))
c1 = 1 - 8 * c3
c2 = 0.5 * (1 - 6 * c3)
c4 = 2 * c3
self.trotter_layer(c3 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
self.trotter_layer(c4 * dt, Hmpo, Operators, param,
quenchtime, False, qt)
self.trotter_layer(c4 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
self.trotter_layer(c4 * dt, Hmpo, Operators, param,
quenchtime, False, qt)
self.trotter_layer(c2 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
self.trotter_layer(c1 * dt, Hmpo, Operators, param,
quenchtime, False, qt)
self.trotter_layer(c2 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
self.trotter_layer(c4 * dt, Hmpo, Operators, param,
quenchtime, False, qt)
self.trotter_layer(c4 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
self.trotter_layer(c4 * dt, Hmpo, Operators, param,
quenchtime, False, qt)
self.trotter_layer(c3 * dt, Hmpo, Operators, param,
quenchtime, True, qt)
else:
raise Exception('Trotter order not implemented!')
self._trotter_finalize(qt, Hmpo, param, quenchtime)
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
[docs] def apply_op2(self, Op, ii, ll_1_term=False, perm=False):
"""
Apply an operator for two neighboring sites to the quantum system.
Operator can either be for pure states or density matrix (Liouville
space).
**Arguments**
Op : np 2d array
containing the operator to be applied, e.g. the propagator of
a Hamiltonian or Liouvillian.
ii : int
Apply term to the sites ii and (ii + 1).
ll_1_term : bool, optional
Boolean to indicate the term ii=L for periodic boundary
conditions.
Default to ``False``.
"""
if(self.SymmSec is None):
self._apply_op2_nosymm(Op, ii, ll_1_term, perm)
else:
self._apply_op2_symm(Op, ii, ll_1_term, perm)
return
[docs] def apply_op(self, Op, idx):
"""
Apply an operator or gate for any number of sites to the quantum
systems. Operatorcan either be for pure states or for density
matrices, but is defined in the Hilbert space for both of them
**Arguments**
Op : np 2d array
containing the operator to be applied, e.g. the propagator of
a Hamiltonian or Liouvillian.
idx : list of integers
Apply term to the sites specified here.
"""
if(self.SymmSec is None):
self._apply_op_nosymm(Op, idx)
else:
self._apply_op_symm(Op, idx)
return
[docs] def trotter_Ham_2site(self, Hmpo, Operators, param, quenchtime, even):
"""
Build the Hamiltonian for 2 neighboring site containing only local
and bond rules.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Defines the Hamiltonnian of the system.
Operators : instance of :py:class:`Ops.OperatorList`
Contains the operators used to build the Hamiltonian.
param : dictionary
Setup for the simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
even : bool
Boolean to switch between even (``True``) and odd (``False``)
layers in the Trotter decomposition.
"""
keys = ['exponential', 'product', 'FiniteFunction', 'MBString',
'InfiniteFunction']
for key in keys:
if(len(Hmpo.data[key]) > 0):
raise Exception(key + ' term not supported in Trotter '
+ 'decomposition!.')
ld = Operators['I'].shape[0]
ll = param['L']
loop = [0, 0]
if(even):
loop[0] = 1
loop[1] = (ll - 2) if(ll%2 == 0) else (ll - 1)
else:
loop[0] = 0
loop[1] = (ll - 2) if(ll%2 == 1) else (ll - 1)
nsite = len(Hmpo.data['site'])
nbond = len(Hmpo.data['bond'])
for jj in range(loop[0], loop[1], 2):
if(jj in self.lastprop):
yield self.lastprop[jj]
continue
Op2 = np.zeros((ld**2, ld**2), dtype=np.complex128)
for ii in range(nsite):
hp = Hmpo.data['site'].hparam(ii, param, quenchtime)
if(Hmpo.pbc):
hp *= 0.5
else:
hp[1:-1] = 0.5 * np.array(hp[1:-1])
Op = Operators[Hmpo.data['site']['Op'][ii]]
Op2 += np.kron(hp[jj] * Op, np.eye(ld))
Op2 += np.kron(np.eye(ld), hp[jj + 1] * Op)
# Loop over bond terms
for ii in range(nbond):
hp = Hmpo.data['bond'].hparam(ii, param, quenchtime)
Op = [Operators[Hmpo.data['bond']['Op'][ii][0]],
Operators[Hmpo.data['bond']['Op'][ii][1]]]
Op2 += np.kron(hp[jj] * Op[0], Op[1])
# Return operator
self.lastprop[jj] = (Op2, jj, False)
yield (Op2, jj, False)
if(Hmpo.pbc and even):
if(ll % 2 != 0):
raise Exception("Periodic boundary only with even System size!")
if(ll - 1 in self.lastprop):
yield self.lastprop[ll - 1]
else:
Op2 = np.zeros((ld**2, ld**2), dtype=np.complex128)
for ii in range(nsite):
hp = Hmpo.data['site'].hparam(ii, param, quenchtime)
hp *= 0.5
Op = Operators[Hmpo.data['site']['Op'][ii]]
Op2 += np.kron(hp[ll - 1] * Op, np.eye(ld))
Op2 += np.kron(np.eye(ld), hp[0] * Op)
# Loop over bond terms
for ii in range(nbond):
hp = Hmpo.data['bond'].hparam(ii, param, quenchtime)
Op = [Operators[Hmpo.data['bond']['Op'][ii][0]],
Operators[Hmpo.data['bond']['Op'][ii][1]]]
Op2 += np.kron(hp[ll - 1] * Op[0], Op[1])
# Return operator
self.lastprop[ll - 1] = (Op2, ll - 1, True)
yield (Op2, ll - 1, True)
[docs] def save_propagator_spectrum(self, Mat, param, quenchtime, tol=1e-13):
"""
Save the eigenvalues of the spectrum of the Hamiltonian or
Liouville operator.
**Arguments**
Mat : 2d matrix
Contains the Hamiltonian or Liouville operator with the
eigenvalues to be determined.
param : dictionary
Contains the simulation setup for one simulation.
quenchtime : tuple of integer and float
Contains the number of the corresponding quench and the
current time of the simulation.
tol : float, optional
Tolerance for the check of hermiticity.
Default to 1e-13.
**Details**
"""
global get_spectrum
if(not get_spectrum): return
if(Mat.shape[0] <= 2**12):
# Do not use ARPACK for "small" matrices
Mat = Mat.toarray()
herm = ((np.abs(Mat - Mat.transpose().conj())).sum() < tol)
if(isinstance(Mat, np.ndarray)):
# Dense matrix as numpy 2d array
if(herm):
vals, vecs = nla.eigh(Mat)
else:
vals, vecs = nla.eig(Mat)
else:
# Assume it is a sparse matrix
kk = 6
while(kk < 0.5 * Mat.shape[0]):
if(herm):
vals, vecs = spla.eigsh(Mat, which='SM')
else:
vals, vecs = spla.eigs(Mat, which='SM')
if(np.abs(vals[0] - vals[-1]) > 1e-12): break
kk += 6
# Construct a filename
base = param['Output_Directory'] + param['job_ID'] + param['unique_ID']
base += '_t%3.6f'%(quenchtime[1])
fvals = base + 'eigvals.npy'
fvecs = base + 'eigvecs.npy'
np.save(fvals, vals)
np.save(fvecs, vecs)
return
# ==============================================================================
# METHODS FOR PURE STATES
# ==============================================================================
[docs]class PureState(_QuantumSystem):
"""
Pure state vector for the evolution of closed systems.
**Arguments**
Psi : 1d numpy array
ld : int
local dimension in the Hilbert space
ll : int
number of sites in the system
SymmSec : None or instance of :py:class:`exactdiag.SymmetrySector`
contains the setup of the symmetry if present in the system.
**Variables**
dd : int
dimension of Psi
rhos : dict
Storing reduced density matrices for measurements.
lastprop : None, 2d numpy array, or dictionary.
Last propagator used for time independent quenches either as
numpy 2d array for the matrix exponential or as dictionary
for the Trotter decomposition.
bidx : dictionary
Used to store the indices of the subblocks with the same
symmetry numbers.
"""
def __init__(self, Psi, ld, ll, SymmSec):
self.Psi = Psi
self.ld = ld
self.ll = ll
self.SymmSec = SymmSec
self.dd = Psi.shape[0]
self.rhos = {}
self.lastprop = None
self.bidx = None
return
[docs] def save(self, flnm):
np.save(flnm, self.Psi)
# --------------------- DENSITY MATRICES -----------------------------------
[docs] def _traceout_nosymm(self, trce, keep):
"""
Trace over a set of sites specified via their indices for states
without a symmetry.
**Arguments**
trce : list of ints
Trace out over the indices contain in this list.
keep : list of ints
Keep the sites whose indices are listed here. Together trce
and keep should list all indices exactly one time.
"""
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
if(isinstance(self.ld, int)):
self.Psi = np.reshape(self.Psi, [self.ld] * self.ll)
else:
self.Psi = np.reshape(self.Psi, self.ld)
rho = np.tensordot(self.Psi, np.conj(self.Psi), (trce, trce))
if(isinstance(self.ld, int)):
self.Psi = np.reshape(self.Psi, self.ld**self.ll)
rho = np.reshape(rho, [self.ld**(len(keep))] * 2)
rho = DensityMatrix(rho, self.ld, self.ll - len(trce), None)
else:
self.Psi = np.reshape(self.Psi, np.prod(self.ld))
rho = np.reshape(rho, [np.prod(self.ld[keep])] * 2)
rho = DensityMatrix(rho, self.ld[keep], self.ll - len(trce), None)
return rho
[docs] def _traceout_symm2nosymm(self, trce, keep, *args, **kwargs):
"""
Trace over a set of sites specified via their indices for states
with symmetry returning a density matrix not considering the
symmetry.
**Arguments**
trce : list of ints
Trace out over the indices contain in this list.
keep : list of ints
Keep the sites whose indices are listed here. Together trce
and keep should list all indices exactly one time.
``*args``, ``**kwargs``
keeping syntax with corresponding method for density matrices.
"""
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
if(isinstance(self.ld, int)):
# Dimension of reduced density matrix
dim = self.ld**len(keep)
# Get indicees of subsystem kept
basek = np.zeros(self.ll, dtype=np.float64)
basek[keep] = self.ld**(linspace(len(keep) - 1, 0,
len(keep), dtype=np.float64))
idx1 = np.array(np.dot(self.SymmSec.sectorbase, basek),
dtype=np.int32)
# Get indices of subsystem traced out
shape = (self.ld**len(trce), 1)
if(shape[0] < 2**31):
# Construct indices on complete space
baset = np.zeros(self.ll, dtype=np.float64)
baset[trce] = self.ld**(linspace(len(trce) - 1, 0,
len(trce), dtype=np.float64))
idx2 = np.array(np.dot(self.SymmSec.sectorbase, baset),
dtype=np.int32)
else:
# Construct indices on subspace
TrcSec = self.SymmSec.get_subspace(trce)
idx2 = np.zeros((len(self.SymmSec)), dtype=np.int32)
for kk in range(len(self.SymmSec)):
idx2[kk] = TrcSec[tuple(self.SymmSec.sectorbase[kk, trce])]
shape = (len(TrcSec), 1)
del TrcSec
# Store sparse vectors in dictionary
Phi = {}
Phic = {}
if(shape[0] >= 2**31):
# Just in case - avoid breakdown of scipy
raise Exception("Complete Hilbert space exceeds 2^31")
for kk in range(dim):
idx = (idx1 == kk)
ptr = np.array([0, np.sum(idx)], np.int32)
Phi[kk] = sp.csc_matrix((self.Psi[idx], idx2[idx], ptr),
shape=shape)
Phic[kk] = sp.csc_matrix(Phi[kk].transpose().conj())
# Calculate reduced density matrix
rho = np.zeros((dim, dim), dtype=self.Psi.dtype)
for kk in range(dim - 1):
rho[kk, kk] = Phic[kk].dot(Phi[kk])[0, 0]
for jj in range(kk + 1, dim):
rho[kk, jj] = Phic[jj].dot(Phi[kk])[0, 0]
rho[jj, kk] = np.conj(rho[kk, jj])
rho[dim - 1, dim - 1] = 1.0 - np.trace(rho)
rho = DensityMatrix(rho, self.ld, self.ll - len(trce), None)
"""
if(False):#ld**ll / Psi.shape[0] < 10):
# Find contracting indices via booleans on complete space
# .......................................................
masks = {}
for kk in range(self.ld):
masks[kk] = np.where(idx1 == kk)[0]
#masks[kk] = idx1 == kk
rho = np.zeros((self.ld, self.ld), dtype=self.Psi.dtype)
for kk in range(self.ld - 1):
rho[kk, kk] = np.sum(np.conj(self.Psi[masks[kk]])
* self.Psi[masks[kk]])
veck = np.zeros((self.ld**(self.ll - 1)), dtype=np.bool)
veck[idx2[masks[kk]]] = True
for jj in range(kk + 1, self.ld):
vecj = np.zeros((self.ld**(self.ll - 1)), dtype=np.bool)
vecj[idx2[masks[jj]]] = True
vecj = np.logical_and(veck, vecj)
rho[kk, jj] = np.sum(np.conj(self.Psi[vecj[idx2]]) *
self.Psi[vecj[idx2]])
rho[jj, kk] = np.conj(rho[kk, jj])
rho[self.ld - 1, self.ld - 1] = 1.0 - np.trace(rho)
"""
else:
# Dimension of reduced density matrix
dim = np.prod(self.ld[keep])
# Get indices of subsystem kept
basek = np.zeros(self.ll, dtype=np.float64)
tmp = np.cumprod(self.ld[keep[::-1]])[::-1] / self.ld[keep[-1]]
basek[keep] = tmp
idx1 = np.array(np.dot(self.SymmSec.sectorbase, basek),
dtype=np.int32)
# Get indices of subsytem traced out
shape = (np.prod(self.ld[trce]), 1)
if(shape[0] < 2**31):
baset = np.zeros(self.ll, dtype=np.float64)
tmp = np.cumprod(self.ld[trce[::-1]])[::-1] / self.ld[trce[-1]]
baset[trce] = tmp
idx2 = np.array(np.dot(self.SymmSec.sectorbase, baset),
dtype=np.int32)
else:
# Construct indices on subspace
TrcSec = self.SymmSec.get_subspace(trce)
idx2 = np.zeros((len(self.SymmSec)), dtype=np.int32)
for kk in range(len(self.SymmSec)):
idx2[kk] = TrcSec[tuple(self.SymmSec.sectorbase[kk, trce])]
shape = (len(TrcSec), 1)
del TrcSec
# Store sparse vectors in dictionary
Phi = {}
Phic = {}
if(shape[0] >= 2**31):
# Just in case - avoid breakdown in scipy
raise Exception("Complete Hilbert space exceeds 2^31")
for kk in range(dim):
idx = (idx1 == kk)
ptr = np.array([0, np.sum(idx)], np.int32)
Phi[kk] = sp.csc_matrix((self.Psi[idx], idx2[idx], ptr),
shape=shape)
Phic[kk] = sp.csc_matrix(Phi[kk].transpose().conj())
# Calculate reduced density matrix
rho = np.zeros((dim, dim), dtype=self.Psi.dtype)
for kk in range(dim - 1):
rho[kk, kk] = Phic[kk].dot(Phi[kk])[0, 0]
for jj in range(kk + 1, dim):
rho[kk, jj] = Phic[jj].dot(Phi[kk])[0, 0]
rho[jj, kk] = np.conj(rho[kk, jj])
rho[dim - 1, dim - 1] = 1.0 - np.trace(rho)
rho = DensityMatrix(rho, self.ld[keep], self.ll - len(trce), None)
return rho
[docs] def _traceout_symm(self, trce, keep, SubSec=None):
"""
Trace over a set of sites specified via their indices for states
with symmetry returning a density matrix considering the symmetry
in the larger Hilbert space.
**Arguments**
trce : list of ints
Trace out over the indices contain in this list.
keep : list of ints
Keep the sites whose indices are listed here. Together trce
and keep should list all indices exactly one time.
SubSec : instance of :py:class:`exactdiag.SymmetrySector`, optional
If symmetry sector is already calculated, it can be passed
with this argument. Otherwise pass ``None``.
Default to ``None``.
"""
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
if(SubSec is None):
# Get dimension of reduced density matrix
SubSec = self.SymmSec.get_subspace(keep)
dim = len(SubSec)
if(isinstance(self.ld, int)):
# Get indices of subsystem kept
idx1 = np.zeros((len(self.SymmSec)), dtype=np.int32)
for kk in range(len(self.SymmSec)):
idx1[kk] = SubSec[self.SymmSec.sectorbase[kk, keep]]
# Get indices of subsystem traced out
shape = (self.ld**len(trce), 1)
if(shape[0] < 2**31):
# Construct indices on complete space
baset = np.zeros(self.ll, dtype=np.float64)
baset[trce] = self.ld**(linspace(len(trce) - 1, 0,
len(trce), dtype=np.float64))
idx2 = np.array(np.dot(self.SymmSec.sectorbase, baset),
dtype=np.int32)
else:
# Construct indices on subspace
TrcSec = self.SymmSec.get_subspace(trce)
idx2 = np.zeros((len(self.SymmSec)), dtype=np.int32)
for kk in range(len(self.SymmSec)):
idx2[kk] = TrcSec[tuple(self.SymmSec.sectorbase[kk, trce])]
shape = (len(TrcSec), 1)
del TrcSec
# Store sparse vectors in dictionary
Phi = {}
Phic = {}
if(shape[0] >= 2**31):
# Just in case - avoid breakdown in scipy
raise Exception("Complete Hilbert space exceeds 2^31")
for kk in range(dim):
idx = (idx1 == kk)
ptr = np.array([0, np.sum(idx)], np.int32)
Phi[kk] = sp.csc_matrix((self.Psi[idx], idx2[idx], ptr),
shape=shape)
Phic[kk] = sp.csc_matrix(Phi[kk].transpose().conj())
# Calculate reduced density matrix
rho = np.zeros((dim, dim), dtype=self.Psi.dtype)
for kk in range(dim - 1):
rho[kk, kk] = Phic[kk].dot(Phi[kk])[0, 0]
for jj in range(kk + 1, dim):
rho[kk, jj] = Phic[jj].dot(Phi[kk])[0, 0]
rho[jj, kk] = np.conj(rho[kk, jj])
rho[dim - 1, dim - 1] = 1.0 - np.trace(rho)
rho = DensityMatrix(rho, self.ld, self.ll - len(trce), SubSec)
else:
# Get indices of subsystem kept
idx1 = np.zeros((dim), dtype=np.int32)
for kk in range(len(self.SymmSec)):
idx1[kk] = SubSec[self.SymmSec.sectorbase[kk, keep]]
# Get indices of subsystem traced out
shape = (np.prod(self.ld[trce]), 1)
if(shape[0] < 2**31):
baset = np.zeros(self.ll, dtype=np.float64)
tmp = np.cumprod(self.ld[trce[::-1]])[::-1] / self.ld[trce[-1]]
baset[trce] = tmp
idx2 = np.array(np.dot(self.SymmSec.sectorbase, baset),
dtype=np.int32)
else:
SubTrc = self.SymmSec.get_subspace(trce)
idx2 = np.zeros((len(self.SymmSec)), dtype=np.int32)
for kk in range(len(self.SymmSec)):
idx2[kk] = SubTrc[tuple(self.SymmSec.sectorbase[kk, trce])]
shape = (len(SubTrc), 1)
del SubTrc
# Store sparse vectors in dictionary
Phi = {}
Phic = {}
if(shape[0] >= 2**31):
# Just in case - avoid breakdown in scipy
raise Exception("Complete Hilbert space exceeds 2^31")
for kk in range(dim):
idx = (idx1 == kk)
ptr = np.array([0, np.sum(idx)], np.int32)
Phi[kk] = sp.csc_matrix((self.Psi[idx], idx2[idx], ptr),
shape=shape)
Phic[kk] = sp.csc_matrix(Phi[kk].transpose().conj())
# Calculate reduced density matrix
rho = np.zeros((dim, dim), dtype=self.Psi.dtype)
for kk in range(dim - 1):
rho[kk, kk] = Phic[kk].dot(Phi[kk])[0, 0]
for jj in range(kk + 1, dim):
rho[kk, jj] = Phic[jj].dot(Phi[kk])[0, 0]
rho[jj, kk] = np.conj(rho[kk, jj])
rho[dim - 1, dim - 1] = 1.0 - np.trace(rho)
rho = DensityMatrix(rho, self.ld[keep], self.ll - len(trce), SubSec)
return rho
# ------------------------ ENTROPY -----------------------------------------
[docs] def entropy(self):
"""
Calculate the von-Neumann entropy of the system (always zero for the
pure state).
"""
return 0.0
# ---------------- Calculating bond entropies ------------------------------
[docs] def get_bond_entropies(self):
"""
Calculate the bond entropies for every possible splitting along the
chain. Iterator returns ll + 1 numbers.
"""
be = np.zeros((self.ll + 1))
# Left bipartition
for elem in self._bond_entropy_to_left():
be[elem[0]] = elem[1]
# Right bipartion
for elem in self._bond_entropy_to_right():
be[elem[0]] = elem[1]
for ii in range(self.ll + 1):
yield be[ii]
[docs] def _bond_entropy_to_left(self):
"""
Start from ii = ll / 2 and calculate bond entropies going to the left.
"""
ii = self.ll // 2 - 1
nn = ii + 1
keep = list(range(ii + 1))
rho = self.reduced(keep)
for kk in range(ii + 1, 0, -1):
be = rho.entropy()
if(nn > 1): rho = rho.traceout_last()
nn -= 1
yield (kk, be)
[docs] def _bond_entropy_to_right(self):
"""
Start from ii = ll / 2 + 1 and calculate bond entropies going to the
left.
"""
ii = self.ll // 2 + 1
nn = self.ll - ii
trce = list(range(ii))
rho = self.traceout(trce)
for kk in range(ii, self.ll):
be = rho.entropy()
if(nn > 1): rho = rho.traceout([0])
nn -= 1
yield (kk, be)
[docs] def _get_bond_entropy_nosymm(self, ii):
"""
Calculate the bond entropy for a state of splitting (0, ..., ii),
(ii + 1, ..., L - 1) in a system where no symmetries are used. For
description of other arguments look into ``get_bond_entropies``.
**Arguments**
ii : int
Specifies the index of the splitting, where ii is the last
site included in the left bipartition.
**Notes**
Could use scipy.linalg.blas.zherk or scipy.linalg.blas.dsyrk
"""
left = ii + 1
right = self.ll - left
if(left > right):
# Build reduced density matrix rho^dagger of right bipartition
# [ii+1, ..., ll]
cut = self.ll - ii - 1
tmp = np.reshape(self.Psi, (self.ld**(ii + 1),
self.ld**(self.ll - ii - 1)))
rho = np.dot(np.transpose(np.conj(tmp)), tmp)
else:
# Build reduced density matrix rho of left bipartition
# [1, ..., ii]
cut = ii + 1
tmp = np.reshape(self.Psi, (self.ld**(ii + 1),
self.ld**(self.ll - ii - 1)))
rho = np.dot(tmp, np.transpose(np.conj(tmp)))
return DensityMatrix(rho, self.ld, cut, None)
[docs] def _get_bond_entropy_symm(self, ii):
"""
Calculate the bond entropy for a state of splitting (0, ..., ii),
(ii + 1, ..., L - 1) in a system with symmetries. For description of
arguments look into ``get_bond_entropies``.
**Arguments**
ii : int
Specifies the index of the splitting, where ii is the last
site included in the left bipartition.
"""
left = ii + 1
right = self.ll - left
if(left > right):
# Build reduced density matrix of right bipartition
# [ii+1, ..., ll-1]
perm = []
for kk in range(ii + 1, self.ll):
perm.append(kk)
cut = len(perm)
for kk in range(ii + 1):
perm.append(kk)
else:
# Build reduced density matrix of left bipartition [0, ..., ii]
perm = list(range(self.ll))
cut = ii + 1
newbase = self.SymmSec.permute(perm)
base2 = self.ld**(linspace(self.ll - cut - 1, 0, self.ll - cut,
dtype=np.int32))
base2p = self.ld**(linspace(cut - 1, 0, cut, dtype=np.int32))
Phi = sp.lil_matrix((self.ld**cut, self.ld**(self.ll - cut)),
dtype=self.Psi.dtype)
for kk in range(self.Psi.shape[0]):
idx1 = np.sum(base2p * newbase[kk, :cut])
idx2 = np.sum(base2 * newbase[kk, cut:])
Phi[idx1, idx2] = self.Psi[kk]
Phic = sp.csr_matrix(Phi.transpose().conj())
Phi = sp.csr_matrix(Phi)
rho = Phi.dot(Phic).toarray()
if(np.abs(1 - np.trace(rho)) > 1e-11):
raise ValueError('Normalization violated with '
+ '%3.6e!'%(np.abs(1 - np.trace(rho))))
return DensityMatrix(rho, self.ld, cut, None)
[docs] def _get_bond_entropy_symm_sub(self, ii):
"""
"""
left = ii + 1
right = self.ll - left
if(left > right):
# Build reduced density matrix of right bipartition
# [ii+1, ..., ll-1]
perm = []
for kk in range(ii + 1, self.ll):
perm.append(kk)
cut = len(perm)
for kk in range(ii + 1):
perm.append(kk)
else:
# Build reduced density matrix of left bipartition [0, ..., ii]
perm = list(range(self.ll))
cut = ii + 1
newbase = self.SymmSec.permute(perm)
base = self.ld**(linspace(self.ll - cut - 1, 0, self.ll - cut,
dtype=np.int32))
SubSec = self.SymmSec.get_subspace(perm[:cut])
Phi = sp.lil_matrix((len(SubSec), self.ld**(self.ll - cut)),
dtype=self.Psi.dtype)
for kk in range(self.Psi.shape[0]):
idx1 = SubSec[self.SymmSec.sectorbase[kk, perm[:cut]]]
idx2 = np.sum(base * newbase[kk, cut:])
Phi[idx1, idx2] = self.Psi[kk]
Phic = sp.csr_matrix(Phi.transpose().conj())
Phi = sp.csr_matrix(Phi)
rho = Phi.dot(Phic).toarray()
if(np.abs(1 - np.trace(rho)) > 1e-11):
raise ValueError('Normalization violated with ' +
'%3.6e!'%(np.abs(1 - np.trace(rho))))
return DensityMatrix(rho, self.ld, cut, SubSec)
# ----------------------------- Dot products -------------------------------
[docs] def norm(self):
"""
Returns the norm of the wave function.
"""
return np.dot(np.conj(self.Psi), self.Psi)
[docs] def dot(self, *args):
"""
General dot product split into different cases.
"""
raise NotImplementedError('General dot for PureState')
[docs] def dot_vec(self, Vec):
"""
Return overlap of wavefunction with other state vector
:math:`\langle self | Vec \\rangle`.
**Arguments**
Vec : np 1d array or instance of :py:class:`exactdiag.PureState`
Other vector in the overlap with the PureState.
"""
if(isinstance(Vec, PureState)):
return np.dot(np.conj(self.Psi), Vec.Psi)
else:
return np.dot(np.conj(self.Psi), Vec)
[docs] def dot_hmpo_vec(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a Hamiltonian represented as MPO / MPO rule set with a vector.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
defines the Hamiltonian through a rule set.
Operators : instance or :py:class:`Ops.OperatorList`
containing all operators needed defined on a local Hilbert space.
param : dictionary
setup for the simulation
quenchtime : tuple or None
If tuple containing index of quench and present time to obtain
coupling.
maxdim : int
maximal dimension of the Hilbert space
"""
if(self.SymmSec is None):
return self._dot_hmpo_vec_nosymm(Hmpo, Operators, param,
quenchtime, maxdim)
return self._dot_hmpo_vec_symm(Hmpo, Operators, param,
quenchtime, maxdim)
[docs] def _dot_hmpo_vec_nosymm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a Hamiltonian represented as MPO / MPO rule set with a vector.
For arguments see :py:func:`exactdiag.PureState.dot_hmpo_vec`. Intended
for internal use.
"""
# Get shortcuts to system size and local dimension of Hilbert space
ll = param['L']
ld = Operators['I'].shape[0]
dim = ld**ll
shape = [ld] * ll
self.Psi = np.reshape(self.Psi, shape)
if(self.Psi.dtype == np.complex128):
Phi = np.zeros((shape), dtype=np.complex128)
else:
Phi = np.zeros((shape), dtype=np.float64)
# Loop over local terms
for ii in range(len(Hmpo.data['site'])):
hp = Hmpo.data['site'].hparam(ii, param, quenchtime)
Op = Operators[Hmpo.data['site']['Op'][ii]]
for jj in range(ll):
perm = list(range(1, jj + 1)) + [0] + \
list(range(jj + 1, ll))
Phi += np.transpose(np.tensordot(hp[jj] * Op, self.Psi,
([1], [jj])), perm)
# Loop over bond terms
for ii in range(len(Hmpo.data['bond'])):
hp = Hmpo.data['bond'].hparam(ii, param, quenchtime)
Op = [Operators[Hmpo.data['bond']['Op'][ii][0]],
Operators[Hmpo.data['bond']['Op'][ii][1]]]
for jj in range(ll - 1):
perm = list(range(2, jj + 2)) + [1, 0] + \
list(range(jj + 2, ll))
tmp = np.tensordot(hp[jj] * Op[0], self.Psi,
([1], [jj]))
Phi += np.transpose(np.tensordot(Op[1], tmp,
([1], [jj + 1])), perm)
if(Hmpo.pbc):
perm = [0] + list(range(2, ll)) + [1]
tmp = np.tensordot(hp[ll - 1] * Op[0], self.Psi,
([1], [ll - 1]))
Phi += np.transpose(np.tensordot(Op[1], tmp, ([1], [1])),
perm)
# Loop over exponential terms
for ii in range(len(Hmpo.data['exponential'])):
hp = Hmpo.data['exponential'].hparam(ii, param, quenchtime)
Op = [Operators[Hmpo.data['exponential']['Op'][ii][0]],
Operators[Hmpo.data['exponential']['Op'][ii][1]],
Operators[Hmpo.data['exponential']['Op'][ii][2]]]
dp = Hmpo.data['exponential']['dp'][ii]
# Check if phase is not identity
phase = (Hmpo.data['exponential']['Op'][ii][2] != 'I')
if(phase):
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
perm = list(range(dist + 1, jj + dist + 1)) + \
list(range(dist + 1))[::-1] + \
list(range(jj + dist + 1, ll))
tmp = np.tensordot(hp[jj] * Op[0], self.Psi,
([1], [jj]))
for nn in range(dist - 1):
tmp = np.tensordot(dp * Op[2], tmp,
([1], [jj + 1 + nn]))
Phi += np.transpose(np.tensordot(Op[1], tmp,
([1], [jj + dist])),
perm)
else:
# No phase - identities
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
perm = list(range(2, jj + 2)) + [1] \
+ list(range(jj + 2, jj + dist + 1)) + [0] \
+ list(range(jj + dist + 1, ll))
tmp = np.tensordot(hp[jj] * Op[0], self.Psi,
([1], [jj]))
tmp = np.tensordot(dp**(dist - 1) * Op[1], tmp,
([1], [jj + dist]))
Phi += np.transpose(tmp, perm)
if(len(Hmpo.data['product']) > 0):
raise NotImplementedError('Product terms not implemented for ED.')
if(len(Hmpo.data['FiniteFunction']) > 0):
raise NotImplementedError('FiniteFunction terms not implemented' +
' for ED.')
# Loop over InfiniteFunction terms
for ii in range(len(Hmpo.data['InfiniteFunction'])):
hp = Hmpo.data['InfiniteFunction'].hparam(ii, param, quenchtime)
Opl = Operators[Hmpo.data['InfiniteFunction']['Op'][ii][0]]
Opr = Operators[Hmpo.data['InfiniteFunction']['Op'][ii][1]]
Opp = Operators[Hmpo.data['InfiniteFunction']['Op'][ii][2]]
# Check if phase is not identity
phase = Hmpo.data['InfiniteFunction']['Op'][ii][2] != 'I'
if(phase):
for jj in range(ll - 1):
for hh in range(jj + 1, ll):
dist = hh - jj
sc = Hmpo.data['InfiniteFunction'].yvals[ii][dist - 1]
sc *= hp[jj]
perm = list(range(dist + 1, jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(jj + dist + 1, ll))
tmp = np.tensordot(sc * Opl, self.Psi, ([1], [jj]))
for nn in range(dist - 1):
tmp = np.tensordot(Opp, tmp, ([1], [jj + 1 + nn]))
Phi += np.transpose(np.tensordot(Opr, tmp,
([1], [jj + dist])),
perm)
else:
# No phase - identities
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
sc = Hmpo.data['InfiniteFunction'].yvals[ii][dist - 1]
sc *= hp[jj]
perm = list(range(2, jj + 2)) + [1] \
+ list(range(jj + 2, jj + dist + 1)) + [0] \
+ list(range(jj + dist + 1, ll))
tmp = np.tensordot(sc * Opl, self.Psi, ([1], [jj]))
Phi += np.transpose(np.tensordot(Opr, tmp,
([1], [jj + dist])),
perm)
# Loop over MBString terms
for ii in range(len(Hmpo.data['MBString'])):
hp = Hmpo.data['MBString'].hparam(ii, param, quenchtime)
Op = []
for opstr in Hmpo.data['MBString']['Op'][ii]:
Op.append(Operators[opstr])
nk = len(Hmpo.data['MBString']['Op'][ii])
for jj in range(ll - nk + 1):
# Last site with operator
kk = jj + nk - 1
# Adapt from exponential rule
dist = kk - jj
perm = list(range(dist + 1, jj + dist + 1)) + \
list(range(dist + 1))[::-1] + \
list(range(jj + dist + 1, ll))
tmp = np.tensordot(hp[jj] * Op[0], self.Psi,
([1], [jj]))
for kk in range(1, nk - 1):
tmp = np.tensordot(Op[kk], tmp, ([1], [jj + kk]))
Phi += np.transpose(np.tensordot(Op[-1], tmp,
([1], [jj + nk - 1])),
perm)
#if(Hmpo.data['lind1']):
# raise NotImplementedError('lind1 terms not implemented for ED.')
self.Psi = np.reshape(self.Psi, (dim))
return np.reshape(Phi, [dim])
[docs] def _dot_hmpo_vec_symm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a Hamiltonian represented as MPO / MPO rule set with a vector.
The vector is defined on the symmetry sector specified. For arguments
see :py:func:`exactdiag.PureState.dot_hmpo_vec`. Intended for internal
use.
"""
spdim = len(self.SymmSec)
dtype = self.Psi.dtype
Phi = np.zeros((spdim), dtype=self.Psi.dtype)
ld = Operators['I'].shape[0]
ll = param['L']
keylist = ['site', 'bond', 'exponential', 'MBString', 'InfiniteFunction']
for key in keylist:
if(len(Hmpo.data[key]) == 0): continue
for k1, k2, elem in Hmpo.data[key].build_symm(param, self.SymmSec,
quenchtime, ll,
Hmpo.pbc):
Phi[k1] += elem * self.Psi[k2]
if(len(Hmpo.data['product']) > 0):
raise NotImplementedError('Product terms not implemented for ED.')
if(len(Hmpo.data['FiniteFunction']) > 0):
raise NotImplementedError('FiniteFunction terms not implemented' +
' for ED.')
#if(Hmpo.data['lind1']):
# raise NotImplementedError('lind1 terms not implemented for ED.')
return Phi
[docs] def dot_effhmpo_vec(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply the effective Hamiltonian represented as MPO / MPO rule
set with a vector.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
defines the effective Hamiltonian through a rule set. The
open system part contributes -0.5j Ldagger L (where the imaginary
unit is cancelled out when mulitplying with (-1j) in the time
evolution.
Operators : instance or :py:class:`Ops.OperatorList`
containing all operators needed defined on a local Hilbert space.
param : dictionary
setup for the simulation
quenchtime : tuple or None
If tuple containing index of quench and present time to obtain
coupling.
maxdim : int
maximal dimension of the Hilbert space
"""
if(self.SymmSec is None):
return self._dot_effhmpo_vec_nosymm(Hmpo, Operators, param,
quenchtime, maxdim)
return self._dot_effhmpo_vec_symm(Hmpo, Operators, param,
quenchtime, maxdim)
[docs] def _dot_effhmpo_vec_nosymm(self, Hmpo, Operators, param, quenchtime,
maxdim):
"""
Multiply the effective Hamiltonian represented as MPO / MPO rule set
with a vector. For arguments see
:py:func:`exactdiag.PureState.dot_effhmpo_vec`. Intended for internal
use.
"""
# Get all terms from the Hamiltonian
Phi = self._dot_hmpo_vec_nosymm(self, Hmpo, Operators, param,
quenchtime, maxdim)
# Shortcuts
ll = param['L']
ld = Operators['I'].shape[0]
dim = ld**ll
shape = [ld] * ll
self.Psi = np.reshape(self.Psi, shape)
Phi = np.reshape(Phi, shape)
if(Phi.dtype != np.complex128): Phi = Phi.astype(np.complex128)
# Loop over Lindblad terms
for ii in range(len(Hmpo.data['lind1'])):
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
Op = Operators[Hmpo.data['lind1']['Op'][ii]]
Op = -0.5j * Op.conj().transpose().dot(Op)
for jj in range(ll):
perm = list(range(1, jj + 1)) + [0] \
+ list(range(jj + 1, ll))
Phi += np.transpose(np.tensordot(hp[jj] * Op, self.Psi,
([1], [jj])), perm)
self.Psi = np.reshape(self.Psi, [dim])
return np.reshape(Phi, [dim])
[docs] def _dot_effhmpo_vec_symm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply the effective Hamiltonian represented as MPO / MPO rule set
with a vector. The vector is defined on the symmetry sector specified.
For arguments see :py:func:`exactdiag.PureState.dot_effhmpo_vec`.
Intended for internal use.
"""
# Get the vector containing all Hamiltonian terms
Phi = self._dot_hmpo_vec_symm(Hmpo, Operators, param, quenchtime,
maxdim)
# Shortcuts
spdim = len(self.SymmSec)
ld = Operators['I'].shape[0]
ll = param['L']
if(Phi.dtype != np.complex128): Phi = Phi.astype(np.complex128)
keylist = ['lind1']
for key in keylist:
if(len(Hmpo.data[key]) == 0): continue
for k1, k2, elem in Hmpo.data[key].build_symm(param, self.SymmSec,
quenchtime, ll,
Hmpo.pbc):
Phi[k1] += elem * self.Psi[k2]
return Phi
[docs] def meas_mpo(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Measure an MPO defined as rule set, e.g. the Hamiltonian.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
defines the Hamiltonian through a rule set.
Operators : instance or :py:class:`Ops.OperatorList`
containing all operators needed defined on a local Hilbert space.
param : dictionary
setup for the simulation
quenchtime : tuple or None
If tuple containing index of quench and present time to obtain
coupling.
maxdim : int
maximal dimension of the Hilbert space
"""
Phi = self.dot_hmpo_vec(Hmpo, Operators, param, quenchtime,
maxdim)
en = self.dot_vec(Phi)
del Phi
#Ham = build_hamiltonian(Hmpo, Operators, param, SymmSec, sparse,
# maxdim, quenchtime=(ii, time))
#en2 = np.real(np.dot(np.conj(Psi), Ham.dot(Psi)))
return en
# ----------------------------- Propagation --------------------------------
[docs] def propagate_expm(self, dt, Hmpo, param, sparse, maxdim, quenchtime, qt):
"""
Propagate a quantum state where the propagator is built from a sparse
or dense Hamiltonian.
**Arguments**
dt : float
size of the time step
Hmpo : instance of :py:class:`MPO.MPO`
defines the Hamiltonian through a rule set.
param : dictionary
setup for the simulation
sparse : bool
Flag on the use of sparse matrices. If ``True``, sparse matrices
will be used.
maxdim : int
maximal dimension of the Hilbert space
quenchtime : tuple
Tuple containing index of quench and present time to obtain
coupling.
qt : bool
flag for quantum trajectories.
"""
if(qt):
self.propagate_expm_qt(dt, Hmpo, param, sparse, maxdim, quenchtime)
return
if(self.lastprop is None):
expm = spla.expm if(sparse) else sla.expm
self.lastprop = Hmpo.build_hamiltonian(param, self.SymmSec,
sparse, maxdim,
quenchtime=quenchtime)
self.lastprop = expm(-1j * dt * self.lastprop)
self.Psi = self.lastprop.dot(self.Psi)
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
[docs] def propagate_expm_qt(self, dt, Hmpo, param, sparse, maxdim, quenchtime):
"""
Propagate a quantum state where the propagator is built from the
effective Hamiltonian.
**Arguments**
dt : float
size of the time step
Hmpo : instance of :py:class:`MPO.MPO`
defines the Hamiltonian through a rule set.
param : dictionary
setup for the simulation
sparse : bool
Flag on the use of sparse matrices. If ``True``, sparse matrices
will be used.
maxdim : int
maximal dimension of the Hilbert space
quenchtime : tuple or None
If tuple containing index of quench and present time to obtain
coupling.
"""
try:
tmp = self.qtnorm
except AttributeError:
self.qtnorm = self.norm()
self.qtr = np.random.rand(1)[0]
if(self.lastprop is None):
expm = spla.expm if(sparse) else sla.expm
self.lastprop = Hmpo.build_effhamiltonian(param, self.SymmSec,
sparse, maxdim,
quenchtime)
self.lastprop = expm(-1j * dt * self.lastprop)
self.Psi = self.lastprop.dot(self.Psi)
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
self.quantum_jump(Hmpo, param, quenchtime)
return
[docs] def quantum_jump(self, Hmpo, param, quenchtime):
"""
Apply a Lindblad operator for the quantum trajectory evolution. This
function works for any time evolution method with an (effective)
Hamiltonian (which are matrix exponential, Trotter, Krylov).
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
contains the MPO and Lindblad operators.
param : dict
contains the simulation setup.
quenchtime : tuple of (int, float)
information about the number of the quench and the present time to
gain the coupling constants for the Lindblad operators.
"""
norm = self.norm()
self.Psi /= np.sqrt(norm)
self.qtnorm *= norm
if(self.qtnorm < self.qtr):
if(self.SymmSec is None):
# Apply quantum jump (No symmetry)
self._quantum_jump_nosymm(Hmpo, param, quenchtime)
else:
# Apply quantum jump( Symmetry present)
self._quantum_jump_symm(Hmpo, param, quenchtime)
# Reset probabilities
self.qtnorm = 1.0
self.qtr = np.random.rand(1)[0]
return
[docs] def _quantum_jump_nosymm(self, Hmpo, param, quenchtime):
"""
Apply a quantum jump for the quantum trajectory evolution in a system
without any symmetry present. For arguments look into documentation
of :py:func:`exactdiag.PureState.quantum_jump`. Intended for internal
use.
"""
# Shortcuts
ll = param['L']
ld = Hmpo.Operators['I'].shape[0]
# Reshape wave function into tensor
self.Psi = np.reshape(self.Psi, [ld] * ll)
# Setup for probability
nlind1 = len(Hmpo.data['lind1']['hp'])
weight = np.zeros((ll * nlind1))
mapping = {}
kk = 0
for ii in range(nlind1):
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
Op = Hmpo.Operators[Hmpo.data['lind1']['Op'][ii]]
for jj in range(ll):
Phi = np.tensordot(hp[jj] * Op,
self.Psi, ([1], [jj])).reshape([ld**ll])
dot = np.abs(np.sum(Phi.conj() * Phi))
weight[kk] = dot
mapping[kk] = (ii, jj)
kk += 1
# Find which jump is to be applied
weight = np.cumsum(weight)
weight = weight / weight[-1]
rand = np.random.rand(1)[0]
# Default if rand = 1.0
idx = weight.shape[0]
for kk in range(weight.shape[0]):
if(rand < weight[kk]):
idx = kk
break
ii, jj = mapping[idx]
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
Op = Hmpo.Operators[Hmpo.data['lind1']['Op'][ii]]
# Apply quantum jump
perm = list(range(1, jj + 1)) + [0] + list(range(jj + 1, ll))
self.Psi = np.tensordot(hp[jj] * Op,
self.Psi, ([1], [jj])).transpose(perm)
# Reshape wave function into vector
self.Psi = np.reshape(self.Psi, [ld**ll])
# Renorm
self.Psi = self.Psi / np.sqrt(self.norm())
return
[docs] def _quantum_jump_symm(self, Hmpo, param, quenchtime):
"""
Apply a quantum jump for the quantum trajectory evolution in a system
with one or multiple symmetries present. For arguments look into
documentation of :py:func:`exactdiag.PureState.quantum_jump`. Intended
for internal use.
"""
# Shortcuts
ll = param['L']
ld = Hmpo.Operators['I'].shape[0]
# Array for wave function L | psi >
Phi = np.zeros((self.Psi.shape[0]), dtype=self.Psi.dtype)
# Setup for probability
nlind1 = len(Hmpo.data['lind1']['hp'])
weight = np.zeros((ll * nlind1))
mapping = {}
kk = 0
for ii in range(nlind1):
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
Op = Hmpo.Operators[Hmpo.data['lind1']['Op'][ii]]
colidx = {}
for jj in range(ld):
colidx[jj] = []
for kk in range(ld):
if(Op[jj, kk] != 0.0): colidx[jj].append(kk)
for jj in range(ll):
Phi[:] = 0.0
for k1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(k1)
idx2 = deepcopy(idx1)
for kk in colidx[idx1[jj]]:
idx2[jj] = kk
try:
k2 = self.SymmSec[tuple(idx2)]
except KeyError:
continue
Phi[k1] += hp[jj] * Op[idx1[jj], kk] * \
self.Psi[k2]
dot = np.real(np.sum(Phi.conj() * Phi))
weight[kk] = dot
mapping[kk] = (ii, jj)
kk += 1
# Find which jump is to be applied
weight = np.cumsum(weight)
weight = weight / weight[-1]
rand = np.random.rand(1)[0]
# Default if rand = 1.0
idx = weight.shape[0]
for kk in range(weight.shape[0]):
if(rand < weight[kk]):
idx = kk
break
ii, jj = mapping[idx]
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
Op = Hmpo.Operators[Hmpo.data['lind1']['Op'][ii]]
colidx = {}
for k1 in range(ld):
colidx[k1] = []
for k2 in range(ld):
if(Op[k1, k2] != 0.0): colidx[k1].append(k2)
Phi[:] = 0.0
for k1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(k1)
idx2 = deepcopy(idx1)
for kk in colidx[idx1[jj]]:
idx2[jj] = kk
try:
k2 = self.SymmSec[tuple(idx2)]
except KeyError:
continue
Phi[k1] += hp[jj] * Op[idx1[jj], kk] * \
self.Psi[k2]
# Copy and renorm
self.Psi = Phi
self.Psi = self.Psi / np.sqrt(self.norm())
return
[docs] def propagate_krylov(self, dt, Hmpo, Operators, param, sparse, maxdim,
quenchtime, qt):
"""
Propagate a quantum state taking the exponential in the Krylov subspace.
**Arguments**
dt : float
time step
Hmpo : MPO
Rule set describing the Hamiltonian.
Operators : instance of :py:class:`Ops.OperatorList`
contains all the operators to setup the Hamiltonian
param : dict
dictionary with simulation setup.
sparse : bool
Flag if sparse matrices should be used. Flag only used for
small matrices when complete Hamiltonian is built.
maxdim : int
maximal dimension of the Hilbert space.
quenchtime : tuple or None
contains number of quench and present time.
qt : bool
flag for quantum trajectories.
**Details**
The following modes are available for the Krylov time evolution,
which is set in the convergence criteria for the quenches, i.e.,
:py:class:`convergence.KrylovConvParam` with the variable
``edlibmode``:
* 0 : Built Hamiltonian as matrix and use built-in scipy method.
* 1 : Built Hamiltonian as matrix and use EDLib Krylov.
* 2 : Use dot product of MPO and state vector and EDLib Krylov.
* 3 : As 2, but save Krylov vectors to hard disk to reduce RAM memory
usage.
* -1 : default. For :math`dim < 2^{10}` take edlibmode=0, then
for :math:`dim < 2^{15}` use edlibmode=2, otherwise edlibmode=3
Since the Hamiltonian is hermitian, Krylov uses the Lanczos
implementation of the algorithm.
"""
if(qt):
self.propagate_krylov_qt(dt, Hmpo, Operators, param, sparse, maxdim,
quenchtime)
return
ConvParam = param['Quenches'].data[quenchtime[0]]
ConvParam = ConvParam['ConvergenceParameters'].data[0]
mode = ConvParam['edlibmode']
if(mode < 0):
if(self.Psi.shape[0] < 2**10):
# Complete H, built-in scipy
mode = 0
elif(self.Psi.shape[0] < 2**10):
# Complete H, EDLib Krylov (not used as default)
mode = 1
elif(self.Psi.shape[0] < 2**15):
# MPO, EDLib Krylov
mode = 2
else:
# MPO, EDLib Krylov, write Krylov vectors to hard drive
mode = 3
if(mode == 0):
# Complete H, built-in scipy
# --------------------------
if(self.lastprop is None):
self.lastprop = Hmpo.build_hamiltonian(param, self.SymmSec,
True, maxdim,
quenchtime=quenchtime)
self.Psi = spla.expm_multiply(-1j * dt * self.lastprop, self.Psi)
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
nvec = ConvParam['max_num_lanczos_iter']
tol = ConvParam['lanczos_tol']
numzero = 1e-16
alpha = np.zeros((nvec), dtype=np.float64)
beta = np.zeros((nvec + 1), dtype=np.float64)
if((mode == 1) or (mode == 2)):
# EDLib Krylov (Complete H, MPO dot product)
# ------------
base = np.zeros((self.Psi.shape[0], nvec), dtype=np.complex128)
beta[0] = np.sqrt(np.real(np.sum(np.conj(self.Psi) * self.Psi)))
if((mode == 1) and (self.lastprop is None)):
self.lastprop = Hmpo.build_hamiltonian(param, self.SymmSec,
sparse, maxdim,
quenchtime=quenchtime)
nn = nvec
for jj in range(nvec):
self.Psi = self.Psi / beta[jj]
base[:, jj] = self.Psi
if(mode == 1):
HPsi = self.lastprop.dot(base[:, jj])
else:
HPsi = self.dot_hmpo_vec(Hmpo, Operators, param,
quenchtime,
maxdim).astype(np.complex128)
alpha[jj] = np.real(np.sum(np.conj(self.Psi) * HPsi))
HPsi = HPsi - alpha[jj] * self.Psi
# Reorthogonalize
for ii in range(jj):
HPsi -= np.sum(np.conj(base[:, ii]) * HPsi) * base[:, ii]
beta[jj + 1] = np.sqrt(np.real(np.sum(np.conj(HPsi) * HPsi)))
prop = expm_symm_tridiag(alpha[:jj+1], beta[1:jj+1], -1j * dt)
est = 2.0 * np.abs(beta[jj]) * np.abs(prop[jj, 0])
if(est < tol):
nn = jj
del HPsi
break
if(abs(beta[jj + 1]) < 10 * numzero):
nn = jj
del HPsi
break
self.Psi = HPsi
self.Psi = prop[0, 0] * base[:, 0]
for jj in range(1, nn):
self.Psi += prop[jj, 0] * base[:, jj]
else:
# MPO, EDLib Krylov, write Krylov vectors to hard drive
# -----------------------------------------------------
# Save vectors to hard disk to save memory
fbase = param['Write_Directory'] + param['job_ID'] + \
param['unique_ID'] + 'Krylovbase_'
beta[0] = np.sqrt(np.real(np.sum(np.conj(self.Psi) * self.Psi)))
if((beta[0] - 1.0) > 1e-12):
print("Normalization:", beta[0])
raise Exception("Normalization violated.")
nn = nvec
for jj in range(nvec):
self.Psi = self.Psi / beta[jj]
np.save(fbase + str(jj) + '.npy', self.Psi)
HPsi = self.dot_hmpo_vec(Hmpo, Operators, param, quenchtime,
maxdim).astype(np.complex128)
alpha[jj] = np.real(np.sum(np.conj(self.Psi) * HPsi))
HPsi = HPsi - alpha[jj] * self.Psi
# Reorthogonalize
for ii in range(jj):
self.Psi = np.load(fbase + str(ii) + '.npy')
HPsi -= np.sum(np.conj(self.Psi) * HPsi) * self.Psi
beta[jj + 1] = np.sqrt(np.real(np.sum(np.conj(HPsi) * HPsi)))
prop = expm_symm_tridiag(alpha[:jj+1], beta[1:jj+1], -1j * dt)
est = 2.0 * np.abs(beta[jj]) * np.abs(prop[jj, 0])
if(est < tol):
nn = jj
del HPsi
break
if(abs(beta[jj + 1]) < 10 * numzero):
nn = jj
del HPsi
break
self.Psi = HPsi
self.Psi = np.load(fbase + str(0) + '.npy')
self.Psi = prop[0, 0] * self.Psi
for jj in range(1, nn):
flnm = fbase + str(jj) + '.npy'
tmp = np.load(flnm)
remove_file(flnm)
self.Psi += prop[jj, 0] * tmp
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
[docs] def propagate_krylov_qt(self, dt, Hmpo, Operators, param, sparse, maxdim,
quenchtime):
"""
Propagate a quantum state taking the exponential in the Krylov subspace.
Although this is on pure states, we need the Krylov-Arnoldi version to
account for non-hermiticity of the expression.
**Arguments**
dt : float
time step
Hmpo : MPO
Rule set describing the Hamiltonian.
Operators : instance of :py:class:`Ops.OperatorList`
contains all the operators to setup the Hamiltonian
param : dict
dictionary with simulation setup.
sparse : bool
Flag if sparse matrices should be used. Flag only used for
small matrices when complete Hamiltonian is built.
maxdim : int
maximal dimension of the Hilbert space.
quenchtime : tuple or None
contains number of quench and present time.
**Details**
The same modes for the Krylov algorithm are available as described in
:py:func:`exactdiag.PureState.propagate_krylov`. Since the
effective Hamiltonian is not hermitian, the Arnoldi version of the
Krylov algorithm is used.
"""
try:
tmp = self.qtnorm
except AttributeError:
self.qtnorm = self.norm()
self.qtr = np.random.rand(1)[0]
ConvParam = param['Quenches'].data[quenchtime[0]]
ConvParam = ConvParam['ConvergenceParameters'].data[0]
mode = ConvParam['edlibmode']
if(mode < 0):
if(self.Psi.shape[0] < 2**10):
# Complete H, built-in scipy
mode = 0
elif(self.Psi.shape[0] < 2**10):
# Complete H, EDLib Krylov (not used as default)
mode = 1
elif(self.Psi.shape[0] < 2**15):
# MPO, EDLib Krylov
mode = 2
else:
# MPO, EDLib Krylov, write Krylov vectors to hard drive
mode = 3
if(mode == 0):
# Complete H, built-in scipy
# --------------------------
if(self.lastprop is None):
self.lastprop = Hmpo.build_effhamiltonian(param, self.SymmSec,
True, maxdim,
quenchtime=quenchtime)
self.Psi = spla.expm_multiply(-1j * dt * self.lastprop, self.Psi)
self.quantum_jump(Hmpo, param, quenchtime)
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
nvec = ConvParam['max_num_lanczos_iter']
tol = ConvParam['lanczos_tol']
numzero = 1e-16
hmat = np.zeros((nvec, nvec), dtype=np.complex128)
norm = np.sqrt(np.real(np.sum(np.conj(self.Psi) * self.Psi)))
beta = deepcopy(norm)
# Set nn in case if-break conditions are never reached
nn = nvec
if((mode == 1) or (mode == 2)):
# EDLib Krylov (Complete H, MPO dot product)
# ------------
base = np.zeros((nvec, self.Psi.shape[0]), dtype=np.complex128)
if((mode == 1) and (self.lastprop is None)):
self.lastprop = Hmpo.build_effhamiltonian(param, self.SymmSec,
sparse, maxdim,
quenchtime=quenchtime)
nn = nvec
for jj in range(nvec):
self.Psi = self.Psi / norm
base[jj, :] = self.Psi
if(mode == 1):
self.Psi = self.lastprop.dot(self.Psi)
else:
HPsi = self.dot_effhmpo_vec(Hmpo, Operators, param,
quenchtime, maxdim)
for ii in range(jj):
hmat[ii, jj] = np.sum(np.conj(self.Psi) * base[ii, :])
self.Psi -= hmat[ii, jj] * base[jj, :]
norm = np.sqrt(np.real(np.sum(np.conj(self.Psi) * self.Psi)))
hmat[jj + 1, jj] = norm
prop = sla.expm(-1j * dt * hmat[:jj + 1, :jj + 1])
est = 2.0 * norm * np.abs(prop[jj, 0])
if(est < tol):
nn = jj
break
if(abs(beta[jj + 1]) < 10 * numzero):
nn = jj
break
self.Psi = HPsi
self.Psi = prop[0, 0] * base[:, 0]
for jj in range(1, nn):
self.Psi += prop[jj, 0] * base[:, jj]
else:
# MPO, EDLib Krylov, write Krylov vectors to hard drive
# -----------------------------------------------------
# Save vectors to hard disk to save memory
fbase = param['Write_Directory'] + param['job_ID'] + \
param['unique_ID'] + 'Krylovbase_'
for jj in range(nvec):
self.Psi = self.Psi / norm
np.save(fbase + str(jj) + '.npy', self.Psi)
self.Psi = self.dot_effhmpo_vec(Hmpo, Operators, param,
quenchtime, maxdim)
for ii in range(jj):
HPsi = np.load(fbase + str(ii) + '.npy')
hmat[ii, jj] = np.sum(np.conj(self.Psi) * HPsi)
self.Psi -= hmat[ii, jj] * HPsi
norm = np.sqrt(np.real(np.sum(np.conj(self.Psi) * self.Psi)))
hmat[jj + 1, jj] = norm
prop = sla.expm(-1j * dt * hmat[:jj + 1, :jj + 1])
est = 2.0 * norm * np.abs(prop[jj, 0])
if(est < tol):
nn = jj
del HPsi
break
if(abs(beta[jj + 1]) < 10 * numzero):
nn = jj
del HPsi
break
self.Psi = HPsi
self.Psi = np.load(fbase + str(0) + '.npy')
self.Psi = beta * prop[0, 0] * self.Psi
for jj in range(1, nn):
flnm = fbase + str(jj) + '.npy'
tmp = np.load(flnm)
remove_file(flnm)
self.Psi += beta * prop[jj, 0] * tmp
self.quantum_jump(Hmpo, param, quenchtime)
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
[docs] def propagate_gates(self, func, args=(), return_copy=False):
"""
Propagate a quantum state with a set of given gates. At present only
2-site gates are supported.
**Arguments**
func : iterator (``__iter__``)
return gates and left most site index where gate should be
applied.
args : tuple, optional
arguments for call to the iterator.
Default to ()
return_copy : boolean, optional
If ``False``, then the state of the class is modified/propagated.
The state remains unchanged and a copy of the state after applying
the sequence of gates is returned for ``True``.
Default to ``False``.
"""
if(return_copy):
psicopy = deepcopy(self.Psi)
self._trotter_initialize(False)
for gate, idx in func(*args):
if(gate.shape[0] != self.ld**2):
raise Exception('Only two site gates supported at present.')
self.apply_op2(gate, idx, perm=True)
self._trotter_finalize(False, None, None, None)
if(return_copy):
self.Psi, psicopy = psicopy, self.Psi
return psicopy
return
# ------------------- Trotter decomposition functions ----------------------
[docs] def _trotter_initialize(self, qt):
"""
The initialization transfers the vector into a rank-L tensor if
no symmetries are present. Furthermore the norm and random number
for the quantum trajectories are created before the first time
step.
**Arguments**
qt : bool
``True`` if evolution with quantum trajectories, otherwise
``False``.
"""
if(qt):
try:
tmp = self.qtnorm
except AttributeError:
self.qtnorm = self.norm()
self.qtr = np.random.rand(1)[0]
if(self.SymmSec is None):
# Only wave vectors without symmetry are reshaped
self.Psi = np.reshape(self.Psi, [self.ld] * self.ll)
return
[docs] def _trotter_finalize(self, qt, Hmpo, param, quenchtime):
"""
At the end of the Trotter step, the rank-L tensor has to be
converted back into the state vector if no symmetry is present.
For quantum trajectories it is checked if a jump has to be applied.
**Arguments**
qt : bool
``True`` if evolution with quantum trajectories, otherwise
``False``.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time.
"""
if(self.SymmSec is None):
# Only wave vectors without symmetry were reshaped
self.Psi = np.reshape(self.Psi, (self.ld**self.ll))
if(qt): self.quantum_jump(Hmpo, param, quenchtime)
return
[docs] def trotter_layer(self, dt, Hmpo, Operators, param, quenchtime,
even, qt):
"""
Apply an even or odd layer of the Trotter decomposition to the
state vector.
**Arguments**
dt : float
time step/discretization of the simulation.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
Operators : instance of :py:class:`Ops.OperatorList`
Contains the operators used to build the Hamiltonian.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time.
even : bool
Boolean to switch between even (``True``) and odd (``False``)
layers in the Trotter decomposition.
qt : bool
Flag if using Schrodinger equation (``False``) or quantum
trajectories (``True``).
"""
if(qt):
# Consider the effective Hamiltonian
self.trotter_layer_qt(dt, Hmpo, Operators, param, quenchtime, even)
return
if((self.bidx is None) and (not self.SymmSec is None)):
self.bidx = idx_symmetry_blocks(Operators, param)
for elem in self.trotter_Ham_2site(Hmpo, Operators, param, quenchtime,
even):
# Build local propagator and apply
if((elem[1], dt) in self.lastprop):
Op2 = self.lastprop[(elem[1], dt)]
else:
Op2 = self._matrix_exp_trotter(dt * elem[0])
self.lastprop[(elem[1], dt)] = Op2
self.apply_op2(Op2, elem[1], ll_1_term=elem[2])
self.trotter_layer_perm(even, param['L'], Hmpo.pbc)
return
[docs] def trotter_layer_qt(self, dt, Hmpo, Operators, param, quenchtime, even):
"""
Trotter layer under the effective Hamiltonian of the Lindblad master
equation representing the hermitian part for an even or odd layer.
**Arguments**
dt : float
time step/discretization of the simulation.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class.
Operators : instance of :py:class:`Ops.OperatorList`
Contains the operators used to build the Hamiltonian.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time.
even : bool
Boolean to switch between even (``True``) and odd (``False``)
layers in the Trotter decomposition.
"""
ld = Operators['I'].shape[0]
ll = param['L']
if((self.bidx is None) and (not self.SymmSec is None)):
self.bidx = idx_symmetry_blocks(Operators, param)
for elem in self.trotter_Ham_2site(Hmpo, Operators, param, quenchtime,
even):
jj = elem[1]
if((jj, dt) in self.lastprop):
Prop = self.lastprop[(jj, dt)]
else:
Ham = elem[0].astype(np.complex128)
jjp1 = (jj + 1)%ll
for ii in range(len(Hmpo.data['lind1']['hp'])):
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
if(Hmpo.pbc):
hp *= 0.5
else:
hp[1:-1] = 0.5 * np.array(hp[1:-1])
# imag j disappears with -1j in exponential
Op = Operators[Hmpo.data['lind1']['Op'][ii]]
Op = -0.5j * np.dot(Op.conj().transpose(), Op)
Ham += np.kron(hp[jj] * Op, np.eye(ld))
Ham += np.kron(np.eye(ld), hp[jjp1] * Op)
Prop = self._matrix_exp_trotter(dt * Ham)
self.lastprop[(jj, dt)] = Prop
self.apply_op2(Prop, jj, elem[2])
self.trotter_layer_perm(even, ll, Hmpo.pbc)
return
[docs] def trotter_layer_perm(self, even, ll, pbc):
"""
Return permutation to be applied after a complete layer of two-site
Trotter terms.
**Arguments**
even : bool
Indicates if the even or odd layer of the Trotter decomposition
was applied.
ll : int
System size.
pbc : bool
Indicates if the Hamiltonian has periodic boundary conditions or
not.
"""
if(self.SymmSec is not None): return
# Boolean for odd-even system size
lleven = (ll%2 == 0)
notpbc = not pbc
# Eight cases
if(notpbc and lleven and even):
# Even layer, even system size, no PBC
nn = ll // 2 - 1
idx = [[2 * ii, 2 * ii + 1] for ii in range(nn - 1, -1, -1)]
tmp = [elem for sublist in idx for elem in sublist]
perm = [ll - 2] + tmp + [ll - 1]
elif(notpbc and lleven):
# Odd layer, even system size, no PBC
nn = ll // 2
idx = [[2 * ii, 2 * ii + 1] for ii in range(nn - 1, -1, -1)]
perm = [elem for sublist in idx for elem in sublist]
elif(notpbc and even):
# Even layer, odd system size, no PBC
nn = ll // 2
idx = [[2 * ii, 2 * ii + 1] for ii in range(nn - 1, -1, -1)]
perm = [ll - 1] + [elem for sublist in idx for elem in sublist]
elif(notpbc):
# Odd layer, odd system size, no PBC
nn = ll // 2
idx = [[2 * ii, 2 * ii + 1] for ii in range(nn - 1, -1, -1)]
perm = [elem for sublist in idx for elem in sublist] + [ll - 1]
elif(lleven and even):
# Even layer, even system size, PBC
nn = ll // 2
idx = [[2 * ii, 2 * ii + 1] for ii in range(nn - 1, 0, -1)]
perm = [1] + [elem for sublist in idx for elem in sublist] + [0]
elif(lleven):
# Odd layer, even system size, PBC (not affected by PCB right now)
nn = ll // 2
idx = [[2 * ii, 2 * ii + 1] for ii in range(nn - 1, -1, -1)]
perm = [elem for sublist in idx for elem in sublist]
else:
# Odd layer, odd system size, PBC
# Even layer, odd system size, PBC
raise Exception("Periodic boundary only with even system size!")
self.Psi = np.transpose(self.Psi, perm)
return
[docs] def _matrix_exp_trotter(self, Ham):
"""
Calculate the matrix exponential considering block diagonal structure
if symmetries are present.
**Arguments**
Ham : np 2d array
matrix representing the Hamiltonian including the time step. To
be exponentiated.
**Details**
Functions accesses ``self.bidx``, a dictionary of tuples or ``None``.
Each set of tuples contains the indices for the same symmetry sector
and the corresponding site indices.
"""
if(self.bidx is None):
# No symmetries
return sla.expm(-1j * Ham)
# Symmetries using bidx
dic = {}
for key in self.bidx.keys():
Sub = sla.expm(-1j * Ham[key, :][:, key])
for i1 in range(len(key)):
j1 = self.bidx[key][i1]
for i2 in range(len(key)):
j2 = self.bidx[key][i2]
if(Sub[i1, i2] == 0.0): continue
try:
dic[j1].append((j2[0], j2[1], Sub[i1, i2]))
except KeyError:
dic[j1] = [(j2[0], j2[1], Sub[i1, i2])]
return dic
[docs] def _apply_op2_nosymm(self, Op, ii, ll_1_term, perm):
"""
Apply a two-site operator to the state vector on a specified site.
**Arguments**
Op : np 2d array
The operator acting on two sites represented as d^2 x d^2 matrix.
ii : int
Apply term to the sites ii and (ii + 1).
ll_1_term : bool, optional
Boolean to indicate the term ii=L for periodic boundary
conditions.
"""
Op = np.reshape(Op, (self.ld, self.ld, self.ld, self.ld))
if(ll_1_term and perm):
# 1st and last site have indices 0 and ll - 1
perm = [1] + list(range(2, self.ll)) + [0]
c1 = self.ll - 1
c2 = 0
elif(ll_1_term):
# PBC setting for even system size, application of all other
# terms before. 1st and last site have indices ll - 2, ll - 1
c1 = self.ll - 1
c2 = self.ll - 2
else:
c1 = ii
c2 = ii + 1
self.Psi = np.tensordot(Op, self.Psi, ([2, 3], [c1, c2]))
if(not perm): return
if(not ll_1_term):
perm = list(range(2, ii + 2)) + [0, 1] + \
list(range(ii + 2, self.ll))
self.Psi = np.transpose(self.Psi, perm)
return
[docs] def _apply_op2_symm(self, Op, ii, ll_1_term, *args):
"""
Apply a two-site operator to the state vector on a specified site.
**Arguments**
Op : np 2d array
The operator acting on two sites represented as d^2 x d^2 matrix.
ii : int
Apply term to the sites ii and (ii + 1).
ll_1_term : bool, optional
Boolean to indicate the term ii=L for periodic boundary
conditions.
"""
Phi = np.zeros((self.Psi.shape[0]), dtype=np.complex128)
if(ll_1_term):
c1 = 0
c2 = self.ll - 1
else:
c1 = ii
c2 = ii + 1
for k1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(k1)
idx2 = deepcopy(idx1)
try:
elems = Op[(idx1[c1], idx1[c2])]
except KeyError:
continue
for kl, kr, elem in elems:
idx2[c1] = kl
idx2[c2] = kr
try:
k2 = self.SymmSec[tuple(idx2)]
except KeyError:
continue
Phi[k1] += elem * self.Psi[k2]
self.Psi = Phi
return
[docs] def _apply_op_nosymm(self, Op, idx):
"""
Apply an operator/gate on an arbitrary number of sites.
**Arguments**
Op : np 2d array
The operator acting on the site represented as matrix.
idx : list of integers
The sites to be contracted with the operator.
"""
# Number of sites acted on
nn = len(idx)
# Get shape of tensor representation of Op
shape = tuple([self.ld] * 2 * nn)
Op_ = np.reshape(Op, shape)
# Get list of contraction indices for Op
contr = list(range(nn, 2 * nn))
self.Psi = np.tensordot(Op_, self.Psi, (contr, idx))
perm = []
pos = nn
for ii in range(self.ll):
if(ii in idx):
for jj in range(nn):
if(ii == idx[jj]):
perm.append(jj)
break
else:
perm.append(pos)
pos += 1
self.Psi = np.transpose(self.Psi, perm)
return
[docs] def _apply_op_symm(self, Op, idx):
"""
Apply an operator/gate on an arbitrary number of sites to a
symmetry conserving state.
**Arguments**
Op : np 2d array
The operator acting on the site represented as matrix.
idx : list of integers
The sites to be contracted with the operator.
"""
# 1) Extract operator (pre-calculating this would certainly
# be more efficient than doing it every time)
dic = {}
for key in self.bidx.keys():
Sub = Op[key, :][:, key]
for i1 in range(len(key)):
j1 = self.bidx[key][i1]
for i2 in range(len(key)):
j2 = self.bidx[key][i2]
if(Sub[i1, i2] == 0.0): continue
val = tuple(list(j2) + [Sub[i1, i2]])
try:
dic[j1].append(val)
except KeyError:
dic[j1] = [val]
# 2) Start matrix-vector multiplication
Phi = np.zeros((self.Psi.shape[0]), dtype=np.complex128)
for k1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(k1)
idx2 = deepcopy(idx1)
try:
elems = dic[tuple(idx1[idx])]
except KeyError:
continue
for entry in elems:
idx2[idx] = entry[:len(idx)]
elem = entry[-1]
try:
k2 = self.SymmSec[tuple(idx2)]
except KeyError:
continue
Phi[k1] += elem * self.Psi[k2]
self.Psi = Phi
return
[docs] def densitymatrix(self):
"""
Convert the pure state into a density matrix.
"""
return DensityMatrix(self.Psi, self.ld, self.ll, self.SymmSec)
[docs] def distance(self, State, dist):
"""
Measure the distance between the pure state and another quantum
state.
**Arguments**
State : 1d or 2d numpy array
Represents another pure state (1d) or a density matrix (2d).
distance : str
Define the type of distance measure. Either ``Trace`` or
``Infidelity``.
"""
return qi.distance(self.Psi, State, dist)
[docs] def fidelity(self, Phi):
"""
Return the fidelity between this state Psi and Phi defined as
:math:`|\langle \psi | \phi \\rangle |^2`.
**Arguments**
Phi : np 1d array or instance of :py:class:`exactdiag.PureState`
Other vector in the overlap with the PureState.
"""
if(isinstance(Phi, PureState)):
overlap = self.dot_vec(Phi)
return np.conj(overlap) * overlap
elif(isinstance(Phi, DensityMatrix)):
return 1 - self.distance(Phi.rho, 'Infidelity')
elif(len(Phi.shape) == 1):
overlap = self.dot_vec(Phi)
return np.conj(overlap) * overlap
elif(len(Phi.shape) == 2):
return 1 - self.distance(Phi, 'Infidelity')
else:
raise Exception("Unknown type of quantum state.")
[docs] def purity(self):
"""
The purity of a pure quantum state is always 1.
"""
return 1.0
# ==============================================================================
# METHODS FOR OPEN SYSTEMS / DENSITY MATRICES
# ==============================================================================
[docs]class DensityMatrix(_QuantumSystem):
"""
Density matrix for evolution of open systems.
**Arguments**
Arr : 1d or 2d np array
If 1d array, a state vector Psi is assumend and the density matrix
is built as :math:`| \psi \\rangle \langle \psi |`. If 2d array, array
is assumed to be the density matrix.
ld : int
local dimension of the Hilbert space for each site.
ll : int
number of sites in the system.
SymmSec: None or instance of :py:class:`exactdiag.SymmetrySector`
contains the setup for symmetries. If no symmetry is present, pass
``None``.
**Variables**
D : int
Dimension of the Hilbert space
rho : 2d numpy array
Density matrix
rhos : dict
Storing reduced density matrices
lastprop : None, 2d numpy array, or dictionary.
Last propagator used for time independent quenches. If 2d numpy array,
the propagator is stored as matrix. If dictionary, the local
propagators for the Trotter decomposition are stored.
bidx : dictionary
Used to store the indices of the subblocks with the same
symmetry numbers.
"""
def __init__(self, Arr, ld, ll, SymmSec):
self.ld = ld
self.ll = ll
self.SymmSec = SymmSec
self.D = Arr.shape[0]
if(len(Arr.shape) == 1):
# State vector - outer product
self.rho = np.dot(np.reshape(Arr, (self.D, 1)),
np.reshape(np.conj(Arr), (1, self.D)))
elif(len(Arr.shape) == 2):
# Density matrix
self.rho = Arr
else:
raise ValueError('Higher than rank-2 tensor for density matrix.')
self.rhos = {}
self.lastprop = None
self.bidx = None
[docs] def save(self, flnm):
np.save(flnm, self.rho)
# --------------------- DENSITY MATRICES -----------------------------------
[docs] def _traceout_nosymm(self, trce, keep):
"""
Trace over a set of sites specified via their indices for states
without a symmetry.
**Arguments**
trce : list of ints
Trace out over the indices contain in this list.
keep : list of ints
Keep the sites whose indices are listed here. Together trce
and keep should list all indices exactly one time.
"""
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
if(isinstance(self.ld, int)):
shape1 = [self.ld] * (2 * self.ll)
d1 = self.ld**(len(keep))
d2 = self.ld**(len(trce))
perm = keep + list(np.array(keep) + self.ll) \
+ trce + list(np.array(trce) + self.ll)
tmp = np.reshape(np.transpose(np.reshape(self.rho, shape1),
perm), (d1, d1, d2**2))
mask = np.cumsum((d2 + 1) * np.ones(d2), dtype=np.int32) - (d2 + 1)
rho = np.sum(tmp[:, :, mask], axis=2)
# Alternative
#np.sum(np.diag(tmp[kk, jj, :, :]))
return DensityMatrix(rho, self.ld, len(keep), None)
else:
shape1 = self.ld + self.ld
d1 = np.prod(self.ld[keep])
d2 = np.prod(self.ld[trce])
perm = keep + list(np.array(keep) + self.ll) \
+ trce + list(np.array(keep) + self.ll)
tmp = np.reshape(np.transpose(np.reshape(self.rho, shape1),
perm), (d1, d1, d2**2))
mask = np.cumsum((d2 + 1) * np.ones(d2), dtype=np.int32) - (d2 + 1)
rho = np.sum(tmp[:, :, mask], axes=2)
return DensityMatrix(rho, self.ld[keep], len(keep), None)
[docs] def _traceout_symm2nosymm(self, trce, keep, SubSec=None):
"""
Trace over a set of sites specified via their indices for states
with symmetry returning a density matrix not considering the
symmetry.
**Arguments**
trce : list of ints
Trace out over the indices contain in this list.
keep : list of ints
Keep the sites whose indices are listed here. Together trce
and keep should list all indices exactly one time.
**Details**
This functions is unefficient when tracing over single sites (or
small subsystems).
"""
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
if(SubSec is None):
SubSec = self.SymmSec.get_subspace(keep)
TrcSec = self.SymmSec.get_subspace(trce)
if(isinstance(self.ld, int)):
d1 = self.ld**(len(keep))
rho = np.zeros((d1, d1), dtype=self.rho.dtype)
# Identify position in reduced rho
base = self.ld**np.linspace(len(keep) - 1, 0, len(keep))
# Identifying element in rho
idx1 = np.zeros((self.ll), dtype=np.int32)
idx2 = np.zeros((self.ll), dtype=np.int32)
for kk in range(len(SubSec)):
idx1[keep] = SubSec(kk)
kkk = int(np.sum(base * SubSec(kk)))
for jj in range(len(SubSec)):
idx2[keep] = SubSec(jj)
jjj = int(np.sum(base * SubSec(jj)))
for ii in range(len(TrcSec)):
try:
idx1[trce] = TrcSec(ii)
idx2[trce] = TrcSec(ii)
i1 = self.SymmSec[list(idx1)]
i2 = self.SymmSec[list(idx2)]
rho[kkk, jjj] += self.rho[i1, i2]
except KeyError:
pass
return DensityMatrix(rho, self.ld, len(keep), None)
else:
d1 = np.prod(self.ld[keep])
rho = np.zeros((d1, d2), dtype=self.rho.dtype)
# Identify position in reduced rho
base = np.cumprod(self.ld[keep::-1])[::-1] // self.ld[keep[-1]]
# Identifying element in rho
idx1 = np.zeros((self.ll), dtype=np.int32)
idx2 = np.zeros((self.ll), dtype=np.int32)
for kk in range(len(SubSec)):
idx1[keep] = SubSec(kk)
kkk = np.sum(base * SubSec(kk))
for jj in range(len(SubSec)):
idx2[keep] = SubSec(jj)
jjj = np.sum(base * SubSec(jj))
for ii in range(len(TrcSec)):
try:
idx1[trce] = TrcSec(ii)
idx2[trce] = TrcSec(ii)
i1 = self.SymmSec[list(idx1)]
i2 = self.SymmSec[list(idx2)]
rho[kkk, jjj] += self.rho[i1, i2]
except KeyError:
pass
return DensityMatrix(rho, self.ld[keep], len(keep), None)
[docs] def _traceout_symm(self, trce, keep, SubSec=None):
"""
Trace over a set of sites specified via their indices for states
with symmetry returning a density matrix considering the symmetry
in the larger Hilbert space.
**Arguments**
trce : list of ints
Trace out over the indices contain in this list.
keep : list of ints
Keep the sites whose indices are listed here. Together trce
and keep should list all indices exactly one time.
SubSec : instance of :py:class:`exactdiag.SymmetrySector`
If symmetry sector is already calculated, it can be passed
with this argument. Otherwise pass ``None``.
Default to ``None``.
"""
if(keep is None):
mask = np.ones((self.ll), dtype=np.bool)
mask[trce] = False
keep = list(np.array(list(range(self.ll)))[mask])
if(SubSec is None):
SubSec = self.SymmSec.get_subspace(keep)
TrcSec = self.SymmSec.get_subspace(trce)
d1 = len(SubSec)
d2 = len(TrcSec)
# Identifying element in rho
idx1 = np.zeros((self.ll), dtype=np.int32)
idx2 = np.zeros((self.ll), dtype=np.int32)
rho = np.zeros((d1, d1), dtype=self.rho.dtype)
for kk in range(d1):
idx1[keep] = SubSec(kk)
for jj in range(d1):
idx2[keep] = SubSec(jj)
for ii in range(d2):
idx1[trce] = TrcSec(ii)
idx2[trce] = TrcSec(ii)
try:
i1 = self.SymmSec[list(idx1)]
i2 = self.SymmSec[list(idx2)]
rho[kk, jj] += self.rho[i1, i2]
except KeyError:
pass
if(isinstance(self.ld, int)):
return DensityMatrix(rho, self.ld, len(keep), SubSec)
else:
return DensityMatrix(rho, self.ld[keep], len(keep), SubSec)
# ------------------------ ENTROPY -----------------------------------------
[docs] def entropy(self):
"""
Calculate the von-Neumann entropy of the density matrix.
"""
numzero = 1e-15
vals = sla.eigh(self.rho, eigvals_only=True)
be = - np.sum(vals[vals > numzero] * np.log(vals[vals > numzero]))
return be
# ---------------------- PARTIAL TRACE FOR BOND ENTROPY --------------------
[docs] def traceout_first(self):
"""
Trace-out the first site in the density matrix.
"""
if(self.SymmSec is None):
return self._traceout_first_nosymm()
SubSec = self.SymmSec.get_subspace(list(range(1, self.ll)))
if(len(SubSec) > 0.75 * self.ld**self.ll):
return self._traceout_first_symm2nosymm(SubSec)
return self._traceout_first_symm(SubSec)
[docs] def _traceout_first_nosymm(self):
"""
Return the reduced density matrix tracing over the first site in the
system. This functions treats density matrices without symmetry.
"""
ld = self.ld
dim = ld**(self.ll - 1)
rho2nn = np.zeros((dim, dim), dtype=self.rho.dtype)
tmp = np.reshape(np.transpose(np.reshape(self.rho, (ld, dim, ld, dim)),
(1, 3, 0, 2)), (dim, dim, ld**2))
mask = np.cumsum((ld + 1) * np.ones(ld, dtype=int)) - (ld + 1)
for kk in range(dim):
for jj in range(dim):
rho2nn[kk, jj] = np.sum(tmp[kk, jj, mask])
if(np.trace(rho2nn) > 1.0 + 1e-11):
raise ValueError('Normalization violated with %3.6e!'%(np.trace(rho2nn)))
return DensityMatrix(rho2nn, ld, (self.ll - 1), None)
[docs] def _traceout_first_symm(self, SubSec):
"""
Return the reduced density matrix tracing over the first site in the
system. This function treats density matrices with symmetry and
calculates the new density matrix in a symmetry adapted space.
**Arguments**
SubSec : instance of :py:class:`exactdiag.SymmetrySector`
This is the symmetry sector of the system from site 2 to site
L (end of the system). It is calculated before hand in order
to estimate if a symmetry adapted basis is appropriate.
"""
dim = len(SubSec)
rho2nn = np.zeros((dim, dim), dtype=self.rho.dtype)
for kk in range(len(SubSec)):
for jj in range(len(SubSec)):
for ii in range(self.ld):
try:
i1 = self.SymmSec[[ii] + list(SubSec(kk))]
i2 = self.SymmSec[[ii] + list(SubSec(jj))]
rho2nn[kk, jj] += self.rho[i1, i2]
except KeyError:
pass
if(np.trace(rho2nn) > 1.0 + 1e-11):
raise ValueError('Normalization violated with %3.6e!'%(np.trace(rho2nn)))
return DensityMatrix(rho2nn, self.ld, (self.ll - 1), SubSec)
[docs] def _traceout_first_symm2nosymm(self, SubSec):
"""
Return the reduced density matrix tracing over the first site in the
system. This function treats density matrices with symmetry and
calculates the new density matrix in the complete Hilbert space.
**Arguments**
SubSec : instance of :py:class:`exactdiag.SymmetrySector`
This is the symmetry sector of the system from site 2 to site
L (end of the system). It is calculated before hand in order
to estimate if a symmetry adapted basis is appropriate.
"""
ld = self.ld
dim = ld**(self.ll - 1)
rho2nn = np.zeros((dim, dim), dtype=self.rho.dtype)
base = ld**np.linspace(self.ll - 2, 0, self.ll - 1)
for kk in range(len(SubSec)):
for jj in range(len(SubSec)):
for ii in range(ld):
try:
i1 = self.SymmSec[[ii] + list(SubSec(kk))]
i2 = self.SymmSec[[ii] + list(SubSec(jj))]
kkk = np.sum(base * SubSec(kk))
jjj = np.sum(base * SubSec(jj))
rho2nn[kkk, jjj] += self.rho[i1, i2]
except KeyError:
pass
if(np.abs(1 - np.trace(rho2nn)) > 1e-11):
raise ValueError('Normalization violated with '
+ '%3.6e!'%(np.abs(1 - np.trace(rho2nn))))
return DensityMatrix(rho2nn, ld, (self.ll - 1), None)
[docs] def traceout_last(self):
"""
Trace-out the last site of the density matrix.
"""
if(self.SymmSec is None):
return self._traceout_last_nosymm()
SubSec = self.SymmSec.get_subspace(list(range(0, self.ll - 1)))
if(len(SubSec) > 0.75 * self.ld**self.ll):
return self._traceout_last_symm2nosymm(SubSec)
return self._traceout_last_symm(SubSec)
[docs] def _traceout_last_nosymm(self):
"""
Return the reduced density matrix tracing over the last site in the
system. This functions treats density matrices without symmetry.
"""
ld = self.ld
dim = ld**(self.ll - 1)
rho2nn = np.zeros((dim, dim), dtype=self.rho.dtype)
tmp = np.reshape(np.transpose(np.reshape(self.rho, (dim, ld, dim, ld)),
(0, 2, 1, 3)), (dim, dim, ld**2))
mask = np.cumsum((ld + 1) * np.ones(ld, dtype=int)) - (ld + 1)
for kk in range(dim):
for jj in range(dim):
rho2nn[kk, jj] = np.sum(tmp[kk, jj, mask])
if(np.abs(1 - np.trace(rho2nn)) > 1e-11):
raise ValueError('Normalization violated with '
+ '%3.6e!'%(np.abs(1 - np.trace(rho2nn))))
return DensityMatrix(rho2nn, ld, (self.ll - 1), None)
[docs] def _traceout_last_symm(self, SubSec):
"""
Return the reduced density matrix tracing over the last site in the
system. This function treats density matrices with symmetry and
calculates the new density matrix in a symmetry adapted space.
**Arguments**
SubSec : instance of :py:class:`exactdiag.SymmetrySector`
This is the symmetry sector of the system from site 1 to site
(L - 1) (L is the end of the system). It is calculated before
hand in order to estimate if a symmetry adapted basis is
appropriate.
"""
dim = len(SubSec)
rho2nn = np.zeros((dim, dim), dtype=self.rho.dtype)
for kk in range(len(SubSec)):
for jj in range(len(SubSec)):
for ii in range(self.ld):
try:
i1 = self.SymmSec[list(SubSec(kk)) + [ii]]
i2 = self.SymmSec[list(SubSec(jj)) + [ii]]
rho2nn[kk, jj] += self.rho[i1, i2]
except KeyError:
pass
if(np.abs(1 - np.trace(rho2nn)) > 1e-11):
raise ValueError('Normalization violated with '
+ '%3.6e!'%(np.abs(1 - np.trace(rho2nn))))
return DensityMatrix(rho2nn, self.ld, (self.ll - 1), SubSec)
[docs] def _traceout_last_symm2nosymm(self, SubSec):
"""
Return the reduced density matrix tracing over the last site in the
system. This function treats density matrices with symmetry and
calculates the new density matrix in the complete Hilbert space.
**Arguments**
SubSec : instance of :py:class:`exactdiag.SymmetrySector`
This is the symmetry sector of the system from site 1 to site
(L - 1) (L is the end of the system). It is calculated before
hand in order to estimate if a symmetry adapted basis is
appropriate.
"""
ld = self.ld
dim = ld**(self.ll - 1)
rho2nn = np.zeros((dim, dim), dtype=self.rho.dtype)
base = ld**np.linspace(self.ll - 2, 0, self.ll - 1)
for kk in range(len(SubSec)):
for jj in range(len(SubSec)):
for ii in range(ld):
try:
i1 = self.SymmSec[list(SubSec(kk)) + [ii]]
i2 = self.SymmSec[list(SubSec(jj)) + [ii]]
kkk = np.sum(base * SubSec(kk))
jjj = np.sum(base * SubSec(jj))
rho2nn[kkk, jjj] += self.rho[i1, i2]
except KeyError:
pass
if(np.abs(1 - np.trace(rho2nn)) > 1e-11):
raise ValueError('Normalization violated with '
+ '%3.6e!'%(np.abs(1 - np.trace(rho2nn))))
return DensityMatrix(rho2nn, ld, (self.ll - 1), None)
[docs] def meas_mpo(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Measure a MPO using the rule set.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class or MPO representing
a measurement.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
maxdim : int
maximal dimension of the Hilbert space. Only referenced in ...
"""
Hrho = self.dot_hmpo_mat(Hmpo, Operators, param, quenchtime, maxdim)
return np.trace(Hrho)
[docs] def dot_hmpo_mat(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a Hamiltonian represented as MPO / MPO rule set with a
density matrix.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class or MPO representing
a measurement.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
maxdim : int
maximal dimension of the Hilbert space. Only referenced in ...
"""
if(self.SymmSec is None):
return self._dot_hmpo_mat_nosymm(Hmpo, Operators, param,
quenchtime, maxdim)
return self._dot_hmpo_mat_symm(Hmpo, Operators, param,
quenchtime, maxdim)
[docs] def _dot_hmpo_mat_nosymm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a MPO which can represent the Hamiltonian or a measurement
with the density matrix. For arguments look into documention of
:py:func:`exactdiag.DensityMatrix.dot_hmpo_mat`. This is the version
for density matrices without symmetry adapted basis. Intended for
internal use.
"""
# Get shortcuts to system size and local dimension of Hilbert space
ll = param['L']
ld = Operators['I'].shape[0]
dim = ld**ll
shape = [ld] * ll + [dim]
self.rho = np.reshape(self.rho, shape)
Phi = np.zeros((shape), dtype=self.rho.dtype)
# Loop over local terms
for ii in range(len(Hmpo.data['site'])):
hp = Hmpo.data['site'].hparam(ii, param, quenchtime)
Op = Operators[Hmpo.data['site']['Op'][ii]]
for jj in range(ll):
perm = list(range(1, jj + 1)) + [0] + \
list(range(jj + 1, ll + 1))
Phi += np.transpose(np.tensordot(hp[jj] * Op, self.rho,
([1], [jj])), perm)
# Loop over bond terms
for ii in range(len(Hmpo.data['bond'])):
hp = Hmpo.data['bond'].hparam(ii, param, quenchtime)
Op = [Operators[Hmpo.data['bond']['Op'][ii][0]],
Operators[Hmpo.data['bond']['Op'][ii][1]]]
for jj in range(ll - 1):
perm = list(range(2, jj + 2)) + [1, 0] + \
list(range(jj + 2, ll + 1))
tmp = np.tensordot(hp[jj] * Op[0], self.rho,
([1], [jj]))
Phi += np.transpose(np.tensordot(Op[1], tmp,
([1], [jj + 1])), perm)
if(Hmpo.pbc):
perm = [0] + list(range(2, ll)) + [1, ll]
tmp = np.tensordot(hp[ll - 1] * Op[0], self.rho,
([1], [ll - 1]))
Phi += np.transpose(np.tensordot(Op[1], tmp, ([1], [1])),
perm)
# Loop over exponential terms
for ii in range(len(Hmpo.data['exponential'])):
hp = Hmpo.data['exponential'].hparam(ii, param, quenchtime)
Op = [Operators[Hmpo.data['exponential']['Op'][ii][0]],
Operators[Hmpo.data['exponential']['Op'][ii][1]],
Operators[Hmpo.data['exponential']['Op'][ii][2]]]
dp = Hmpo.data['exponential']['dp'][ii]
# Check if phase is not identity
phase = (Hmpo.data['exponential']['Op'][ii][2] != 'I')
if(phase):
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
perm = list(range(dist + 1, jj + dist + 1)) + \
list(range(dist + 1))[::-1] + \
list(range(jj + dist + 1, ll + 1))
tmp = np.tensordot(hp[jj] * Op[0], self.rho,
([1], [jj]))
for nn in range(dist - 1):
tmp = np.tensordot(dp * Op[2], tmp,
([1], [jj + 1 + nn]))
Phi += np.transpose(np.tensordot(Op[1], tmp,
([1], [jj + dist])),
perm)
else:
# No phase - identities
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
perm = list(range(2, jj + 2)) + [1] \
+ list(range(jj + 2, jj + dist + 1)) + [0] \
+ list(range(jj + dist + 1, ll + 1))
tmp = np.tensordot(hp[jj] * Op[0], self.rho,
([1], [jj]))
tmp = np.tensordot(dp**(dist - 1) * Op[1], tmp,
([1], [jj + dist]))
Phi += np.transpose(tmp, perm)
if(len(Hmpo.data['product']) > 0):
raise NotImplementedError('Product terms not implemented for ED.')
if(len(Hmpo.data['FiniteFunction']) > 0):
raise NotImplementedError('FiniteFunction terms not implemented' +
' for ED.')
# Loop over InfiniteFunction terms
key = 'InfiniteFunction'
for ii in range(len(Hmpo.data[key])):
hp = Hmpo.data[key].hparam(ii, param, quenchtime)
Opl = Operators[Hmpo.data[key]['Op'][ii][0]]
Opr = Operators[Hmpo.data[key]['Op'][ii][1]]
Opp = Operators[Hmpo.data[key]['Op'][ii][2]]
# Check if phase is not identity
phase = Hmpo.data[key]['Op'][ii][2] != 'I'
if(phase):
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
sc = hp[jj] * Hmpo.data[key].yvals[ii][dist - 1]
perm = list(range(dist + 1, jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(jj + dist + 1, ll + 1))
tmp = np.tensordot(sc * Opl, self.rho, ([1], [jj]))
for nn in range(dist - 1):
tmp = np.tensordot(Opp, tmp, ([1], [jj + 1 + nn]))
Phi += np.transpose(np.tensordot(Opr, tmp,
([1], [jj + dist])),
perm)
else:
# No phase - identities
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
sc = hp[jj] * Hmpo.data[key].yvals[ii][dist - 1]
perm = list(range(2, jj + 1)) + [1] \
+ list(range(jj + 2, jj + dist + 1)) + [0] \
+ list(range(jj + dist + 1, ll + 1))
tmp = np.tensordot(sc * Opl, self.rho, ([1], [jj]))
Phi += np.transpose(np.tensordot(Opr, tmp,
([1], [jj + dist])),
perm)
# Loop over MBString terms
for ii in range(len(Hmpo.data['MBString'])):
hp = Hmpo.data['MBString'].hparam(ii, param, quenchtime)
Op = []
for opstr in Hmpo.data['MBString']['Op'][ii]:
Op.append(Operators[opstr])
nk = len(Hmpo.data['MBString']['Op'][ii])
for jj in range(ll - nk + 1):
# Last site with operator
kk = jj + nk - 1
# Adapt from exponential rule
dist = kk - jj
perm = list(range(dist + 1, jj + dist + 1)) + \
list(range(dist + 1))[::-1] + \
list(range(jj + dist + 1, ll + 1))
tmp = np.tensordot(hp[jj] * Op[0], self.rho,
([1], [jj]))
for kk in range(1, nk - 1):
tmp = np.tensordot(Op[kk], tmp, ([1], [jj + kk]))
Phi += np.transpose(np.tensordot(Op[-1], tmp,
([1], [jj + nk - 1])),
perm)
#if(Hmpo.data['lind1']):
# raise NotImplementedError('lind1 terms not implemented for ED.')
self.rho = np.reshape(self.rho, (dim, dim))
return np.reshape(Phi, (dim, dim))
[docs] def _dot_hmpo_mat_symm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a MPO which can represent the Hamiltonian or a measurement
with the density matrix. For arguments look into documention of
:py:func:`exactdiag.DensityMatrix.dot_hmpo_mat`. This is the version
for density matrices with symmetry adapted basis.
"""
spdim = len(self.SymmSec)
dtype = self.rho.dtype
Phi = np.zeros((spdim, spdim), dtype=self.rho.dtype)
ld = Operators['I'].shape[0]
ll = param['L']
keylist = ['site', 'bond', 'exponential', 'MBString', 'InfiniteFunction']
for key in keylist:
if(len(Hmpo.data[key]) == 0): continue
for k1, k2, elem in Hmpo.data[key].build_symm(param, self.SymmSec,
quenchtime, ll,
Hmpo.pbc):
Phi[k1, :] += elem * self.rho[k2, :]
if(len(Hmpo.data['product']) > 0):
raise NotImplementedError('Product terms not implemented for ED.')
if(len(Hmpo.data['FiniteFunction']) > 0):
raise NotImplementedError('FiniteFunction terms not implemented' +
' for ED.')
#if(Hmpo.data['lind1']):
# raise NotImplementedError('lind1 terms not implemented for ED.')
return Phi
[docs] def dot_hmpo_vec(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a Liouville operator in the Liouville space represented
as MPO / MPO rule set with a density matrix represented as superket.
The field ``self.rho`` has to be converted to the superket before
calling the subroutine. Result is returned as superket as well.
**Arguments**
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class or MPO representing
a measurement.
Operators : instance of :py:class:`Ops.OperatorList` (or dictionary)
Operators defined used in the simulation.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
maxdim : int
maximal dimension of the Hilbert space. Only referenced in ...
"""
if(self.SymmSec is None):
return self._dot_hmpo_vec_nosymm(Hmpo, Operators, param,
quenchtime, maxdim)
return self._dot_hmpo_vec_symm(Hmpo, Operators, param,
quenchtime, maxdim)
[docs] def _dot_hmpo_vec_nosymm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a MPO which represents the Liouville operator with the density
matrix represented as superket vector. For arguments look into the
documentation of :py:func:`exactdiag.DensityMatrix.dot_hmpo_vec`. This
is the version for density matrices without symmetries. Intended for
internal use.
"""
# Get shortcuts to system size and local dimension, dimension, shape
ll = param['L']
ld = Operators['I'].shape[0]
dim = ld**(2 * ll)
shape = [ld] * (2 * ll)
# Reshape superket into tensor, get empty array (always complex)
if(len(self.rho.shape) > 1): raise ValueError('Not a superket!')
self.rho = np.reshape(self.rho, shape)
phi = np.zeros((shape), dtype=np.complex128)
# Loop over local terms
for ii in range(len(Hmpo.data['site'])):
hp = Hmpo.data['site'].hparam(ii, param, quenchtime)
Op = Operators[Hmpo.data['site']['Op'][ii]]
# Terms of - 1j H x Id
for jj in range(ll):
perm = list(range(1, jj + 1)) + [0] \
+ list(range(jj + 1, 2 * ll))
phi -= np.transpose(np.tensordot(1j * hp[jj] * Op,
self.rho, ([1], [jj])), perm)
# Terms of 1j Id * H^T
for jj in range(ll):
perm = list(range(1, ll + jj + 1)) + [0] \
+ list(range(ll + jj + 1, 2 * ll))
phi += np.transpose(np.tensordot(1j * hp[jj] \
* Op.transpose(),
self.rho, ([1], [ll + jj])),
perm)
# Loop over bond terms
for ii in range(len(Hmpo.data['bond'])):
hp = Hmpo.data['bond'].hparam(ii, param, quenchtime)
Op0 = Operators[Hmpo.data['bond']['Op'][ii][0]]
Op1 = Operators[Hmpo.data['bond']['Op'][ii][1]]
# Terms of -1j H x Id
for jj in range(ll - 1):
perm = list(range(2, jj + 2)) + [1, 0] \
+ list(range(jj + 2, 2 * ll))
tmp = np.tensordot(1j * hp[jj] * Op0, self.rho,
([1], [jj]))
phi -= np.transpose(np.tensordot(Op1, tmp, ([1], [jj + 1])),
perm)
if(Hmpo.pbc):
perm = [0] + list(range(2, ll)) + [1] \
+ list(range(ll + 1, 2 * ll))
tmp = np.tensordot(1j * hp[ll - 1] * Op0, self.rho,
([1], [ll - 1]))
phi -= np.transpose(np.tensordot(Op1, tmp, ([1], [1])),
perm)
# Terms of 1j Id x H^T
for jj in range(ll - 1):
perm = list(range(2, ll + jj + 2)) + [1, 0] \
+ list(range(ll + jj + 2, 2 * ll))
tmp = np.tensordot(1j * hp[jj] * Op0.transpose(),
self.rho, ([1], [ll + jj]))
phi += np.transpose(np.tensordot(Op1.transpose(), tmp,
([1], [ll + jj + 1])), perm)
if(Hmpo.pbc):
perm = list(range(2, ll + 2)) + [0] \
+ list(range(ll + 2, 2 * ll)) + [1]
tmp = np.tensordot(1j * hp[ll - 1] *
Op0.conj().transpose(), self.rho,
([1], [2 * ll - 1]))
phi += np.transpose(np.tensordot(Op1.transpose(),
tmp, ([1], [ll + 1])),
perm)
# Loop over exponential terms
for ii in range(len(Hmpo.data['exponential'])):
hp = Hmpo.data['exponential'].hparam(ii, param, quenchtime)
Op0 = Operators[Hmpo.data['exponential']['Op'][ii][0]]
Op1 = Operators[Hmpo.data['exponential']['Op'][ii][1]]
Opp = Operators[Hmpo.data['exponential']['Op'][ii][2]]
dp = Hmpo.data['exponential']['dp'][ii]
# Check if phase is not identity
phase = (Hmpo.data['exponential']['Op'][ii][2] != 'I')
if(phase):
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
# Terms of -1j H x Id
perm = list(range(dist + 1, jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * hp[jj] * Op0, self.rho,
([1], [jj]))
for nn in range(dist - 1):
tmp = np.tensordot(dp * Opp, tmp,
([1], [jj + 1 + nn]))
phi -= np.transpose(np.tensordot(Op1, tmp,
([1], [jj + dist])),
perm)
# Terms of 1j I x H^T
perm = list(range(dist + 1, ll + jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(ll + jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * hp[jj]
* Op0.transpose(), self.rho,
([1], [ll + jj]))
for nn in range(dist - 1):
tmp = np.tensordot(dp * Opp.transpose(),
tmp, ([1], [ll + jj + 1 + nn]))
phi += np.transpose(np.tensordot(Op1.transpose(),
tmp,
([1], [ll + jj + dist])),
perm)
else:
# No phase - identities
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
# Terms of -1j H x Id
perm = list(range(2, jj + 2)) + [1] \
+ list(range(jj + 2, jj + dist + 1)) + [0] \
+ list(range(jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * hp[jj] * Op0, self.rho,
([1], [jj]))
tmp = np.tensordot(dp**(dist - 1) * Op1, tmp,
([1], [jj + dist]))
phi -= np.transpose(tmp, perm)
# Terms of 1j Id x H^t
perm = list(range(2, ll + jj + 2)) + [1] \
+ list(range(ll + jj + 2, ll + jj + dist + 1)) \
+ [0] + list(range(ll + jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * hp[jj]
* Op0.transpose(), self.rho,
([1], [ll + jj]))
tmp = np.tensordot(dp**(dist - 1)
* Op1.transpose(),
([1], [ll + jj + dist]))
phi += np.transpose(tmp, perm)
# Loop over infinite terms
key = 'InfiniteFunction'
for ii in range(len(Hmpo.data[key])):
hp = Hmpo.data[key].hparam(ii, param, quenchtime)
Opl = Operators[Hmpo.data[key]['Op'][ii][0]]
Opr = Operators[Hmpo.data[key]['Op'][ii][1]]
Opp = Operators[Hmpo.data[key]['Op'][ii][2]]
yvals = Hmpo.data[key].yvals[ii]
# Check if phase is not identity
phase = Hmpo.data[key]['Op'][ii][2] != 'I'
if(phase):
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
sc = hp[jj] * yvals[dist - 1]
# Terms of -1j H x Id
perm = list(range(dist + 1, jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * sc * Opl, self.rho,
([1], [jj]))
for nn in range(dist - 1):
tmp = np.tensordot(Opp, tmp,
([1], [jj + 1 + nn]))
phi -= np.transpose(np.tensordot(Opr, tmp,
([1], [jj + dist])),
perm)
# Terms of 1j I x H^T
perm = list(range(dist + 1, ll + jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(ll + jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * sc * Opl.transpose(),
self.rho, ([1], [ll + jj]))
for nn in range(dist - 1):
tmp = np.tensordot(Opp.transpose(),
tmp, ([1], [ll + jj + 1 + nn]))
phi += np.transpose(np.tensordot(Opr.transpose(), tmp,
([1], [ll + jj + dist])),
perm)
else:
# No phase - identities
for jj in range(ll - 1):
for kk in range(jj + 1, ll):
dist = kk - jj
sc = hp[jj] * yvals[dist - 1]
# Terms of -1j H x Id
perm = list(range(2, jj + 2)) + [1] \
+ list(range(jj + 2, jj + dist + 1)) + [0] \
+ list(range(jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * sc * Opl, self.rho,
([1], [jj]))
phi -= np.transpose(np.tensordot(Opr, tmp,
([1], [jj + dist])),
perm)
# Terms of 1j Id x H^t
perm = list(range(2, ll + jj + 2)) + [1] \
+ list(range(ll + jj + 2, ll + jj + dist + 1)) \
+ [0] + list(range(ll + jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * sc * Opl.transpose(),
self.rho, ([1], [ll + jj]))
phi += np.transpose(np.tensordot(Opr.transpose(), tmp,
([1], [ll + jj + dist])),
perm)
# Loop over MB String terms
for ii in range(len(Hmpo.data['MBString'])):
hp = Hmpo.data['exponential'].hparam(ii, param, quenchtime)
Op = []
for opstr in Hmpo.data['MBString']['Op'][ii]:
Op.append(Operators[opstr])
nk = len(Hmpo.data['MBString']['Op'][ii])
for jj in range(ll - nk + 1):
dist = nk - 1
# Terms of -1j * H x Id (adapt from exponential rule)
perm = list(range(dist + 1, jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * hp[jj] * Op[0], self.rho,
([1], [jj]))
for kk in range(1, nk - 1):
tmp = np.tensordot(Op[kk], tmp, ([1], [jj + kk]))
phi -= np.transpose(np.tensordot(Op[-1], tmp,
([1], [jj + nk - 1])), perm)
# Terms of 1j Id x H^T (adapt from exponential rule)
perm = list(range(dist + 1, ll + jj + dist + 1)) \
+ list(range(dist + 1))[::-1] \
+ list(range(ll + jj + dist + 1, 2 * ll))
tmp = np.tensordot(1j * hp[jj] * Op[0].transpose(),
self.rho, ([1], [ll + jj]))
for kk in range(1, nk - 1):
tmp = np.tensordot(Op[kk].transpose(), tmp,
([1], [ll + jj + kk]))
phi += np.transpose(np.tensordot(Op[-1].transpose(), tmp,
([1], [ll + jj + nk - 1])),
perm)
if(len(Hmpo.data['product']) > 0):
raise NotImplementedError('Product terms not implemented for ED.')
if(len(Hmpo.data['FiniteFunction']) > 0):
raise NotImplementedError('FiniteFunction terms not implemented'
+ ' for ED.')
for ii in range(len(Hmpo.data['lind1'])):
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
Op = Operators[Hmpo.data['lind1']['Op'][ii]]
LdL = Op.conj().transpose().dot(Op)
for jj in range(ll):
# Terms L x (Ld)^T
perm = list(range(2, jj + 2)) + [1] \
+ list(range(jj + 2, ll + jj + 1)) + [0] \
+ list(range(ll + jj + 1, 2 * ll))
tmp = np.tensordot(hp[jj] * Op, self.rho, ([1], [jj]))
phi += np.transpose(np.tensordot(Op.conj(), tmp,
([1], [ll + jj])), perm)
# Terms -0.5 Ld L x Id
perm = list(range(1, jj + 1)) + [0] \
+ list(range(jj + 1, 2 * ll))
phi -= np.transpose(np.tensordot(0.5 * hp[jj] * LdL,
self.rho, ([1], [jj])), perm)
# Terms -0.5 Id x (Ld L)^T
perm = list(range(1, ll + jj + 1)) + [0] \
+ list(range(ll + jj + 1, 2 * ll))
phi -= np.transpose(np.tensordot(0.5 * hp[jj] * LdL.transpose(),
self.rho, ([1], [ll + jj])),
perm)
self.rho = np.reshape(self.rho, (self.D**2))
return np.reshape(phi, (self.D**2))
[docs] def _dot_hmpo_vec_symm(self, Hmpo, Operators, param, quenchtime, maxdim):
"""
Multiply a MPO which represents the Liouville operator with the density
matrix represented as superket vector. For arguments look into the
documentation of :py:func:`exactdiag.DensityMatrix.dot_hmpo_vec`. This
is the version for density matrices with symmetry adapted basis.
Intended for internal use.
"""
# Reshape density matrix into superket
if(len(self.rho.shape) > 1): raise ValueError('Not a superket!')
phi = np.zeros((self.D**2), dtype=np.complex128)
dim = len(self.SymmSec)
ld = Operators['I'].shape[0]
ll = param['L']
# Hamiltonian parts
keylist = ['site', 'bond', 'exponential', 'MBString', 'InfiniteFunction']
for key in keylist:
if(len(Hmpo.data[key]) == 0): continue
for k1, k2, elem in Hmpo.data[key].build_symm(param, self.SymmSec,
quenchtime, ll,
Hmpo.pbc):
for xx in range(dim):
# Term - 1j H x Id
phi[k1 * dim + xx] -= 1j * elem * self.rho[k2 * dim + xx]
# Term 1j Id x H^T (swap k1 <--> k2)
phi[xx * dim + k2] += 1j * elem * self.rho[xx * dim + k1]
keylist = ['lind1']
for key in keylist:
if(len(Hmpo.data[key]) == 0): continue
for k1, k2, elem in Hmpo.data[key].lbuild_symm(param, self.SymmSec,
quenchtime, ll,
Hmpo.pbc):
phi[k1] += elem * self.rho[k2]
return phi
[docs] def get_bond_entropies(self):
"""
Function to return dummy value for bond entropy. There is no bond
entropy defined right now since the eigenvalues do not hold the
relation any more.
**Details**
We could split via SVD for each bipartition.
"""
for ii in range(self.ll + 1):
yield 0.0
# ----------------------------- PROPAGATION --------------------------------
[docs] def propagate_expm(self, dt, Hmpo, param, sparse, maxdim, quenchtime, qt):
"""
Propagate a state in the Liouville space according to the Lindblad
master equation.
**Arguments**
dt : float
time step for evolution.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class or MPO representing
a measurement.
param : dictionary
dictionary representing a MPS simulation.
sparse : boolean
flag if sparse matrix should be used for building Hamiltonian.
maxdim : int
maximal dimension of the Hilbert space. Only referenced in ...
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
qt : boolean
If ``True``, evolution should be carried out with quantum
trajectories which is impossible for density matrices and
therefore an exception is raised. If ``False``, evolution
in Liouville space is ok.
"""
if(qt): raise Exception("Quantum trajectories in Liouville space!")
if(self.lastprop is None):
expm = spla.expm if(sparse) else sla.expm
self.lastprop = Hmpo.build_liouville(param, self.SymmSec, sparse,
maxdim, quenchtime)
self.lastprop = expm(self.lastprop * dt)
self.rho = np.reshape(self.rho, (self.D**2))
self.rho = self.lastprop.dot(self.rho)
self.rho = np.reshape(self.rho, (self.D, self.D))
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
[docs] def propagate_krylov(self, dt, Hmpo, Operators, param, sparse, maxdim,
quenchtime, qt):
"""
Propagation with Krylov using the Arnoldi algorithm.
**Arguments**
dt : float
time step for evolution.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class or MPO representing
a measurement.
Operators : instance of :py:class:`Ops.OperatorList`
contains the operators defined on the local Hilbert space.
param : dictionary
dictionary representing a MPS simulation.
sparse : boolean
Build Liouville operators as sparse or dense matrix. Only
accessed in edlibmode=1.
maxdim : int
maximal dimension of the Hilbert space. Only referenced in ...
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
qt : boolean
If ``True``, evolution should be carried out with quantum
trajectories which is impossible for density matrices and
therefore an exception is raised. If ``False``, evolution
in Liouville space is ok.
**Details**
The following modes are available for the Krylov time evolution,
which is set in the convergence criteria for the quenches, i.e.,
:py:class:`convergence.KrylovConvParam` with the variable
``edlibmode``:
* 0 : Built Hamiltonian as matrix and use built-in scipy method.
* 1 : Built Hamiltonian as matrix and use EDLib Krylov.
* 2 : Use dot product of MPO and state vector and EDLib Krylov.
* 3 : As 2, but save Krylov vectors to hard disk to reduce RAM memory
usage.
* -1 : default. For :math`dim < 2^{10}` take edlibmode=0, then
for :math:`dim < 2^{15}` use edlibmode=2, otherwise edlibmode=3
Since the Liouville operator is in general not hermitian, Krylov uses
the Arnoldi implementation of the algorithm.
"""
if(qt): raise Exception("Quantum trajectories in Liouville space!")
# Reshape rho into vector
self.rho = np.reshape(self.rho, (self.D**2))
ConvParam = param['Quenches'].data[quenchtime[0]]
ConvParam = ConvParam['ConvergenceParameters'].data[0]
mode = ConvParam['edlibmode']
if(mode < 0):
if(self.rho.shape[0] < 2**10):
# Complete Liou, built-in scipy
mode = 0
elif(self.rho.shape[0] < 2**10):
# Complete Liou, EDLib Krylov (not used as default)
mode = 1
elif(self.rho.shape[0] < 2**15):
# MPO, EDLib Krylov
mode = 2
else:
# MPO, EDLib Krylov, write Krylov vectors to hard drive
mode = 3
if(mode == 0):
# Complete H, built-in scipy
if(self.lastprop is None):
self.lastprop = Hmpo.build_liouville(param, self.SymmSec, True,
maxdim, quenchtime)
self.save_propagator_spectrum(self.lastprop, param, quenchtime)
self.rho = spla.expm_multiply(dt * self.lastprop, self.rho)
# Reshape rho from vector back to matrix
self.rho = np.reshape(self.rho, (self.D, self.D))
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
nvec = ConvParam['max_num_lanczos_iter']
tol = ConvParam['lanczos_tol']
numzero = 1e-16
hmat = np.zeros((nvec, nvec), dtype=np.complex128)
norm = np.sqrt(np.real(np.sum(np.conj(self.rho) * self.rho)))
beta = deepcopy(norm)
# Set nn in case if-break conditions are never reached
nn = nvec
if((mode == 1) or (mode == 2)):
# EDLib Krylov (Complete H, MPO dot product)
# ------------
base = np.zeros((nvec, self.rho.shape[0]), dtype=np.complex128)
if((mode == 1) and (self.lastprop is None)):
self.lastprop = Hmpo.build_liouville(param, self.SymmSec,
sparse, maxdim,
quenchtime)
for jj in range(nvec):
self.rho = self.rho / norm
base[jj, :] = self.rho
if(mode == 1):
self.rho = self.lastprop.dot(self.rho)
else:
self.rho = self.dot_hmpo_vec(Hmpo, Operators, param,
quenchtime, maxdim)
for ii in range(jj):
hmat[ii, jj] = np.sum(np.conj(self.rho) * base[ii, :])
self.rho -= hmat[ii, jj] * base[ii, :]
norm = np.sqrt(np.real(np.sum(np.conj(self.rho) * self.rho)))
hmat[jj + 1, jj] = norm
prop = sla.expm(dt * hmat[:jj + 1, :jj + 1])
est = 2.0 * norm * np.abs(prop[jj, 0])
if(est < tol):
nn = jj
break
if(norm < 10 * numzero):
nn = jj
break
self.rho = beta * prop[0, 0] * base[0, :]
for jj in range(1, nn):
self.rho += beta * prop[jj, 0] * base[jj, :]
else:
# MPO, EDLib Krylov, write Krylov vectors to hard drive
# -----------------------------------------------------
# Save vectors to hard disk to save memory
fbase = param['Write_Directory'] + param['job_ID'] + \
param['unique_ID'] + 'Krylovbase_'
for jj in range(nvec):
self.rho = self.rho / norm
np.save(fbase + str(jj) + '.npy', self.rho)
self.rho = self.dot_hmpo_vec(Hmpo, Operators, param,
quenchtime, maxdim)
for ii in range(jj):
Lrho = np.load(fbase + str(ii) + '.npy')
hmat[ii, jj] = np.sum(np.conj(self.rho) * Lrho)
self.rho -= hmat[ii, jj] * Lrho
norm = np.sqrt(np.real(np.sum(np.conj(self.rho) * self.rho)))
hmat[jj + 1, jj] = norm
prop = sla.expm(dt * hmat[:jj + 1, :jj + 1])
est = 2.0 * norm * np.abs(prop[jj, 0])
if(est < tol):
nn = jj
del Lrho
break
if(norm < 10 * numzero):
nn = jj
del Lrho
break
self.rho = np.load(fbase + str(0) + '.npy')
self.rho = beta * prop[0, 0] * self.rho
for jj in range(1, nn):
tmp = np.load(fbase + str(jj) + '.npy')
remove_file(fbase + str(jj) + '.npy')
self.rho += beta * prop[jj, 0] * tmp
# Reshape rho from vector back to matrix
self.rho = np.reshape(self.rho, (self.D, self.D))
if(param['Quenches'].is_time_dependent(quenchtime[0])):
self.lastprop = None
return
[docs] def propagate_gates(self, *args, **kwargs):
"""
Propagation with gates is not implemented for density matrices
"""
raise NotImplementedError("Gates with density matrices")
# ------------------- Trotter decomposition functions ----------------------
[docs] def _trotter_initialize(self, qt):
"""
Density matrices without symmetry-adapted basis are converted to
a rank-2L tensor.
**Arguments**
qt : boolean
If ``True``, an error is raised because quantum trajectories
on a density matrix are not supported. If ``False``, then ok.
"""
if(qt): raise Exception("Quantum trajectories in Liouville space!")
if(self.SymmSec is None):
# Only density matrices without symmetry are reshaped
self.rho = np.reshape(self.rho, [self.ld] * 2 * self.ll)
return
[docs] def _trotter_finalize(self, *args):
"""
Density matrices without symmetry-adapted basis are converted
from the rank-2L tensor back to their matrix representation.
"""
if(self.SymmSec is None):
# Only density matrices without symmetry were reshaped
self.rho = np.reshape(self.rho, [self.ld**self.ll] * 2)
return
[docs] def trotter_layer(self, dt, Hmpo, Operators, param, quenchtime,
even, *args):
"""
Apply two site operators for the Trotter decomposition for the odd or
even layer.
**Arguments**
dt : float
time step for evolution.
Hmpo : instance of :py:class:`MPO.MPO`
Hamiltonian represented in the MPO class or MPO representing
a measurement.
Operators : instance of :py:class:`Ops.OperatorList`
contains the operators for building the Hamiltonian and
observables.
param : dictionary
dictionary representing a MPS simulation.
quenchtime : tuple of int and float
Storing the number of the quench and the actual time
even : boolean
If ``True``, even layer of operators is applied on sites
(2k , 2k + 1). If ``False``, we use the odd layer.
"""
ld = Operators['I'].shape[0]
ll = param['L']
if((self.bidx is None) and (not self.SymmSec is None)):
self.bidx = Lidx_symmetry_blocks(Operators, param, self.ld)
for elem in self.trotter_Ham_2site(Hmpo, Operators, param, quenchtime,
even):
jj = elem[1]
if((jj, dt) in self.lastprop):
Op2 = self.lastprop[(jj, dt)]
else:
Liou = 1j * (np.kron(-elem[0], np.eye(elem[0].shape[0])) \
+ np.kron(np.eye(elem[0].shape[0]),
elem[0].transpose()))
jjp1 = (jj + 1)%ll
for ii in range(len(Hmpo.data['lind1']['hp'])):
hp = Hmpo.data['lind1'].hparam(ii, param, quenchtime)
if(Hmpo.pbc):
hp *= 0.5
else:
hp[1:-1] = 0.5 * np.array(hp[1:-1])
Op = Operators[Hmpo.data['lind1']['Op'][ii]]
# Terms L x (L^dagger)^T
Liou += np.kron(np.kron(hp[jj] * Op, np.eye(ld)),
np.kron(Op.conj(), np.eye(ld)))
Liou += np.kron(np.kron(hp[jjp1] * np.eye(ld), Op),
np.kron(np.eye(ld), Op.conj()))
# Terms L^dagger L x Id
Tmp = -0.5 * np.dot(Op.transpose().conj(), Op)
Liou += np.kron(hp[jj] * Tmp, np.eye(ld**3))
Liou += np.kron(np.kron(hp[jjp1] * np.eye(ld), Tmp),
np.eye(ld**2))
# Terms Id x (L^dagger L)^T
Tmp = Tmp.transpose()
Liou += np.kron(np.kron(np.eye(ld**2), hp[jj] * Tmp),
np.eye(ld))
Liou += np.kron(np.eye(ld**3), hp[jjp1] * Tmp)
Op2 = self._matrix_exp_trotter(dt * Liou)
self.lastprop[(jj, dt)] = Op2
self.apply_op2(Op2, jj, elem[2])
return
[docs] def _matrix_exp_trotter(self, Liou):
"""
Taking the matrix exponential of the Liouville operator to obtain the
propagator. Without symmetries, the propagator is returned as matrix.
With symmetries, a dictionary is returned where the key is the row
index and the value are the four local indices plus the matrix element.
In the latter case the blockdiagonal structure is used.
**Arguments**
Liou : matrix
a matrix representing the Liouville operator on two sites in the
Liouville space.
**Details**
Functions accesses class variabes ``self.bidx``, a dictionary of tuples
or ``None``. If ``None``, no symmetries are used. If symmetries are
used, bidx is a dictionary where each set of tuples contains the indices
for the same symmetry sector and the corresponding site indices.
"""
if(self.bidx is None):
# No symmetries
return sla.expm(Liou)
# Symmetries using bidx
dic = {}
for key in self.bidx.keys():
Sub = sla.expm(Liou[key, :][:, key])
for i1 in range(len(key)):
j1 = self.bidx[key][i1]
for i2 in range(len(key)):
j2 = self.bidx[key][i2]
if(Sub[i1, i2] == 0.0): continue
try:
dic[(j1)].append((j2[0], j2[1], j2[2], j2[3],
Sub[i1, i2]))
except KeyError:
dic[(j1)] = [(j2[0], j2[1], j2[2], j2[3],
Sub[i1, i2])]
return dic
[docs] def _apply_op2_nosymm(self, Op, ii, ll_1_term, *args):
"""
Apply a two-site operator to the density matrix on a specified site.
This is the version for density matrices without a symmetry adapted
basis.
**Arguments**
Op : np 2d array
The operator acting on two sites represented as d^4 x d^4 matrix.
ii : int
Apply term to the sites ii and (ii + 1).
ll_1_term : bool, optional
Boolean to indicate the term ii=L for periodic boundary
conditions.
"""
shape = [self.ld] * 8
Op = np.reshape(Op, shape)
if(ll_1_term):
perm = [1] + list(range(4, self.ll + 2)) + [0, 3] \
+ list(range(self.ll + 2, 2 * self.ll)) + [2]
cs = [self.ll, 0, 2 * self.ll, self.ll + 1]
else:
perm = list(range(4, ii + 4)) + [0, 1] \
+ list(range(ii + 4, self.ll + ii + 2)) + [2, 3] \
+ list(range(self.ll + ii + 2, 2 * self.ll))
cs = [ii, ii + 1, ii + self.ll, ii + 1 + self.ll]
self.rho = np.tensordot(Op, self.rho, ([4, 5, 6, 7], cs))
self.rho = np.transpose(self.rho, perm)
return
[docs] def _apply_op2_symm(self, Op, ii, ll_1_term, *args):
"""
Apply a two-site operator to the density matrix on a specified site.
This is the version for density matrices with a symmetry adapted
basis.
**Arguments**
Op : np 2d array
The operator acting on two sites represented as d^4 x d^4 matrix.
ii : int
Apply term to the sites ii and (ii + 1).
ll_1_term : bool, optional
Boolean to indicate the term ii=L for periodic boundary
conditions.
"""
Rho = np.zeros(self.rho.shape, dtype=np.complex128)
if(ll_1_term):
c1 = 0
c2 = self.ll - 1
else:
c1 = ii
c2 = ii + 1
for i1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(i1)
idx2 = deepcopy(idx1)
for j1 in range(len(self.SymmSec)):
jdx1 = self.SymmSec(j1)
jdx2 = deepcopy(jdx1)
try:
elems = Op[(idx1[c1], idx1[c2], jdx1[c1], jdx1[c2])]
except KeyError:
continue
for il, ir, jl, jr, elem in elems:
idx2[c1] = il
idx2[c2] = ir
jdx2[c1] = jl
jdx2[c2] = jr
try:
i2 = self.SymmSec[tuple(idx2)]
j2 = self.SymmSec[tuple(jdx2)]
except KeyError:
continue
Rho[i1, j1] += elem * self.rho[i2, j2]
self.rho = Rho
return
[docs] def _apply_op_nosymm(self, Op, idx):
"""
Apply an operator/gate on an arbitrary number of sites. The
gate is defined on Hilbert space, not Liouville space.
**Arguments**
Op : np 2d array
The operator acting on the site represented as matrix.
idx : list of integers
The sited to be contracted with the operator.
"""
# Number of sites acted on
nn = len(idx)
# Get shape and transform to tensor
shape = tuple([self.ld] * 2 * nn)
Op_ = np.reshape(Op, shape)
# Apply Gate x Rho
# ----------------
# Get list of contraction indices for Op
contr = list(range(nn, 2 * nn))
self.Rho = np.tensordot(Op_, self.Rho, (contr, idx))
perm = []
pos = nn
for ii in range(self.ll):
if(ii in idx):
for jj in range(nn):
if(ii == idx[jj]):
perm.append(jj)
break
else:
perm.append(pos)
pos += 1
self.Rho = np.transpose(self.Rho, perm)
# Apply Rho x Gate^Dagger
# -----------------------
idx_ = list(np.array(idx, dtype=np.int32) + self.ll)
contr = list(range(nn))
Op_ = np.reshape(np.transpose(np.conj(Op)), shape)
self.Rho = np.tensordot(self.Rho, Op_, (idx_, contr))
perm = list(self.ll)
pos = self.ll
posc = 2 * self.ll - nn
for ii in range(self.ll):
if(ii in idx):
# Index at the back
for jj in range(nn):
if(ii == idx[jj]):
perm.append(posc + jj)
break
else:
perm.append(pos)
pos += 1
self.Rho = np.transpose(self.Rho, perm)
return
[docs] def _apply_op_symm(self, Op, idx):
"""
Apply an operator/gate on an arbitrary number of sites to a
symmetry conserving density matrix. The gate is defined
on Hilbert space, not Liouville space.
**Arguments**
Op : np 2d array
The operator acting on the site represented as matrix.
idx : list of integers
The sites to be contracted with the operator.
"""
# 1) Extract operator (pre-calculating this would certainly
# be more efficient than doing it every time)
dic = {}
for key in self.bidx.keys():
Sub = Op[key, :][:, key]
for i1 in range(len(key)):
j1 = self.bidx[key][i1]
for i2 in range(len(key)):
j2 = self.bidx[key][i2]
if(Sub[i1, i2] == 0.0): continue
val = tuple(list(j2) + [Sub[i1, i2]])
try:
dic[j1].append(val)
except KeyError:
dic[j1] = [val]
# 2) Multiplication Gate x Rho
Sigma = np.zeros(self.rho.shape, dtype=np.complex128)
for k1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(k1)
idx2 = deepcopy(idx1)
try:
elems = dic[tuple(idx1[idx])]
except KeyError:
continue
for entry in elems:
idx2[idx] = entry[:len(idx)]
elem = entry[-1]
try:
k2 = self.SymmSec[tuple(idx2)]
except KeyError:
continue
Sigma[k1, :] += elem * self.Psi[k2, :]
# 3) Multiplication Rho x Gate^dagger
self.Psi[:, :] = 0.0
for k1 in range(len(self.SymmSec)):
idx1 = self.SymmSec(k1)
idx2 = deepcopy(idx1)
try:
elems = dic[tuple(idx1[idx])]
except KeyError:
continue
for entry in elems:
idx2[idx] = entry[:len(idx)]
elem = entry[-1]
try:
k2 = self.SymmSec[tuple(idx2)]
except KeyError:
continue
self.Psi[:, k1] += np.conj(elem) * Sigma[:, k2]
return
[docs] def distance(self, State, dist):
"""
Measure the distance between the pure state and another quantum
state.
**Arguments**
State : 1d or 2d numpy array
Represents another pure state (1d) or a density matrix (2d).
distance : str
Define the type of distance measure. Either ``Trace`` or
``Infidelity``.
"""
return qi.distance(self.rho, State, dist)
[docs] def fidelity(self, Phi):
"""
Return the fidelity. Not implemented for any density matrix yet.
**Arguments**
Phi : pure state or density matrix
Calculate fidelity to this state. It can be represented as
1d or 2d numpy array, an instance of
:py:class:`exactdiag.PureState`, or as instance of
:py:class:`exactdiag.DensityMatrix`.
**Details**
Implementation for density matrices is missing to avoid calculating
the square root of large matrices.
"""
if(isinstance(Phi, PureState)):
return self._fidelity_with_state(Phi)
elif(isinstance(Phi, DensityMatrix)):
return 1 - self.distance(Phi.rho, 'Infidelity')
elif(len(Phi.shape) == 1):
return self._fidelity_with_state(Phi)
elif(len(Phi.shape) == 2):
return 1 - self.distance(Phi, 'Infidelity')
else:
raise Exception("Unknown type of quantum state.")
return 0.0
[docs] def _fidelity_with_state(self, Phi):
"""
Return the fidelity between the density matrix and a quantum
state. This is possible independent of the size of the density
matrix, since the fidelity simplifies to
:math:`F(\\rho, | \psi \\rangle) =`
:math:`\sqrt{\langle \psi | \\rho | \psi \\rangle}`.
**Arguments**
Phi : instance of :py:class:`exactdiag.PureState` or numpy 1D array
Calculate the fidelity to this state vector.
"""
if(isinstance(Phi, PureState)):
# Instance of PureState
return np.sqrt(np.dot(np.conj(Phi.Psi), self.rho.dot(Phi.Psi)))
# Must be numpy array
return np.sqrt(np.dot(np.conj(Phi), self.rho.dot(Phi)))
[docs] def trace(self):
"""
Return the trace of a density matrix.
"""
return np.trace(self.rho)
[docs] def norm(self):
"""
Return the trace as norm of a density matrix.
"""
return self.trace()
[docs] def purity(self):
"""
The purity of a pure quantum state is always 1.
"""
return np.real(np.trace(self.rho.dot(self.rho)))
[docs]def idx_symmetry_blocks(Operators, param, ll=2):
"""
Generates the indices for blocks conserving the symmetry.
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
contains all operators, amongst them the diagonal generators
of the symmetry.
param : dictionary
containing the string identifiers for the generators of the
U(1) and Z2 symmetries.
ll : int, optional
Number of sites to generate the blocks
Default to 2.
**Details**
If no symmetries are present, the funtion returns ``None``. Otherwise
the following output in form of a dictionary is generated. The keys
of the dictionary are a tuple of indices in the combined Hilbert space
of ll sites with the same quantum number. The value for the corresponding
key is a tupled tuple. The i-th tuple contains the local indices for the
i-th combined index.
"""
ld = Operators['I'].shape[0]
if('Abelian_generators' in param):
GenU1 = param['Abelian_generators']
else:
GenU1 = []
if('Discrete_generators' in param):
GenZ2 = param['Discrete_generators']
else:
GenZ2 = []
nn = len(GenU1) + len(GenZ2)
if(nn == 0):
# Fast return - no symmetries
return None
QN = np.zeros((ld**ll, nn))
ii = 0
# Store U1 quantum numbers
for key in GenU1:
Op = np.diag(Operators[key])
for jj in range(ll):
QN[:, ii] += np.kron(np.kron(np.ones(ld**jj), Op),
np.ones(ld**(ll - 1 - jj)))
ii += 1
# Store Z2 quantum numbers
for key in GenZ2:
Op = np.diag(Operators[key])
for jj in range(ll):
QN[:, ii] += np.kron(np.kron(np.ones(ld**jj), Op),
np.ones(ld**(ll - 1 - jj)))
QN[:, ii] = np.array(np.array(QN[:, ii], dtype=np.int32)%2, dtype=float)
ii += 1
# Dictionary with quantum numbers as key and list of basis indices
# corresponding to index of combined Hilbert space.
dici = {}
# Dictionary with quantum numbers as key and list of indices in local Hilbert
# space as value
dicij = {}
# Indices of local Hilbert space
idx = [0] * ll
for jj in range(QN.shape[0]):
key = tuple(np.round(QN[jj, :], decimals=5))
try:
Li = dici[key]
Li.append(jj)
dici[key] = Li
Li = dicij[key]
Li.append(tuple(idx))
dicij[key] = Li
except KeyError:
dici[key] = [jj]
dicij[key] = [tuple(idx)]
idx[-1] += 1
if(idx[-1] == ld):
for ii in range(2, ll + 1):
idx[ll - ii + 1] = 0
idx[ll - ii] += 1
if(idx[ll - ii] < ld): break
# Bundle lists of identic hashes in nested list
bidx = {}
for key in dici.keys():
bidx[tuple(dici[key])] = tuple(dicij[key])
return bidx
[docs]def Lidx_symmetry_blocks(Operators, param, ld, ll=2):
"""
Function analogous to ``idx_symmetry_block``, but constructing the
mapping for the Liouville space used for density matrices:
**Arguments**
Operators : instance of :py:class:`Ops.OperatorList`
contains all operators, amongst them the diagonal generators
of the symmetry.
param : dictionary
containing the string identifiers for the generators of the
U(1) and Z2 symmetries.
ld : integer
local dimension of the Hilbert space.
ll : int, optional
Number of sites to generate the blocks
Default to 2.
"""
# Get standard mapping and extend to Liouville space from there
bidx = idx_symmetry_blocks(Operators, param, ll=ll)
Lbidx = {}
for k1 in bidx.keys():
for k2 in bidx.keys():
tmpli = []
for elem1 in k1:
for elem2 in k2:
tmpli.append(elem1 * ld**ll + elem2)
tmpval = []
for kk1 in range(len(k1)):
for kk2 in range(len(k2)):
tmpval.append(tuple(list(bidx[k1][kk1])
+ list(bidx[k2][kk2])))
Lbidx[tuple(tmpli)] = tmpval
return Lbidx