"""Dynamics.py
Module to define time evolutions Currently under heavy development.
last revision 090111-MLW
See the manual sec.## for more information about the representation of a dynamical process in MPSPyLib."""
from math import ceil, sqrt, floor
import numpy as np
from scipy import linalg
#from scipy.optimize import leastsq
#from scipy.special import legendre
import random
import string
#import os.path
from os import makedirs
#from .convergence import KrylovConvParam, TEBDConvParam, TDVPConvParam
#from .convergence import LRKConvParam, MpdoplusConvParam, ExpmConvParam
[docs]def CheckOrDefault(p,parameter,default):
"""Set default value in the dictionary, if key not already set
**Arguments**
p : dictionary
dictionary to be checked and, if necessary, modified.
parameter : string, boolean, integer, float, tuple, ...
key in the dictionary
default : any
default value to be set for the key `parameter` if the key does not
yes exist in `p`
"""
if parameter not in p:
p[parameter]=default
[docs]def SetupGLQ(order,nexp=2):
"""Set up the parameters for Gauss-Legendre quadrature and its
Magnus variants. Returned is the 1D array x and the 2D array
g.
**Arguments**
order : integer, either 4 or 6
order of the Magnus variants??
nexp : integer, optional
number of exponentials
for order=2, nexp can either be 2 or 3.
for order=6, nexp must be 5.
"""
if order==2:
#points for 2nd order Gauss-Legendre quadrature
x=np.zeros(1)
x[0]=0.5
#weights for 2nd order Gauss-Legendre quadrature
w=np.zeros(1)
w[0]=1.0
sa=1; sb=1
f=np.zeros((sa,2))
f[0,0]=1.0;
g=np.zeros((sa,sb))
elif order==4:
if nexp==2:
#points for 3rd order Gauss-Legendre quadrature
x=np.zeros(2)
x[0]=0.5-sqrt(3.0/36.0); x[1]=0.5+sqrt(3.0/36.0)
#weights for 3rd order Gauss-Legendre quadrature
w=np.zeros(2)
w[0]=0.5; w[1]=0.5
sa=2; sb=2
f=np.zeros((sa,3))
f[0,0]=0.5; f[0,1]=1.0/3.0
f[1,0]=0.5; f[1,1]=-1.0/3.0
g=np.zeros((sa,sb))
else:
if nexp!=3:
raise Exception("nexp must be 2 or 3 for order=4 in SetupGLQ!")
#points for 3rd order Gauss-Legendre quadrature
x=np.zeros(3)
x[0]=0.5-sqrt(3.0/20.0); x[1]=0.5; x[2]=0.5+sqrt(3.0/20.0)
#weights for 3rd order Gauss-Legendre quadrature
w=np.zeros(3)
w[0]=5.0/18.0; w[1]=4.0/9.0; w[2]=5.0/18.0
sa=3; sb=3
f=np.zeros((sa,sb))
f[0,0]=11.0/40.0; f[0,1]=20.0/87.0; f[0,2]=7.0/50.0
f[1,0]=9.0/20.0; f[1,2]=-7.0/25.0
f[2,0]=11.0/40.0; f[2,1]=-20.0/87.0; f[2,2]=7.0/50.0
g=np.zeros((sa,sb))
elif order==6:
#points for 4th order Gauss-Legendre quadrature
x=np.zeros(4)
x[0]=0.5-sqrt((3.0+2.0*sqrt(6.0/5.0))/28.0); x[1]=0.5-sqrt((3.0-2.0*sqrt(6.0/5.0))/28.0);
x[2]=0.5+sqrt((3.0-2.0*sqrt(6.0/5.0))/28.0); x[3]=0.5+sqrt((3.0+2.0*sqrt(6.0/5.0))/28.0)
#weights for 4th order Gauss-Legendre quadrature
w=np.zeros(4)
w[0]=(18.0-sqrt(30.0))/72.0; w[1]=(18.0+sqrt(30.0))/72.0;
w[2]=(18.0+sqrt(30.0))/72.0; w[3]=(18.0-sqrt(30.0))/72.0
#Coefficients for CF6:5Opt
sa=5; sb=4
f=np.zeros((sa,5))
f[0,0]=0.1714; f[0,1]=0.15409059414309687213; f[0,2]=0.11947178242929061641; f[0,3]=0.07195
f[1,0]=0.37496374319946236513; f[1,1]=0.13813675394387646682;
f[1,2]=-0.13090674649282935743; f[1,3]=-0.21123356253315514306
f[2,0]=1.0-2.0*f[1,0]-2.0*f[0,0]; f[2,1]=0.0; f[2,2]=-2.0*f[1,2]-2.0*f[0,2]; f[2,3]=0.0
f[3,0]=0.37496374319946236513; f[3,1]=-0.13813675394387646682;
f[3,2]=-0.13090674649282935743; f[3,3]=0.21123356253315514306
f[4,0]=0.1714; f[4,1]=-0.15409059414309687213; f[4,2]=0.11947178242929061641; f[4,3]=-0.07195
g=np.zeros((sa,sb))
else:
raise Exception("order="+str(order)+" not recognized in SetupGLQ! Use 4 or 6!")
for i in range(sa):
for m in range(sb):
tmp=0.0
for n in range(1, order // 2 + 2):
tmp+=(2.0*n-1.0)*legendre(n-1)(2*x[m]-1)*f[i,n-1]
g[i,m]=w[m]*tmp
return x,g
[docs]class QuenchList:
"""Class defining a series of time-dependent changes of Hamiltonian
parameters, aka "Quenches".
Once an instance of the Quenchlist class has been created, quenches are
added using the AddQuench method. Quenches are communicated to the Fortran
backend by specifying a QuenchList object as the element of the
dynamicsParameters dictionary with key 'Quenches'
**Arguments**
H : instance of mps.MPO
the MPO defining the system (is equal for all through statics and all
dynamics).
**Variables**
data : list
list where the quenches are added later on. This is not an argument
in the constructor. The list contains one dictionary for each quench.
"""
def __init__(self, H):
self.H = H
self.data=[]
self.unique = str(id(self))
[docs] def __unicode__(self):
"""Return unicode string."""
return "QuenchList_" + str(self.unique)
[docs] def __str__(self):
"""Return string representation"""
return self.__unicode__()
def classname(self):
return 'QuenchList'
[docs] def AddQuench(self, Hparams, time, deltat, funcs,
ConvergenceParameters=None, stepsforoutput=1, timedep=None):
"""Add a quench process to the specification of the dynamics.
**Arguments**
Hparams : string or list of strings
containing the Hamiltonian parameters to be quenched. The
name of the Hamiltonians parameters are identified via strings.
time : float
total time for this quench
deltat : float
infinitesimal time step of the quench
funcs : callable or list of callables
functional forms describing the time-dependence of the Hamitlonian
parameters in Hparams. The number of functions must correspond to
the number of `Hparams`
ConvergenceParameters : instance of ConvParam for time evolultions
Specifies the time evolution method and its convergence
parameters.
* :py:class:`convergence.KrylovConvParam`: slow but reliable.
* :py:class:`convergence.TEBDConvParam` : TEBD routines based on
the Sornborger-Stewart decompositions (PRA 60, 1956 (1999)),
which only work with nearest-neighbor Hamiltonians. Either 2nd or
4th order.
* :py:class:`convergence.TDVPConvParam`: TDVP_2 is the
second-order time-dependent variational principle method of
arXiv:1408.5056, which *in principle* works with any
Hamiltonian.
* :py:class:`convergence.LRKConvParam` : not implemented are the
LRK_2 and LRK_4 are second and fourth-order local Runge-Kutta
methods (PRB 91, 165112 (2015)). LRK_2 is the default, as it
works with any Hamiltonian and is fairly reliable.
Default: instance of :py:class:`convergence.KrylovConvParam` is
generated during the function containing the default values.
stepsforoutput : integer, optional
allows for coarser measurement
default to 1
timedep : boolean, optional
This flag, right now only used for the EDLib, signals that
the couplings in the quench are time dependent (``True``),
are time independent (``False``) or the time dependence of the
coupling should be determined thourhg the length of the
``Hparams`` arguments (``None``). This flag might be used for
a sudden quench where parameters change in regard to previous
setting, but not over time.
Default to ``None``.
"""
Quench={}
if isinstance(Hparams,list):
lenHparams=len(Hparams)
else:
lenHparams=1
Hparams=[Hparams]
if isinstance(funcs,list):
lenfuncs=len(funcs)
else:
lenfuncs=1
funcs=[funcs]
if lenHparams!=lenfuncs:
raise Exception("The length of the Hparams and funcs arrays " +
"in AddQuench should be the same!")
if ConvergenceParameters is None:
# Choose default convergence parameters
ConvergenceParameters = KrylovConvParam()
else:
# Check if valid instance of class
if((not isinstance(ConvergenceParameters, KrylovConvParam))
and (not isinstance(ConvergenceParameters, TEBDConvParam))
and (not isinstance(ConvergenceParameters, TDVPConvParam))
and (not isinstance(ConvergenceParameters, LRKConvParam))
and (not isinstance(ConvergenceParameters, MpdoplusConvParam))
and (not isinstance(ConvergenceParameters, ExpmConvParam))):
raise ValueError('Convergence parameter class for time ' +
'evolution not known.')
# Check that the parameters in Hparams exist in the MPO template
hparams = self.H.getHparams()
for Hparam in Hparams:
if Hparam not in hparams:
raise Exception("Hamiltonian parameter " + str(Hparam) +
" passed in to AddQuench not found in " +
"MPO template!")
# Wishlist: determine if function is constant? If so, don't use
# magnus, just use ordinary exponential
Quench['Hparams']=Hparams
Quench['funcs']=funcs
Quench['stepsforoutput']=stepsforoutput
Quench['time']=time
Quench['deltat']=deltat
Quench['ConvergenceParameters']=ConvergenceParameters
Quench['tdep'] = timedep
self.data.append(Quench)
return
[docs] def is_time_dependent(self, ii):
"""
Check if a quench is time dependent. Returns ``True`` for time dependent
quenches and ``False`` for time independent quenches.
**Arguments**
ii : int
Check for quench ii if it is time dependent. Time dependent quenches
have at least one entry in the Hparams to be modified if flag
``timedep`` was not passed to quench.
"""
if(self.data[ii]['tdep'] is not None): return self.data[ii]['tdep']
return len(self.data[ii]['Hparams']) > 0
[docs] def write(self, filestub):
"""
Write the information about the quenches.
**Arguments**
filestub : string
filename where the information are written
"""
# # Filenames, create subdirectory
# writedir = p['Write_Directory']
# totalID = p['job_ID'] + p['unique_ID']
# subdir = p['Output_Directory'] + totalID + 'ObsOutDynamics'
# if(not os.path.exists(subdir)):
# makedirs(subdir)
# if(isinstance(p['Quenches'], QuantumCircuits)): return
# Get fundamental data
PostProcess = self.H.PostProcess
nquenches = len(self.data)
# p['nquenches'] = nquenches
if not PostProcess:
# files = []
MetaFileName = filestub + 'Dynamics.dat'
# p['int_files'].append(MetaFileName)
# files.append(MetaFileName)
MetaFile = open(MetaFileName, 'w')
# write out quenches
MetaFile.write('%16i'%(nquenches)+'\n')
#Write out what method is being used
for ii in range(nquenches):
mkey = self.data[ii]['ConvergenceParameters'].get_method()
MetaFile.write('%16i'%(mkey) + '\n')
MetaFile.close()
# if('QT' in p): qtcopy(MetaFileName, p['QT'])
#if(not PostProcess):
# ...
self.H.IndexHParams()
thehparams = self.H.getHparams()
# p['nmeasurements'] = []
time = 0.0
for ii in range(nquenches):
# Order is contained in convergence parameters
myTO = self.data[ii]['ConvergenceParameters'].get_TO()
splitTO = myTO.split('_')
# Commutator-free magnus expansion setup
# --------------------------
order = int(splitTO[0])
if(len(splitTO) > 1):
nexp = int(splitTO[1])
x, g = SetupGLQ(order, nexp)
else:
x, g = SetupGLQ(order)
if(order == 2):
nexp = 1
elif(order == 6):
nexp = 5
# Set up GL parameters
sa = g.shape[0]
sb = g.shape[1]
scalestr = ''
# Get the "scale factors" s_i=\sum_{m}g_{i,m} to scale the time
# step by. This ensures the proper behavior of all terms which
# remain constant. For those which are not constant, divide
# through by the scale factor
scalefac = np.zeros(sa)
for k in range(sa):
for m in range(sb):
scalefac[k] += g[k,m]
scalestr += '%30.15E'%(scalefac[k])
myHparams = self.data[ii]['Hparams']
ftmp = filestub + 'tMPS_' + str(ii + 1) + '.dat'
self.data[ii]['ConvergenceParameters'].write(ftmp)
# if(not PostProcess and ('QT' in p)): qtcopy(ftmp, p['QT'])
nmeasurements = 0
nsteps = int(floor(self.data[ii]['time'] / self.data[ii]['deltat']))
extradt = self.data[ii]['time'] - nsteps * self.data[ii]['deltat']
hasextradt = (extradt > 0.000001)
# Number of parameters being changed
nchanged = len(myHparams)
if not PostProcess:
ThisQuenchname = filestub + 'Quench_' + str(ii + 1) + '.dat'
# p['int_files'].append(ThisQuenchname)
# files.append(ThisQuenchname)
ThisQuench = open(ThisQuenchname, 'w')
ThisQuench.write('%16i'%(nsteps)
+ '%16i'%(self.data[ii]['stepsforoutput'])
+ '%30.15E\n'%(self.data[ii]['deltat']))
ThisQuench.write('%16i\n'%(nexp))
if(hasextradt):
ThisQuench.write('%16i%30.15E\n'%(1, extradt))
else:
ThisQuench.write('%16i%30.15E\n'%(0, extradt))
ThisQuench.write(scalestr + '\n')
ThisQuench.write('%16i\n'%(nchanged))
# Indices of parameters being changed
chngstr = ''
for hparam in myHparams:
chngstr += '%16i'%(H.IntHparam[hparam])
ThisQuench.write(chngstr + '\n')
# Write regular time steps
# ------------------------
# Add first time step already for Trotter decomposition
for j in range(1, nsteps + 1):
# Write parameters for the time evolution
# .......................................
hstr = ''
for k in range(sa):
cntr = 0
for hparams in myHparams:
h = 0.0 * self.data[ii]['funcs'][cntr](0.0)
if isinstance(h, (list,np.ndarray)):
for m in range(sb):
h += g[k,m] * self.data[ii]['funcs'][cntr](time +
x[m] * self.data[ii]['deltat']) / scalefac[k]
for l in range(len(h)):
hstr += '%30.15E'%(h[l])
hstr += '\n'
else:
for m in range(sb):
h += g[k,m] * self.data[ii]['funcs'][cntr](time +
x[m] * self.data[ii]['deltat']) / scalefac[k]
hstr += '%30.15E\n'%(h)
cntr += 1
if not PostProcess:
ThisQuench.write(hstr)
time += self.data[ii]['deltat']
# Write parameters for energy measurement
# .......................................
if(j%self.data[ii]['stepsforoutput'] == 0):
hstr = ''
# write actual h value at the given time for energy evaluation
cntr = 0
for hparams in myHparams:
h = self.data[ii]['funcs'][cntr](time)
if isinstance(h, (list,np.ndarray)):
# Space dependent case
for l in range(len(h)):
hstr += '%30.15E'%(h[l])
hstr += '\n'
else:
# Space independent case
hstr += '%30.15E\n'%(h)
cntr += 1
if(not PostProcess):
ThisQuench.write(hstr)
nmeasurements += 1
ThisQuench.flush()
# Write extra time step
# ---------------------
if hasextradt:
# Write parameters for the time evolution
# .......................................
hstr = ''
for k in range(sa):
cntr = 0
for hparams in myHparams:
h = 0.0 * self.data[ii]['funcs'][cntr](0.0)
if isinstance(h, (list, np.ndarray)):
for m in range(sb):
h += g[k,m] * self.data[ii]['funcs'][cntr](time +
x[m] * self.data[ii]['deltat']) / scalefac[k]
for l in range(len(h)):
hstr += '%30.15E'%(h[l])
hstr += '\n'
else:
for m in range(sb):
h += g[k,m] * self.data[ii]['funcs'][cntr](time +
x[m] * self.data[ii]['deltat']) / scalefac[k]
hstr += '%30.15E\n'%(h)
time += extradt
if not PostProcess:
ThisQuench.write(hstr)
# Write parameters for the energy measurement
# ...........................................
hstr = ''
cntr = 0
for hparams in myHparams:
h = self.data[ii]['funcs'][cntr](time)
if isinstance(h, (list, np.ndarray)):
for l in range(len(h)):
hstr += '%30.15E'%(h[l])
hstr += '\n'
else:
hstr += '%30.15E\n'%(h)
if(not PostProcess):
ThisQuench.write(hstr)
cntr += 1
nmeasurements += 1
# p['nmeasurements'].append(nmeasurements)
if not PostProcess:
ThisQuench.close()
# if('QT' in p): qtcopy(ThisQuenchname, p['QT'])
return
[docs]class QuantumCircuits():
"""
Class defining the time evolution with quantum circuits or gates.
**Variables**
subcircuit : dict
A dictionary storing subcircuits identified via keys.
evolution : list
Evolution contains a list of keys, where each key must refer
to one subcircuit to be applied.
"""
def __init__(self):
self.subcircuit = {}
self.evolution = []
[docs] def classname(self):
"""
Return the string ``QuantumCircuits``, which is the name of the
class.
"""
return 'QuantumCircuits'
[docs] def __len__(self):
"""
Return the length of a quantum circuit identified by the number
of subcircuits applied.
"""
return len(self.evolution)
[docs] def __setitem__(self, key, value):
"""
Add or renew a subcircuit with dictionary syntax.
**Arguments**
key :
Identifier for the subcircuit being added or renewed.
value : :py:class:`Dynamics.Subcircuit`
The subcircuit to be added. If not of this class, it must
provide similar class methods.
"""
self.subcircuit[key] = value
[docs] def __getitem__(self, key):
"""
Get a subcircuit.
**Arguments**
key :
Identifier for the subcircuit.
"""
return self.subcircuit[key]
[docs] def add(self, key):
"""
Add a subcircuit to the evolution.
**Arguments**
key :
Identifier for the subcircuit to be applied.
"""
self.evolution.append(key)
[docs] def __iadd__(self, key):
"""
Add a subcircuit to the evolution.
**Arguments**
key :
Identifier for the subcircuit to be applied.
"""
self.add(key)
return self
[docs] def __iter__(self):
"""
Iterate over the subcircuits for the time evolution return
the subcircuit.
"""
for ii in range(len(self)):
yield self.subcircuit[self.evolution[ii]]
[docs]class Subcircuit():
"""
Class defining one application pattern for a circuit to a
set of sites.
**Arguments**
Ops : instance of :py:class:`Ops.OperatorList`
Operators containing the quantum gates specified via their
string identifier.
**Variables**
gates : list of string
The string identifier for the gates in the order in which
they should be applied. First string is applied first.
sites : listed list of integers
Each element in the list is a list of integers. The integers
specify to which sites the gate on the corresponding position'
in the gates list is applied to.
time : list of floats
The time steps associate with the corresponding item in the
gates list.
"""
def __init__(self, Ops):
self.Ops = Ops
self.gates = []
self.sites = []
self.time = []
[docs] def add(self, gate, sites, time=1.0):
"""
Add another gate to the sub-circuit.
**Arguments**
gate : str
String identifier for the gate.
sites : list of integers
Specifies the sites where the gate should be applied
time : float, optional
Time associated with the gate.
Default to 1.0
"""
self.gates.append(gate)
self.sites.append(sites)
self.time.append(time)
return
[docs] def __iadd__(self, args):
"""
Short-cut to add a gate.
**Arguments**
args : tuple or list
If two entries, then the string identifier for the gate and
the list of integers where the gate is applied. If three
entries, the last entry is the time.
"""
if(len(args) == 2):
self.add(args[0], args[1])
elif(len(args) == 3):
self.add(args[0], args[1], time=args[2])
else:
raise Exception("Unknown length of arguments")
return self
[docs] def __iter__(self):
"""
Iterator returns the operator as 2d array, the list of sites,
and the time.
"""
for ii in range(len(self.gates)):
yield self.Ops[self.gates[ii]], self.sites[ii], self.time[ii]
[docs]def swp_subcircuit(Opl, Opb, Opr, ll, Operators):
"""
Build a subcircuit for a three qubit gate in the SWP configuration.
**Arguments**
Opl : string
String identifier for the two qubit gate for the left boundary.
Opb : string
String identifier for the three qubit gate in the bulk of the system.
Opr : string
String identifier for the two qubit gate for the right boundary.
ll : int
Number of sites in the system.
Operators : instance of :py:class:`Ops.OperatorList`
Operators containing the quantum gates specified via their
string identifier.
"""
QCsub = Subcircuit(Operators)
QCsub += (Opl, [0, 1], 0.0)
for ii in range(ll - 2):
QCsub += (Opb, [ii, ii + 1, ii + 2], 0.0)
QCsub += (Opr, [ll - 2, ll - 1], 1.0)
return QCsub
[docs]def alt_subcircuit(Opl, Opb, Opr, ll, Operators):
"""
Build a subcircuit for a three qubit gate in the ALT configuration.
**Arguments**
Opl : string
String identifier for the two qubit gate for the left boundary.
Opb : string
String identifier for the three qubit gate in the bulk of the system.
Opr : string
String identifier for the two qubit gate for the right boundary.
ll : int
Number of sites in the system.
Operators : instance of :py:class:`Ops.OperatorList`
Operators containing the quantum gates specified via their
string identifier.
"""
QCsub = Subcircuit(Operators)
# First layer
# -----------
QCsub += (Opl, [0, 1], 1.0)
ii = 3
while(ii < ll):
QCsub += (Opb, [ii - 2, ii - 1, ii], 0.0)
ii += 2
if(ll % 2 == 1):
QCsub += (Opr, [ll - 2, ll - 1], 0.0)
# Second layer
# ------------
ii = 2
while(ii < ll):
QCsub += (Opb, [ii - 2, ii - 1, ii], 0.0)
ii += 2
if(ll % 2 == 0):
QCsub += (Opr, [ll - 2, ll - 1], 0.0)
return QCsub
[docs]def blk_subcircuit(Opl, Opb, Opr, ll, Operators):
"""
Build a subcircuit for a three qubit gate in the BLK configuration.
**Arguments**
Opl : string
String identifier for the two qubit gate for the left boundary.
Opb : string
String identifier for the three qubit gate in the bulk of the system.
Opr : string
String identifier for the two qubit gate for the right boundary.
ll : int
Number of sites in the system.
Operators : instance of :py:class:`Ops.OperatorList`
Operators containing the quantum gates specified via their
string identifier.
"""
QCsub = Subcircuit(Operators)
# First layer
# -----------
QCsub += (Opl, [0, 1], 1.0)
ii = 4
while(ii < ll):
QCsub += (Opb, [ii - 2, ii - 1, ii], 0.0)
ii += 3
if(ll % 3 == 1):
QCsub += (Opr, [ll - 2, ll - 1], 0.0)
# Second layer
# ------------
ii = 2
while(ii < ll):
QCsub += (Opb, [ii - 2, ii - 1, ii], 0.0)
ii += 3
if(ll % 3 == 2):
QCsub += (Opr, [ll - 2, ll - 1], 0.0)
# Third layer
# -----------
ii = 3
while(ii < ll):
QCsub += (Opb, [ii - 2, ii - 1, ii], 0.0)
ii += 3
if(ll % 3 == 0):
QCsub += (Opr, [ll - 2, ll - 1], 0.0)
return QCsub
[docs]def WriteDynamics(H, p, qtcopy, PostProcess=False):
"""Write out the dynamics to a Fortran-readable file.
**Arguments**
H : instance of `mps.MPO`
contains the information about the system Hamiltonian
p : dictionary
containing the static parameters of a simulation as `job_ID`,
`unique_ID`, etc.
qtcopy : callable
copies files for quantum trajectories. Takes two arguments: filename
and number of trajectories (eqauls to the number of copies needed).
PostProcess : boolean, optional
`False` : write output file used in fortran simulation
`True` : used within python postprocessing script
default to `False`
**Details**
The first file (meta file) should be structured in the following way:
1) Number of quenches `n`
2) `n` lines with the corresponding quench data method e.g. `krylov`
The file for each quench contains the following information
1) number of time steps, step between measurements, time step dt
2) number of exponentials (for time-ordering)
3) 0 or 1, extradt (float, time steps size, only Krylov)
4) scale factors (time-ordering)
5) number of changed Hamiltonian parameters (integer)
6) indices of the changed Hamiltonian parameter
Looping over the times steps:
7a) Hamiltonian : looping over the decomposition, then one line
for each changing Hamiltonian parameter with one float (space
independent) or L floats (space dependent)
8) Measurements : for every changed parameter one
line either containing one float (space independent) or couplings for
all sites.
For the extradt (one time step):
9a) Hamiltonian for extradt: same string as for Hamiltonian
10) Measurements for extra dt: one float or line
of floats.
"""
# Filenames, create subdirectory
writedir = p['Write_Directory']
totalID = p['job_ID'] + p['unique_ID']
subdir = p['Output_Directory'] + totalID + 'ObsOutDynamics'
if(not os.path.exists(subdir)):
makedirs(subdir)
if(isinstance(p['Quenches'], QuantumCircuits)): return
# Get fundamental data
nquenches = len(p['Quenches'].data)
p['nquenches'] = nquenches
Quenches = p['Quenches'].data
# Check for errors
# ----------------
if not PostProcess:
files=[]
MetaFileName=writedir+totalID+'Dynamics.dat'
p['int_files'].append(MetaFileName)
files.append(MetaFileName)
MetaFile=open(MetaFileName,'w')
#write out quenches
MetaFile.write('%16i'%(nquenches)+'\n')
#Write out what method is being used
for i in range(nquenches):
mkey = Quenches[i]['ConvergenceParameters'].get_method()
MetaFile.write('%16i'%(mkey) + '\n')
MetaFile.close()
if('QT' in p): qtcopy(MetaFileName, p['QT'])
H.IndexHParams()
thehparams=H.getHparams()
p['nmeasurements']=[]
time = 0.0
for i in range(nquenches):
# Order is contained in convergence parameters
myTO = Quenches[i]['ConvergenceParameters'].get_TO()
splitTO = myTO.split('_')
# Commutator-free magnus expansion setup
# --------------------------
order = int(splitTO[0])
if(len(splitTO)>1):
nexp = int(splitTO[1])
x, g = SetupGLQ(order, nexp)
else:
x, g = SetupGLQ(order)
if order==2:
nexp=1
elif order==6:
nexp=5
# Set up GL parameters
sa = g.shape[0]
sb = g.shape[1]
scalestr = ''
# Get the "scale factors" s_i=\sum_{m}g_{i,m} to scale the time
# step by. This ensures the proper behavior of all terms which
# remain constant. For those which are not constant, divide
# through by the scale factor
scalefac = np.zeros(sa)
for k in range(sa):
for m in range(sb):
scalefac[k]+=g[k,m]
scalestr+='%30.15E'%(scalefac[k])
myHparams=Quenches[i]['Hparams']
ftmp = writedir + totalID + 'tMPS_' + str(i + 1) + '.dat'
Quenches[i]['ConvergenceParameters'].write(ftmp)
if(not PostProcess and ('QT' in p)): qtcopy(ftmp, p['QT'])
nmeasurements=0
nsteps=int(floor(Quenches[i]['time']/Quenches[i]['deltat']))
print('nsteps', nsteps)
extradt=Quenches[i]['time']-nsteps*Quenches[i]['deltat']
if extradt>0.000001:
hasextradt=True
else:
hasextradt=False
#Number of parameters being changed
nchanged=len(myHparams)
if not PostProcess:
ThisQuenchname=writedir+totalID+'Quench_'+str(i+1)+'.dat'
p['int_files'].append(ThisQuenchname)
files.append(ThisQuenchname)
ThisQuench=open(ThisQuenchname,'w')
ThisQuench.write('%16i'%(nsteps)+'%16i'%(Quenches[i]['stepsforoutput'])+'%30.15E'%(Quenches[i]['deltat'])+'\n')
ThisQuench.write('%16i'%(nexp)+'\n')
if(hasextradt):
ThisQuench.write('%16i'%(1)+'%30.15E'%(extradt)+'\n')
else:
ThisQuench.write('%16i'%(0)+'%30.15E'%(extradt)+'\n')
ThisQuench.write(scalestr+'\n')
ThisQuench.write('%16i'%(nchanged)+'\n')
#Indices of parameters being changed
chngstr=''
for hparam in myHparams:
chngstr+='%16i'%(H.IntHparam[hparam])
ThisQuench.write(chngstr+'\n')
# Write regular time steps
# ------------------------
# Add first time step already for Trotter decomposition
for j in range(1,nsteps+1):
# Write parameters for the time evolution
# .......................................
hstr=''
for k in range(sa):
cntr = 0
for hparams in myHparams:
h = 0.0*Quenches[i]['funcs'][cntr](0.0)
if isinstance(h, (list,np.ndarray)):
for m in range(sb):
h += g[k,m] * Quenches[i]['funcs'][cntr](time +
x[m] * Quenches[i]['deltat']) / scalefac[k]
for l in range(len(h)):
hstr += '%30.15E'%(h[l])
hstr += '\n'
else:
for m in range(sb):
h += g[k,m] * Quenches[i]['funcs'][cntr](time +
x[m] * Quenches[i]['deltat']) / scalefac[k]
hstr += '%30.15E'%(h) + '\n'
cntr += 1
if not PostProcess:
ThisQuench.write(hstr)
time+=Quenches[i]['deltat']
# Write parameters for energy measurement
# .......................................
if j%Quenches[i]['stepsforoutput']==0:
hstr=''
#write actual h value at the given time for energy evaluation
cntr=0
for hparams in myHparams:
h=Quenches[i]['funcs'][cntr](time)
if isinstance(h,(list,np.ndarray)):
# Space dependent case
for l in range(len(h)):
hstr+='%30.15E'%(h[l])
hstr+='\n'
else:
# Space independent case
hstr += '%30.15E'%(h) + '\n'
cntr+=1
if(not PostProcess):
ThisQuench.write(hstr)
nmeasurements+=1
# Write extra time step
# ---------------------
if hasextradt:
# Write parameters for the time evolution
# .......................................
hstr=''
for k in range(sa):
cntr = 0
for hparams in myHparams:
h = 0.0 * Quenches[i]['funcs'][cntr](0.0)
if isinstance(h, (list,np.ndarray)):
for m in range(sb):
h += g[k,m] * Quenches[i]['funcs'][cntr](time +
x[m] * Quenches[i]['deltat']) / scalefac[k]
for l in range(len(h)):
hstr += '%30.15E'%(h[l])
hstr += '\n'
else:
for m in range(sb):
h += g[k,m] * Quenches[i]['funcs'][cntr](time +
x[m] * Quenches[i]['deltat']) / scalefac[k]
hstr += '%30.15E'%(h) + '\n'
time += extradt
if not PostProcess:
ThisQuench.write(hstr)
# Write parameters for the energy measurement
# ...........................................
hstr=''
cntr=0
for hparams in myHparams:
h=Quenches[i]['funcs'][cntr](time)
if isinstance(h,(list,np.ndarray)):
for l in range(len(h)):
hstr+='%30.15E'%(h[l])
hstr+='\n'
else:
hstr+='%30.15E'%(h)+'\n'
if(not PostProcess):
ThisQuench.write(hstr)
cntr+=1
nmeasurements+=1
p['nmeasurements'].append(nmeasurements)
if not PostProcess:
ThisQuench.close()
if('QT' in p): qtcopy(ThisQuenchname, p['QT'])
return