import copy
import numpy as np
from ..tt.TTs import TTs
from .opett import (zGetSegLeftSt, zGetSegLeft, zGetSegRightEn, zGetSegRight, zGetKSt,
zGetKEn, zGetK, zGetS, zGetSK, zGetKS, zQRLeft, zQRRight
)
[docs]
class timeEvolution():
"""Time evolution of the tensor-train density operator via TDVP.
Attributes
----------
dt : float
Step width.
numCore : int
Number of cores.
segArray : list
Segment tensors used for updating each core.
ZNRK : int
Number of stages in the Runge-Kutta scheme.
zA : numpy.ndarray
Runge-Kutta coefficient array A.
zB : numpy.ndarray
Runge-Kutta coefficient array B.
zC : numpy.ndarray
Runge-Kutta coefficient array C.
"""
def __init__(self, TTsIni: TTs, dt: float, isRK13: bool):
"""Initialize the time evolution object.
Parameters
----------
TTsIni : TTs.TTs
Initialized MPS and MPO.
dt : float
Step width for forward/backward time integration.
isRK13 : bool
Runge-Kutta method selection.
``True``: 13-stage 5th-order Runge-Kutta.
``False``: 5-stage 4th-order Runge-Kutta.
"""
self.dt = dt
self.numH = TTsIni.numH
self.getPrefactors = TTsIni.getPrefactors
self.segArray = self.zInitSegment(TTsIni)
if isRK13: # RK13-5
self.ZNRK = 13
self.zA = np.array([
0,
-0.33672143119427413,
-1.2018205782908164,
-2.6261919625495068,
-1.5418507843260567,
-0.2845614242371758,
-0.1700096844304301,
-1.0839412680446804,
-11.61787957751822,
-4.5205208057464192,
-35.86177355832474,
-0.000021340899996007288,
-0.066311516687861348
])
self.zB = np.array([
0.069632640247059393,
0.088918462778092020,
1.0461490123426779,
0.42761794305080487,
0.20975844551667144,
-0.11457151862012136,
-0.01392019988507068,
4.0330655626956709,
0.35106846752457162,
-0.16066651367556576,
-0.0058633163225038929,
0.077296133865151863,
0.054301254676908338
])
self.zC = np.array([
0,
0.069632640247059393,
0.12861035097891748,
0.34083022189561149,
0.54063706308495402,
0.59927749518613931,
0.49382042519248519,
0.48207852767699775,
0.82762865209834452,
0.82923953914857933,
0.67190565554748019,
0.87194975193167848,
0.94930216564503562
])
else: # RK5-4
self.ZNRK = 5
self.zA = np.array([
0,
-567301805773 / 1357537059087,
-2404267990393 / 2016746695238,
-3550918686646 / 2091501179385,
-1275806237668 / 842570457699
])
self.zB = np.array([
1432997174477 / 9575080441755,
5161836677717 / 13612068292357,
1720146321549 / 2090206949498,
3134564353537 / 4481467310338,
2277821191437 / 14882151754819
])
self.zC = np.array([
0,
1432997174477 / 9575080441755,
2526269341429 / 6820363962896,
2006345519317 / 3224310063776,
2802321613138 / 2924317926251
])
[docs]
def zInitSegment(self, TTsIni: TTs):
"""Initialize the segment tensors from the rightmost core.
Parameters
----------
TTsIni : TTs.TTs
Initialized MPS and MPO.
Returns
-------
segArray : list
List of segment tensors used for updating each core.
"""
segArray = [[None for _ in range(TTsIni.numCore)]
for _ in range(self.numH)]
i = TTsIni.numCore-1
for j in range(self.numH):
segArray[j][i] = zGetSegRightEn(TTsIni.rho[i],
TTsIni.H[j, i])
for i in range(TTsIni.numCore-2, 0, -1):
for j in range(self.numH):
segArray[j][i] = zGetSegRight(TTsIni.rho[i],
TTsIni.H[j, i],
segArray[j][i+1])
i = 0
for j in range(self.numH):
segArray[j][i] = np.zeros(TTsIni.rho[i].bondDimL**2
* TTsIni.H[j, i].bondDimL,
dtype = np.complex128)
return segArray
[docs]
def zTTTimeEvo(self, rho, H, time, stepNum):
"""Perform one forward-backward TDVP sweep to advance ``rho`` by ``dt``.
Parameters
----------
rho : numpy.ndarray
1-D array of :class:`~ttheom.tt.tt.zTT` (MPS); overwritten in place.
H : numpy.ndarray
2-D array of :class:`~ttheom.tt.tt.zTT` (MPO) representing the Hamiltonian.
time : float
Current time.
stepNum : int
Current step number.
"""
numCore = len(rho)
intBonds = [0] * (self.numH)
self.stepNum = stepNum
# Forward sweep
i = 0
self.zKRK4St(rho[i], H[:, i], i, time)
rhoTmp, S = zQRLeft(rho[i])
rho[i] = copy.deepcopy(rhoTmp)
for j in range(self.numH):
self.segArray[j][i] = zGetSegLeftSt(rho[i], H[j, i])
intBonds[j] = H[j, i].bondDimR
self.zSRK4(S, i, rho[i].bondDimR, intBonds, time)
for i in range(1, numCore - 1):
zGetSK(S, rho[i])
self.zKRK4(rho[i], H[:, i], i, time)
rhoTmp, S = zQRLeft(rho[i])
rho[i] = copy.deepcopy(rhoTmp)
for j in range(self.numH):
self.segArray[j][i] = zGetSegLeft(
rho[i], H[j, i], self.segArray[j][i - 1])
intBonds[j] = H[j, i].bondDimR
self.zSRK4(S, i, rho[i].bondDimR, intBonds, time)
i = numCore - 1
zGetSK(S, rho[i])
self.zKRK4En(rho[i], H[:, i], i, time)
# Backward sweep
self.zKRK4En(rho[i], H[:, i], i, time + self.dt)
rhoTmp, S = zQRRight(rho[i])
rho[i] = copy.deepcopy(rhoTmp)
for j in range(self.numH):
self.segArray[j][i] = zGetSegRightEn(rho[i], H[j, i])
intBonds[j] = H[j, i].bondDimL
self.zSRK4(S, i - 1, rho[i].bondDimL, intBonds, time + self.dt)
for i in range(numCore - 2, 0, -1):
zGetKS(rho[i], S)
self.zKRK4(rho[i], H[:, i], i, time + self.dt)
rhoTmp, S = zQRRight(rho[i])
rho[i] = copy.deepcopy(rhoTmp)
for j in range(self.numH):
self.segArray[j][i] = zGetSegRight(
rho[i], H[j, i], self.segArray[j][i + 1])
intBonds[j] = H[j, i].bondDimL
self.zSRK4(S, i - 1, rho[i].bondDimL,
intBonds, time + self.dt)
i = 0
zGetKS(rho[i], S)
self.zKRK4St(rho[i], H[:, i], i, time + self.dt)
[docs]
def zKRK4St(self, rho, H, coreIdx, time):
"""Apply the Runge-Kutta update to the first (starting) core.
Parameters
----------
rho : tt.zTT
MPS core; overwritten in place.
H : numpy.ndarray
Array of MPO cores.
coreIdx : int
Index of the current core.
time : float
Current time.
"""
dumCore1 = np.zeros(rho.core.shape, dtype=np.complex128)
for i in range(self.ZNRK):
timeTmp = time + self.dt * self.zC[i]
preFact = self.getPrefactors(self.dt, timeTmp, self.stepNum)
dumCore1 *= self.zA[i]
for j in range(self.numH):
dumCore2 = zGetKSt(rho, H[j],
self.segArray[j][coreIdx+1])
dumCore1 += preFact[j] * dumCore2
rho.core += self.zB[i] * dumCore1
[docs]
def zKRK4En(self, rho, H, coreIdx, time):
"""Apply the Runge-Kutta update to the last (ending) core.
Parameters
----------
rho : tt.zTT
MPS core; overwritten in place.
H : numpy.ndarray
Array of MPO cores.
coreIdx : int
Index of the current core.
time : float
Current time.
"""
dumCore1 = np.zeros(rho.core.shape, dtype=np.complex128)
for i in range(self.ZNRK):
timeTmp = time + self.dt * self.zC[i]
preFact = self.getPrefactors(self.dt, timeTmp, self.stepNum)
dumCore1 *= self.zA[i]
for j in range(self.numH):
dumCore2 = zGetKEn(rho, H[j],
self.segArray[j][coreIdx-1])
dumCore1 += preFact[j] * dumCore2
rho.core += self.zB[i] * dumCore1
[docs]
def zKRK4(self, rho, H, coreIdx, time):
"""Apply the Runge-Kutta update to an intermediate core.
Parameters
----------
rho : tt.zTT
MPS core; overwritten in place.
H : numpy.ndarray
Array of MPO cores.
coreIdx : int
Index of the current core.
time : float
Current time.
"""
dumCore1 = np.zeros(rho.core.shape, dtype=np.complex128)
for i in range(self.ZNRK):
timeTmp = time + self.dt * self.zC[i]
preFact = self.getPrefactors(self.dt, timeTmp, self.stepNum)
dumCore1 *= self.zA[i]
for j in range(self.numH):
dumCore2 = zGetK(rho, H[j],
self.segArray[j][coreIdx-1],
self.segArray[j][coreIdx+1])
dumCore1 += preFact[j] * dumCore2
rho.core += self.zB[i] * dumCore1
[docs]
def zSRK4(self, S, coreIdx, rhoR, HRs, time):
"""Apply the Runge-Kutta update to the bond matrix S.
Parameters
----------
S : numpy.ndarray
Bond matrix; overwritten in place.
coreIdx : int
Core index to the left of the bond.
rhoR : int
Right bond dimension of ``rho``.
HRs : list of int
Right bond dimensions of each Hamiltonian term.
time : float
Current time.
"""
dumS1 = np.zeros(S.shape, dtype=np.complex128)
for i in range(self.ZNRK):
timeTmp = time + self.dt * self.zC[i]
preFact = -self.getPrefactors(self.dt, timeTmp, self.stepNum)
dumS1 *= self.zA[i]
for j in range(self.numH):
dumS2 = zGetS(rhoR, HRs[j], S,
self.segArray[j][coreIdx],
self.segArray[j][coreIdx+1])
dumS1 += preFact[j] * dumS2
S += self.zB[i] * dumS1
[docs]
def zRightOrth(rho):
"""Right-orthogonalize all MPS cores (initialization sweep).
Parameters
----------
rho : numpy.ndarray
1-D array of :class:`~ttheom.tt.tt.zTT` MPS cores; modified in place.
"""
numCore = len(rho)
# Start from the last core
i = numCore - 1
rhoTmp, S = zQRRight(rho[i])
rho[i] = copy.deepcopy(rhoTmp)
# Sweep from right to left
for i in range(numCore - 2, 0, -1):
zGetKS(rho[i], S)
rhoTmp, S = zQRRight(rho[i])
rho[i] = copy.deepcopy(rhoTmp)
# Final propagation of S to the first core
i = 0
zGetKS(rho[i], S)