commit a425029cc2f8da8b54478f2b3539d8b4926114a3
parent 43e837ff846bd38b7766877c64515bc8fa8efd9e
Author: Miroslav Urbanek <mu@miroslavurbanek.com>
Date: Wed, 28 Jun 2017 19:43:34 +0200
Add support for two-dimensional models
Two-dimensional states are represented as tree tensor network states
(TTNS). The time-evolution algorithm rearranges the network to apply a
two-site operator to each pair of adjacent sites.
Diffstat:
5 files changed, 1014 insertions(+), 11 deletions(-)
diff --git a/tebdol/mps.lisp b/tebdol/mps.lisp
@@ -47,12 +47,6 @@
1
0)))))
-(defun val (tensor)
- (let ((s (tensor-sectors tensor)))
- (if s
- (aref (sector-array (first s)))
- +blas-0+)))
-
(defun right-singvals-tensor (index)
(functional-tensor
(list
@@ -112,7 +106,7 @@
#'(lambda (s) (declare (ignore s)) 1)))
(setf l (tensor-contract l '(0 1) r '(0 1)))))
(if (= rank (1- size))
- (val l)
+ (tensor-scalar l)
(mpi-send-object l (1+ rank)))))
(defun mps-expectation-value (size rank length mpo submps)
@@ -145,7 +139,7 @@
#'(lambda (s) (declare (ignore s)) 1)))
(setf l (tensor-contract l '(0 1 2) r '(0 1 2)))))
(if (= rank (1- size))
- (val l)
+ (tensor-scalar l)
(mpi-send-object l (1+ rank)))))
(defun save-metadata (dimension length particles pathname)
diff --git a/tebdol/tebd.lisp b/tebdol/tebd.lisp
@@ -1,6 +1,8 @@
(defpackage :tebd
(:use :common-lisp :blas :tensor :mpi :mps)
- (:export :make-mph-propagator
+ (:export :make-single-site-propagator
+ :make-double-site-propagator
+ :make-mph-propagator
:make-uniform-propagator
:tebd-evolve))
diff --git a/tebdol/tebdol.asd b/tebdol/tebdol.asd
@@ -4,7 +4,7 @@
:description "TEBDOL is an implementation of the time-evolving block
decimation algorithm for simulating evolution of ultracold atoms in
optical lattices."
- :version "0.2.0"
+ :version "0.3.0"
:author "Miroslav Urbanek <miroslav.urbanek@mff.cuni.cz>"
:licence "GPL-3+"
:serial t
@@ -18,4 +18,5 @@
(:file "mps")
(:file "part")
(:file "tebd")
- (:file "bhm")))
+ (:file "bhm")
+ (:file "ttns")))
diff --git a/tebdol/tensor.lisp b/tebdol/tensor.lisp
@@ -24,6 +24,8 @@
:tensor-rank
:tensor-dimensions
:tensor-total-size
+ :tensor-scalar
+ :tensor-normalize
:tensor-conjugate
:tensor-contract
:tensor-permute
@@ -142,6 +144,35 @@
(defun tensor-total-size (tensor)
(reduce #'+ (tensor-sectors tensor) :key #'(lambda (s) (array-total-size (sector-array s)))))
+(defun tensor-scalar (tensor)
+ (let ((s (tensor-sectors tensor)))
+ (if s
+ (aref (sector-array (first s)))
+ +blas-0+)))
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+;; tensor-normalize
+
+(defun tensor-norm (tensor)
+ (sqrt
+ (loop
+ for s in (tensor-sectors tensor)
+ sum (expt (array-norm (sector-array s)) 2))))
+
+(defun tensor-normalize (tensor)
+ (let ((n (tensor-norm tensor)))
+ (make-tensor
+ :indices
+ (tensor-indices tensor)
+ :sectors
+ (mapcar
+ #'(lambda (s)
+ (make-sector
+ :numbers (sector-numbers s)
+ :array (array-scalar-/ (sector-array s) n)))
+ (tensor-sectors tensor)))))
+
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; tensor-conjugate
diff --git a/tebdol/ttns.lisp b/tebdol/ttns.lisp
@@ -0,0 +1,975 @@
+(defpackage :ttns
+ (:use :common-lisp :tensor :bhm :tebd)
+ (:export :make-vacuum-ttns
+ :ttns-create-particle
+ :ttns-overlap
+ :compress-ttns
+ :make-ssps
+ :make-dsps
+ :ttns-tebd-evolve
+ :all-pn-operators
+ :all-pn
+ :print-pn))
+
+(in-package :ttns)
+
+(defstruct site
+ id
+ subscript
+ bonds)
+
+(defstruct ttns
+ width
+ height
+ lattice
+ tensors)
+
+(defun site-by-id (ttns id)
+ (find id (ttns-lattice ttns) :key #'site-id :test #'equal))
+
+(defun site-by-bond (ttns site bond)
+ (let ((i (cdr (assoc bond (site-bonds site)))))
+ (if i (svref (ttns-lattice ttns) i))))
+
+(defun tensor-by-site (ttns site)
+ (svref (ttns-tensors ttns) (site-subscript site)))
+
+(defun (setf tensor-by-site) (tensor ttns site)
+ (setf
+ (svref (ttns-tensors ttns) (site-subscript site))
+ tensor))
+
+(defun tensor-at-site (ttns id)
+ (svref (ttns-tensors ttns) (site-subscript (site-by-id ttns id))))
+
+(defun (setf tensor-at-site) (tensor ttns id)
+ (setf
+ (svref (ttns-tensors ttns) (site-subscript (site-by-id ttns id)))
+ tensor))
+
+(defun make-rectangular-lattice (width height)
+ (flet ((subscript (i j)
+ (if (and (>= i 0) (< i width) (>= j 0) (< j height))
+ (+ i (* j width)))))
+ (let ((a (make-array (* width height))))
+ (dotimes (i width a)
+ (dotimes (j height)
+ (let ((s (subscript i j)))
+ (setf (svref a s)
+ (make-site
+ :id (list i j)
+ :subscript s
+ :bonds (list
+ (cons :x- (subscript (1- i) j))
+ (cons :x+ (subscript (1+ i) j))
+ (cons :y- (subscript i (1- j)))
+ (cons :y+ (subscript i (1+ j))))))))))))
+
+(defun make-vacuum-ttns (width height dimension)
+ (let* ((l (make-rectangular-lattice width height))
+ (a (make-array (length l))))
+ (dotimes (i (length l)
+ (make-ttns :width width :height height :lattice l :tensors a))
+ (setf (svref a i)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1)))
+ #'(lambda (s) (declare (ignore s)) 1))))))
+
+(defun identity-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s)
+ (declare (ignore s)) 1)))
+
+(defun xx-connecting-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '(-1) :dimension 1))
+ (list (make-segment :numbers '( 1) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s) (declare (ignore s)) 1)))
+
+(defun yy-connecting-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '(-1) :dimension 1))
+ (list (make-segment :numbers '( 1) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s) (declare (ignore s)) 1)))
+
+(defun xy-connecting-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '(-1) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 1) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s) (declare (ignore s)) 1)))
+
+(defun x-creation-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '(-1) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s)
+ (let ((n (- (first (subscript-numbers (sixth s))))))
+ (sqrt (1+ n))))))
+
+(defun y-creation-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (list (make-segment :numbers '(-1) :dimension 1))
+ (list (make-segment :numbers '( 0) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s)
+ (let ((n (- (first (subscript-numbers (sixth s))))))
+ (sqrt (1+ n))))))
+
+(defun pn-site-operator (dimension)
+ (functional-tensor
+ (list
+ (ket-physical-index dimension)
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (bra-physical-index dimension))
+ #'(lambda (s)
+ (- (first (subscript-numbers (sixth s)))))))
+
+(defun identity-operator (width height dimension)
+ (let* ((l (make-rectangular-lattice width height))
+ (a (make-array (length l)))
+ (o (identity-site-operator dimension)))
+ (dotimes (i (length l)
+ (make-ttns
+ :width width
+ :height height
+ :lattice l
+ :tensors a))
+ (setf (svref a i) o))))
+
+(defun apply-and-fuse (operator ttns)
+ (let* ((l (length (ttns-lattice ttns)))
+ (o (ttns-tensors operator))
+ (p (ttns-tensors ttns))
+ (a (make-array l)))
+ (dotimes (i l
+ (make-ttns
+ :width (ttns-width ttns)
+ :height (ttns-height ttns)
+ :lattice (ttns-lattice ttns)
+ :tensors a))
+ (setf (svref a i)
+ (tensor-fuse
+ (tensor-permute
+ (tensor-contract (svref o i) 5 (svref p i) 0)
+ '(0 1 3 5 7 2 4 6 8))
+ '(1 3 5 7)
+ '(:normal :reverse :normal :reverse :normal))))))
+
+(defun ttns-create-particle (ttns x y dimension)
+ (flet ((f (i j o)
+ (let ((s (site-subscript (site-by-id ttns (list i j))))
+ (p (ttns-tensors ttns)))
+ (setf (svref p s)
+ (tensor-fuse
+ (tensor-permute
+ (tensor-contract o 5 (svref p s) 0)
+ '(0 1 3 5 7 2 4 6 8))
+ '(1 3 5 7)
+ '(:normal :reverse :normal :reverse :normal))))))
+ (loop
+ with o = (xx-connecting-site-operator dimension)
+ for i below x
+ do (f i 0 o))
+ (loop
+ with o = (yy-connecting-site-operator dimension)
+ for i from 1 below y
+ do (f x i o))
+ (if (zerop y)
+ (f x y (x-creation-site-operator dimension))
+ (progn
+ (f x 0 (xy-connecting-site-operator dimension))
+ (f x y (y-creation-site-operator dimension))))))
+
+(defun pn-operator (width height x y dimension)
+ (let ((o (identity-operator width height dimension)))
+ (setf (tensor-at-site o (list x y))
+ (pn-site-operator dimension))
+ o))
+
+(defun mps-boundary-tensor ()
+ (let ((m (list (make-segment :numbers '(0) :dimension 1))))
+ (functional-tensor
+ (list m m m)
+ #'(lambda (s) (declare (ignore s)) 1))))
+
+(defun mps-overlap-boundary-tensor ()
+ (let ((m (list (make-segment :numbers '(0) :dimension 1))))
+ (functional-tensor
+ (list m m)
+ #'(lambda (s) (declare (ignore s)) 1))))
+
+(defun terminate-dummy-index (tensor index)
+ (tensor-contract
+ tensor
+ index
+ (functional-tensor
+ (list (list (make-segment :numbers '(0) :dimension 1)))
+ #'(lambda (s) (declare (ignore s)) 1))
+ 0))
+
+(defun terminate-dummy-indices (tensor indices)
+ (reduce #'terminate-dummy-index indices :initial-value tensor))
+
+(defun ttns-overlap-boundary-tensor (tensor1 tensor2)
+ (functional-tensor
+ (list (second (tensor-indices tensor2))
+ (conjugate-index (second (tensor-indices tensor1))))
+ #'(lambda (s)
+ (if (eql (subscript-subscript (first s))
+ (subscript-subscript (second s)))
+ 1
+ 0))))
+
+(defun ttns-boundary-tensor (tensor1 tensor2)
+ (functional-tensor
+ (list (second (tensor-indices tensor2))
+ (list (make-segment :numbers '(0) :dimension 1))
+ (conjugate-index (second (tensor-indices tensor1))))
+ #'(lambda (s)
+ (if (eql (subscript-subscript (first s))
+ (subscript-subscript (third s)))
+ 1
+ 0))))
+
+(defun ttns-overlap (ttns1 ttns2)
+ (loop
+ with w = (ttns-width ttns1)
+ with h = (ttns-height ttns1)
+ with l = (ttns-overlap-boundary-tensor (tensor-at-site ttns1 '(0 0)) (tensor-at-site ttns2 '(0 0)))
+ for i below w
+ for u = (loop
+ for j from h above 0
+ for v = (mps-overlap-boundary-tensor)
+ then (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i j)) '(1 1)))
+ (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i j)) '(1 1))))
+ (tensor-contract
+ (tensor-conjugate t2)
+ '(0 2)
+ (tensor-contract t1 2 v 1)
+ '(0 2)))
+ finally (return v))
+ do
+ (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i 0)) '(3)))
+ (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i 0)) '(3))))
+ (setf l (tensor-contract
+ (tensor-conjugate t2)
+ '(1 3 0)
+ (tensor-contract l 1 (tensor-contract u 1 t1 3) 2)
+ '(0 1 2))))
+ finally (return (tensor-scalar (tensor-contract l '(0 1) (mps-overlap-boundary-tensor) '(0 1))))))
+
+(defun ttns-expectation (ttns1 ttns2 ttno)
+ (loop
+ with w = (ttns-width ttns1)
+ with h = (ttns-height ttns1)
+ with l = (ttns-boundary-tensor (tensor-at-site ttns1 '(0 0)) (tensor-at-site ttns2 '(0 0)))
+ for i below w
+ for u = (loop
+ for j from h above 0
+ for v = (mps-boundary-tensor)
+ then (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i j)) '(1 1)))
+ (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i j)) '(1 1)))
+ (to (terminate-dummy-indices (tensor-at-site ttno (list i j)) '(1 1))))
+ (tensor-contract
+ (tensor-conjugate t2)
+ '(0 2)
+ (tensor-contract to '(2 3) (tensor-contract t1 2 v 2) '(3 0))
+ '(0 3)))
+ finally (return v))
+ do
+ (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i 0)) '(3)))
+ (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i 0)) '(3)))
+ (to (terminate-dummy-indices (tensor-at-site ttno (list i 0)) '(3))))
+ (setf l (tensor-contract
+ (tensor-conjugate t2)
+ '(0 1 3)
+ (tensor-contract
+ to
+ '(1 3 4)
+ (tensor-contract l 2 (tensor-contract u 2 t1 3) 3)
+ '(1 3 4))
+ '(0 2 3))))
+ finally (return (tensor-scalar (tensor-contract l '(0 1 2) (mps-boundary-tensor) '(0 1 2))))))
+
+(defun normalize-tensor-by-site (ttns site bond)
+ (let ((u (site-by-bond ttns site bond)))
+ (case bond
+ (:x- (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns site) '(0 2 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute l '(0 2 3 4 1)))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns u) 2)
+ '(2 0 1 3 4)))))
+ (:x+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns site) '(0 1 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute l '(0 1 3 4 2)))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns u) 1)
+ '(1 0 2 3 4)))))
+ (:y- (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns site) '(0 1 2 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute l '(0 1 2 4 3)))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns u) 4)
+ '(4 0 1 2 3)))))
+ (:y+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns site) '(0 1 2 3) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns site) l)
+ (setf (tensor-by-site ttns u)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns u) 3)
+ '(3 0 1 2 4))))))))
+
+(defun compress-ttns (ttns dimension)
+ (let ((w (ttns-width ttns))
+ (h (ttns-height ttns)))
+ ;; vertical top to bottom
+ (loop for i below w do
+ (loop for j from (1- h) above 0 do
+ (let* ((u (site-by-id ttns (list i j)))
+ (v (site-by-bond ttns u :y-)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns u) '(3) nil :right)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute r '(3 0 1 2 4)))
+ (setf (tensor-by-site ttns v)
+ (tensor-contract (tensor-by-site ttns v) 4 l 0))))))
+ ;; horizontal left to right
+ (loop for i below (1- w) do
+ (let* ((u (site-by-id ttns (list i 0)))
+ (v (site-by-bond ttns u :x+)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns u) '(0 1 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute l '(0 1 3 4 2)))
+ (setf (tensor-by-site ttns v)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns v) 1)
+ '(1 0 2 3 4))))))
+ ;; horizontal right to left
+ (loop for i from (1- w) downto 0 do
+ ;; vertical bottom to top
+ (loop for j from 0 below (1- h) do
+ (let* ((u (site-by-id ttns (list i j)))
+ (v (site-by-bond ttns u :y+)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns u) '(0 1 2 3) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns v) 3)
+ '(3 0 1 2 4))))))
+ ;; vertical top to bottom with truncation
+ (loop for j from (1- h) above 0 do
+ (let* ((u (site-by-id ttns (list i j)))
+ (v (site-by-bond ttns u :y-)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns u) '(3) dimension :right)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute r '(3 0 1 2 4)))
+ (setf (tensor-by-site ttns v)
+ (tensor-contract (tensor-by-site ttns v) 4 l 0)))))
+ ;; right to left with truncation
+ unless (zerop i) do
+ (let* ((u (site-by-id ttns (list i 0)))
+ (v (site-by-bond ttns u :x-)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns u) '(1) dimension :right)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute r '(1 0 2 3 4)))
+ (setf (tensor-by-site ttns v)
+ (tensor-permute
+ (tensor-contract (tensor-by-site ttns v) 2 l 0)
+ '(0 1 3 4 2))))))
+ (setf (tensor-by-site ttns (site-by-id ttns '(0 0)))
+ (tensor-normalize (tensor-by-site ttns (site-by-id ttns '(0 0)))))))
+
+(defun print-dimensions (ttns)
+ (let ((w (ttns-width ttns))
+ (h (ttns-height ttns)))
+ (loop
+ for j from (1- h) downto 0
+ do
+ (loop
+ for i below w
+ do (format t " ~3D " (fifth
+ (tensor-dimensions
+ (tensor-at-site
+ ttns
+ (list i j))))))
+ (terpri)
+ (loop
+ for i below w
+ for d = (tensor-dimensions (tensor-at-site ttns (list i j)))
+ do (format t "~3D * ~3D " (second d) (third d)))
+ (terpri)
+ (loop
+ for i below w
+ do (format t " ~3D " (fourth
+ (tensor-dimensions
+ (tensor-at-site
+ ttns
+ (list i j))))))
+ (terpri))))
+
+(defun split-index-tensor (tensor index)
+ (let ((i (nth index (tensor-indices tensor))))
+ (functional-tensor
+ (list (conjugate-index i)
+ i
+ (list (make-segment :numbers '(0) :dimension 1)))
+ #'(lambda (s)
+ (if (eql (subscript-subscript (first s))
+ (subscript-subscript (second s)))
+ 1
+ 0)))))
+
+(defun shift-right-bond-up (ttns site maxdim)
+ (let* ((u (tensor-by-site ttns site))
+ (v (tensor-contract u 2 (split-index-tensor u 2) 0)))
+ (multiple-value-bind (l s r) (tensor-decompose v '(3 4) maxdim :right)
+ (declare (ignore s))
+ (let ((w (site-by-bond ttns site :y+)))
+ (setf (tensor-by-site ttns w)
+ (tensor-fuse
+ (tensor-permute
+ (tensor-contract (tensor-by-site ttns w) 3 l 0)
+ '(0 1 2 5 3 4))
+ '(1 2 4 5 6)
+ '(:normal :normal :normal :normal :normal))))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute
+ r
+ '(4 0 1 3 2))))))
+
+(defun shift-left-bond-up (ttns site maxdim)
+ (let* ((u (tensor-by-site ttns site))
+ (v (tensor-contract u 1 (split-index-tensor u 1) 0)))
+ (multiple-value-bind (l s r) (tensor-decompose v '(3 4) maxdim :right)
+ (declare (ignore s))
+ (let ((w (site-by-bond ttns site :y+)))
+ (setf (tensor-by-site ttns w)
+ (tensor-fuse
+ (tensor-permute
+ (tensor-contract (tensor-by-site ttns w) 3 l 0)
+ '(0 1 3 5 2 4))
+ '(1 3 4 5 6)
+ '(:normal :reverse :normal :normal :normal))))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute
+ r
+ '(4 0 2 3 1))))))
+
+(defun shift-right-bond-down (ttns site maxdim)
+ (let* ((u (tensor-by-site ttns site))
+ (v (tensor-contract u 2 (split-index-tensor u 2) 0)))
+ (multiple-value-bind (l s r) (tensor-decompose v '(0 1 5 3) maxdim :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute
+ l
+ '(0 1 2 4 3)))
+ (let ((w (site-by-bond ttns site :y-)))
+ (setf (tensor-by-site ttns w)
+ (tensor-fuse
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns w) 4)
+ '(5 2 0 1 3 4))
+ '(1 2 3 5 6)
+ '(:normal :normal :normal :normal :normal)))))))
+
+(defun shift-left-bond-down (ttns site maxdim)
+ (let* ((u (tensor-by-site ttns site))
+ (v (tensor-contract u 1 (split-index-tensor u 1) 0)))
+ (multiple-value-bind (l s r) (tensor-decompose v '(0 5 1 3) maxdim :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns site)
+ (tensor-permute
+ l
+ '(0 1 2 4 3)))
+ (let ((w (site-by-bond ttns site :y-)))
+ (setf (tensor-by-site ttns w)
+ (tensor-fuse
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns w) 4)
+ '(5 1 0 2 3 4))
+ '(1 2 4 5 6)
+ '(:normal :reverse :normal :normal :normal)))))))
+
+;; time evolution
+
+(defun apply-single-site-operator-x (operator tensor)
+ (tensor-contract operator 1 tensor 0))
+
+(defun apply-double-site-operator-x (operator tensor1 tensor2 maxdim normalization)
+ (multiple-value-bind (l1 s1 r1) (tensor-decompose tensor1 '(1 3 4) nil :left)
+ (declare (ignore s1))
+ (multiple-value-bind (l2 s2 r2) (tensor-decompose tensor2 '(1 0) nil :right)
+ (declare (ignore s2))
+ (let ((u (tensor-contract
+ operator
+ '(1 3)
+ (tensor-contract r1 2 l2 0)
+ '(1 2))))
+ (multiple-value-bind (l s r)
+ (tensor-decompose u '(2 0) maxdim normalization)
+ (declare (ignore s))
+ (let ((v1 (tensor-permute
+ (tensor-contract l1 3 l 0)
+ '(1 3 4 0 2))))
+ (let ((v2 (tensor-permute
+ (tensor-contract r 2 r2 0)
+ '(1 0 2 3 4))))
+ (values v1 v2))))))))
+
+(defun apply-double-site-operator-y (operator tensor1 tensor2 maxdim normalization)
+ (multiple-value-bind (l1 s1 r1) (tensor-decompose tensor1 '(1 2 3) nil :left)
+ (declare (ignore s1))
+ (multiple-value-bind (l2 s2 r2) (tensor-decompose tensor2 '(3 0) nil :right)
+ (declare (ignore s2))
+ (let ((u (tensor-contract
+ operator
+ '(1 3)
+ (tensor-contract r1 2 l2 0)
+ '(1 2))))
+ (multiple-value-bind (l s r)
+ (tensor-decompose u '(2 0) maxdim normalization)
+ (declare (ignore s))
+ (values
+ (tensor-permute
+ (tensor-contract l1 3 l 0)
+ '(1 2 3 0 4))
+ (tensor-permute
+ (tensor-contract r 2 r2 0)
+ '(3 0 1 2 4))))))))
+
+(defun normalize-line (ttns site direction)
+ (loop
+ for u = site then v
+ for v = (site-by-bond ttns u direction)
+ while v
+ do
+ (normalize-tensor-by-site ttns u direction)))
+
+(defun apply-single-column-tebd-operators (ssps ttns x)
+ (loop
+ for u = (site-by-id ttns (list x 0)) then v
+ while u
+ for v = (site-by-bond ttns u :y+)
+ do
+ (setf (tensor-by-site ttns u)
+ (apply-single-site-operator-x (apply #'aref ssps (site-id u)) (tensor-by-site ttns u)))))
+
+(defun apply-double-column-tebd-operators (dsps ttns x maxdim normalization)
+ (loop
+ for u = (site-by-id ttns (list x 0)) then w
+ while u
+ for v = (site-by-bond ttns u :x+)
+ for w = (site-by-bond ttns u :y+)
+ do
+ (multiple-value-bind (l r)
+ (apply-double-site-operator-x
+ (apply #'aref dsps (site-id u))
+ (tensor-by-site ttns u)
+ (tensor-by-site ttns v)
+ maxdim
+ (if w
+ (case normalization
+ (:left :right)
+ (:right :left))
+ normalization))
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v) r))
+ when w
+ do
+ (case normalization
+ (:left
+ (shift-right-bond-up ttns u maxdim)
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns w) '(0 1 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns w)
+ (tensor-permute l '(0 1 3 4 2)))
+ (setf (tensor-by-site ttns v)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns v) 1)
+ '(1 0 2 3 4))))
+ (shift-left-bond-up ttns v maxdim))
+ (:right
+ (shift-left-bond-up ttns v maxdim)
+ (let ((z (site-by-bond ttns w :x+)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns z) '(0 2 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns z)
+ (tensor-permute l '(0 2 3 4 1)))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns u) 2)
+ '(2 0 1 3 4)))))
+ (shift-right-bond-up ttns u maxdim)))))
+
+(defun apply-double-column-tebd-operators-down (dsps ttns x maxdim normalization)
+ (loop
+ with h = (ttns-height ttns)
+ for u = (site-by-id ttns (list x (1- h))) then w
+ while u
+ for v = (site-by-bond ttns u :x+)
+ for w = (site-by-bond ttns u :y-)
+ do
+ (multiple-value-bind (l r)
+ (apply-double-site-operator-x
+ (apply #'aref dsps (site-id u))
+ (tensor-by-site ttns u)
+ (tensor-by-site ttns v)
+ maxdim
+ (if w
+ (case normalization
+ (:left :right)
+ (:right :left))
+ normalization))
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v) r))
+ when w
+ do
+ (case normalization
+ (:left
+ (shift-right-bond-down ttns u maxdim)
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns w) '(0 1 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns w)
+ (tensor-permute l '(0 1 3 4 2)))
+ (setf (tensor-by-site ttns v)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns v) 1)
+ '(1 0 2 3 4))))
+ (shift-left-bond-down ttns v maxdim))
+ (:right
+ (shift-left-bond-down ttns v maxdim)
+ (let ((z (site-by-bond ttns w :x+)))
+ (multiple-value-bind (l s r)
+ (tensor-decompose (tensor-by-site ttns z) '(0 2 3 4) nil :left)
+ (declare (ignore s))
+ (setf (tensor-by-site ttns z)
+ (tensor-permute l '(0 2 3 4 1)))
+ (setf (tensor-by-site ttns u)
+ (tensor-permute
+ (tensor-contract r 1 (tensor-by-site ttns u) 2)
+ '(2 0 1 3 4)))))
+ (shift-right-bond-down ttns u maxdim)))))
+
+(defun apply-even-tebd-operators-x (ssps dsps ttns maxdim)
+ (let ((w (ttns-width ttns))
+ (h (ttns-height ttns)))
+ (loop
+ for i below (1- w) by 2
+ do
+ (apply-double-column-tebd-operators dsps ttns i maxdim :left)
+ (unless (eql i (- w 2))
+ (normalize-line ttns (site-by-id ttns (list (1+ i) (1- h))) :y-)
+ (normalize-tensor-by-site ttns (site-by-id ttns (list (1+ i) 0)) :x+)))
+ (when (oddp w)
+ (apply-single-column-tebd-operators ssps ttns (1- w))
+ (normalize-line ttns (site-by-id ttns (list (1- w) (1- h))) :y-))))
+
+(defun apply-odd-tebd-operators-x (ssps dsps ttns maxdim)
+ (let ((w (ttns-width ttns))
+ (h (ttns-height ttns)))
+ (when (evenp w)
+ (apply-single-column-tebd-operators ssps ttns (1- w))
+ (normalize-line ttns (site-by-id ttns (list (1- w) 0)) :y+))
+ (loop
+ for i from (if (evenp w) (- w 3) (- w 2)) downto 1 by 2
+ do
+ (unless (eql i (- w 2))
+ (normalize-tensor-by-site ttns (site-by-id ttns (list (+ i 2) (1- h))) :x-)
+ (normalize-line ttns (site-by-id ttns (list (1+ i) (1- h))) :y-))
+ (apply-double-column-tebd-operators dsps ttns i maxdim :right)
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i (1- h))) :x-))))
+
+(defun apply-even-tebd-operators-x-down (ssps dsps ttns maxdim)
+ (let ((w (ttns-width ttns))
+ (h (ttns-height ttns)))
+ (when (oddp w)
+ (apply-single-column-tebd-operators ssps ttns (1- w))
+ (normalize-line ttns (site-by-id ttns (list (1- w) (1- h))) :y-))
+ (loop
+ for i from (if (evenp w) (- w 2) (- w 3)) downto 0 by 2
+ do
+ (unless (eql i (- w 2))
+ (normalize-tensor-by-site ttns (site-by-id ttns (list (+ i 2) 0)) :x-)
+ (normalize-line ttns (site-by-id ttns (list (1+ i) 0)) :y+))
+ (apply-double-column-tebd-operators-down dsps ttns i maxdim :right)
+ (unless (zerop i)
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x-)))))
+
+(defun apply-odd-tebd-operators-x-down (ssps dsps ttns maxdim)
+ (let ((w (ttns-width ttns))
+ (h (ttns-height ttns)))
+ (normalize-tensor-by-site ttns (site-by-id ttns (list 0 (1- h))) :x+)
+ (loop
+ for i from 1 below (1- w) by 2
+ do
+ (apply-double-column-tebd-operators-down dsps ttns i maxdim :left)
+ (unless (eql i (- w 2))
+ (normalize-line ttns (site-by-id ttns (list (1+ i) 0)) :y+)
+ (normalize-tensor-by-site ttns (site-by-id ttns (list (1+ i) (1- h))) :x+)))
+ (when (evenp w)
+ (apply-single-column-tebd-operators ssps ttns (1- w))
+ (normalize-line ttns (site-by-id ttns (list (1- w) 0)) :y+))))
+
+(defun apply-even-tebd-operators-y (double ttns maxdim)
+ (loop
+ with w = (ttns-width ttns)
+ with h = (ttns-height ttns)
+ for i below w
+ do
+ (loop
+ for j below (1- h) by 2
+ for u = (site-by-id ttns (list i j))
+ for v = (site-by-bond ttns u :y+)
+ do
+ (multiple-value-bind (l r)
+ (apply-double-site-operator-y
+ double
+ (tensor-by-site ttns u)
+ (tensor-by-site ttns v)
+ maxdim
+ :left)
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v) r))
+ (unless (>= j (- h 3))
+ (normalize-tensor-by-site ttns v :y+)))
+ (loop
+ for j from (1- h) above 0
+ do
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-))
+ (unless (eql i (1- w))
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x+))))
+
+(defun apply-odd-tebd-operators-y (double ttns maxdim)
+ (loop
+ with w = (ttns-width ttns)
+ with h = (ttns-height ttns)
+ for i from (1- w) downto 0
+ do
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :y+)
+ (loop
+ for j from 1 below (1- h) by 2
+ for u = (site-by-id ttns (list i j))
+ for v = (site-by-bond ttns u :y+)
+ do
+ (multiple-value-bind (l r)
+ (apply-double-site-operator-y
+ double
+ (tensor-by-site ttns u)
+ (tensor-by-site ttns v)
+ maxdim
+ :left)
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v) r))
+ (unless (>= j (- h 3))
+ (normalize-tensor-by-site ttns v :y+)))
+ (loop
+ for j from (1- h) above 0
+ do
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-))
+ (unless (zerop i)
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x-))))
+
+(defun apply-even-tebd-operators-y-down (double ttns maxdim)
+ (loop
+ with w = (ttns-width ttns)
+ with h = (ttns-height ttns)
+ for i from (1- w) downto 0
+ do
+ (loop
+ for j below (1- h) by 2
+ for u = (site-by-id ttns (list i j))
+ for v = (site-by-bond ttns u :y+)
+ do
+ (multiple-value-bind (l r)
+ (apply-double-site-operator-y
+ double
+ (tensor-by-site ttns u)
+ (tensor-by-site ttns v)
+ maxdim
+ :left)
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v) r))
+ (unless (>= j (- h 3))
+ (normalize-tensor-by-site ttns v :y+)))
+ (loop
+ for j from (1- h) above 0
+ do
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-))
+ (unless (zerop i)
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x-))
+ ))
+
+(defun apply-odd-tebd-operators-y-down (double ttns maxdim)
+ (loop
+ with w = (ttns-width ttns)
+ with h = (ttns-height ttns)
+ for i below w
+ do
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :y+)
+ (loop
+ for j from 1 below (1- h) by 2
+ for u = (site-by-id ttns (list i j))
+ for v = (site-by-bond ttns u :y+)
+ do
+ (multiple-value-bind (l r)
+ (apply-double-site-operator-y
+ double
+ (tensor-by-site ttns u)
+ (tensor-by-site ttns v)
+ maxdim
+ :left)
+ (setf (tensor-by-site ttns u) l)
+ (setf (tensor-by-site ttns v) r))
+ (unless (>= j (- h 3))
+ (normalize-tensor-by-site ttns v :y+)))
+ (loop
+ for j from (1- h) above 0
+ do
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-))
+ (unless (eql i (1- w))
+ (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x+))))
+
+(defun make-ssps (width height dimension interaction potential step)
+ (loop
+ with a = (make-array (list width height))
+ for i below width
+ do
+ (loop
+ for j below height
+ for h = (bose-hubbard-site-hamiltonian dimension 1 interaction (aref potential i j))
+ do
+ (setf (aref a i j) (make-single-site-propagator h step)))
+ finally (return a)))
+
+(defun make-dsps (width height dimension interaction potential step)
+ (loop
+ with a = (make-array (list (1- width) height))
+ for i below (1- width)
+ do
+ (loop
+ for j below height
+ for f = (bose-hubbard-site-hamiltonian dimension 1 interaction (aref potential i j))
+ for g = (bose-hubbard-site-hamiltonian dimension 1 interaction (aref potential (1+ i) j))
+ do
+ (setf (aref a i j) (make-double-site-propagator f g step)))
+ finally (return a)))
+
+(defun ttns-tebd-evolve (ttns ssps dsps idsp maxdim)
+ (apply-even-tebd-operators-y idsp ttns maxdim)
+ (apply-odd-tebd-operators-y idsp ttns maxdim)
+ (apply-even-tebd-operators-x ssps dsps ttns maxdim)
+ (apply-odd-tebd-operators-x ssps dsps ttns maxdim)
+ (apply-odd-tebd-operators-x-down ssps dsps ttns maxdim)
+ (apply-even-tebd-operators-x-down ssps dsps ttns maxdim)
+ (apply-odd-tebd-operators-y-down idsp ttns maxdim)
+ (apply-even-tebd-operators-y-down idsp ttns maxdim))
+
+(defun all-pn-operators (width height dimension)
+ (loop
+ with a = (make-array (list width height))
+ for i below width
+ do
+ (loop
+ for j below height
+ do
+ (setf (aref a i j)
+ (pn-operator width height i j dimension)))
+ finally
+ (return a)))
+
+(defun all-pn (ttns norm pn)
+ (loop
+ with w = (array-dimension pn 0)
+ with h = (array-dimension pn 1)
+ with a = (make-array (list w h))
+ for i below w
+ do
+ (loop
+ for j below h
+ do
+ (setf (aref a i j)
+ (realpart (/ (ttns-expectation
+ ttns
+ ttns
+ (aref pn i j))
+ norm))))
+ finally
+ (return a)))
+
+(defun print-pn (pn)
+ (loop
+ with w = (array-dimension pn 0)
+ with h = (array-dimension pn 1)
+ for j from (1- h) downto 0
+ do
+ (loop
+ for i below w
+ do
+ (format t
+ ;; "~,3f~:[ ~;~%~]"
+ "~F~:[ ~;~%~]"
+ (aref pn i j)
+ (eql i (1- w))))))
