ttns.lisp (31326B)
1 (defpackage :ttns 2 (:use :common-lisp :tensor :bhm :tebd) 3 (:export :make-vacuum-ttns 4 :ttns-create-particle 5 :ttns-overlap 6 :compress-ttns 7 :make-ssps 8 :make-dsps 9 :ttns-tebd-evolve 10 :all-pn-operators 11 :all-pn 12 :print-pn)) 13 14 (in-package :ttns) 15 16 (defstruct site 17 id 18 subscript 19 bonds) 20 21 (defstruct ttns 22 width 23 height 24 lattice 25 tensors) 26 27 (defun site-by-id (ttns id) 28 (find id (ttns-lattice ttns) :key #'site-id :test #'equal)) 29 30 (defun site-by-bond (ttns site bond) 31 (let ((i (cdr (assoc bond (site-bonds site))))) 32 (if i (svref (ttns-lattice ttns) i)))) 33 34 (defun tensor-by-site (ttns site) 35 (svref (ttns-tensors ttns) (site-subscript site))) 36 37 (defun (setf tensor-by-site) (tensor ttns site) 38 (setf 39 (svref (ttns-tensors ttns) (site-subscript site)) 40 tensor)) 41 42 (defun tensor-at-site (ttns id) 43 (svref (ttns-tensors ttns) (site-subscript (site-by-id ttns id)))) 44 45 (defun (setf tensor-at-site) (tensor ttns id) 46 (setf 47 (svref (ttns-tensors ttns) (site-subscript (site-by-id ttns id))) 48 tensor)) 49 50 (defun make-rectangular-lattice (width height) 51 (flet ((subscript (i j) 52 (if (and (>= i 0) (< i width) (>= j 0) (< j height)) 53 (+ i (* j width))))) 54 (let ((a (make-array (* width height)))) 55 (dotimes (i width a) 56 (dotimes (j height) 57 (let ((s (subscript i j))) 58 (setf (svref a s) 59 (make-site 60 :id (list i j) 61 :subscript s 62 :bonds (list 63 (cons :x- (subscript (1- i) j)) 64 (cons :x+ (subscript (1+ i) j)) 65 (cons :y- (subscript i (1- j))) 66 (cons :y+ (subscript i (1+ j)))))))))))) 67 68 (defun make-vacuum-ttns (width height dimension) 69 (let* ((l (make-rectangular-lattice width height)) 70 (a (make-array (length l)))) 71 (dotimes (i (length l) 72 (make-ttns :width width :height height :lattice l :tensors a)) 73 (setf (svref a i) 74 (functional-tensor 75 (list 76 (ket-physical-index dimension) 77 (list (make-segment :numbers '(0) :dimension 1)) 78 (list (make-segment :numbers '(0) :dimension 1)) 79 (list (make-segment :numbers '(0) :dimension 1)) 80 (list (make-segment :numbers '(0) :dimension 1))) 81 #'(lambda (s) (declare (ignore s)) 1)))))) 82 83 (defun identity-site-operator (dimension) 84 (functional-tensor 85 (list 86 (ket-physical-index dimension) 87 (list (make-segment :numbers '(0) :dimension 1)) 88 (list (make-segment :numbers '(0) :dimension 1)) 89 (list (make-segment :numbers '(0) :dimension 1)) 90 (list (make-segment :numbers '(0) :dimension 1)) 91 (bra-physical-index dimension)) 92 #'(lambda (s) 93 (declare (ignore s)) 1))) 94 95 (defun xx-connecting-site-operator (dimension) 96 (functional-tensor 97 (list 98 (ket-physical-index dimension) 99 (list (make-segment :numbers '(-1) :dimension 1)) 100 (list (make-segment :numbers '( 1) :dimension 1)) 101 (list (make-segment :numbers '( 0) :dimension 1)) 102 (list (make-segment :numbers '( 0) :dimension 1)) 103 (bra-physical-index dimension)) 104 #'(lambda (s) (declare (ignore s)) 1))) 105 106 (defun yy-connecting-site-operator (dimension) 107 (functional-tensor 108 (list 109 (ket-physical-index dimension) 110 (list (make-segment :numbers '( 0) :dimension 1)) 111 (list (make-segment :numbers '( 0) :dimension 1)) 112 (list (make-segment :numbers '(-1) :dimension 1)) 113 (list (make-segment :numbers '( 1) :dimension 1)) 114 (bra-physical-index dimension)) 115 #'(lambda (s) (declare (ignore s)) 1))) 116 117 (defun xy-connecting-site-operator (dimension) 118 (functional-tensor 119 (list 120 (ket-physical-index dimension) 121 (list (make-segment :numbers '(-1) :dimension 1)) 122 (list (make-segment :numbers '( 0) :dimension 1)) 123 (list (make-segment :numbers '( 0) :dimension 1)) 124 (list (make-segment :numbers '( 1) :dimension 1)) 125 (bra-physical-index dimension)) 126 #'(lambda (s) (declare (ignore s)) 1))) 127 128 (defun x-creation-site-operator (dimension) 129 (functional-tensor 130 (list 131 (ket-physical-index dimension) 132 (list (make-segment :numbers '(-1) :dimension 1)) 133 (list (make-segment :numbers '( 0) :dimension 1)) 134 (list (make-segment :numbers '( 0) :dimension 1)) 135 (list (make-segment :numbers '( 0) :dimension 1)) 136 (bra-physical-index dimension)) 137 #'(lambda (s) 138 (let ((n (- (first (subscript-numbers (sixth s)))))) 139 (sqrt (1+ n)))))) 140 141 (defun y-creation-site-operator (dimension) 142 (functional-tensor 143 (list 144 (ket-physical-index dimension) 145 (list (make-segment :numbers '( 0) :dimension 1)) 146 (list (make-segment :numbers '( 0) :dimension 1)) 147 (list (make-segment :numbers '(-1) :dimension 1)) 148 (list (make-segment :numbers '( 0) :dimension 1)) 149 (bra-physical-index dimension)) 150 #'(lambda (s) 151 (let ((n (- (first (subscript-numbers (sixth s)))))) 152 (sqrt (1+ n)))))) 153 154 (defun pn-site-operator (dimension) 155 (functional-tensor 156 (list 157 (ket-physical-index dimension) 158 (list (make-segment :numbers '(0) :dimension 1)) 159 (list (make-segment :numbers '(0) :dimension 1)) 160 (list (make-segment :numbers '(0) :dimension 1)) 161 (list (make-segment :numbers '(0) :dimension 1)) 162 (bra-physical-index dimension)) 163 #'(lambda (s) 164 (- (first (subscript-numbers (sixth s))))))) 165 166 (defun identity-operator (width height dimension) 167 (let* ((l (make-rectangular-lattice width height)) 168 (a (make-array (length l))) 169 (o (identity-site-operator dimension))) 170 (dotimes (i (length l) 171 (make-ttns 172 :width width 173 :height height 174 :lattice l 175 :tensors a)) 176 (setf (svref a i) o)))) 177 178 (defun apply-and-fuse (operator ttns) 179 (let* ((l (length (ttns-lattice ttns))) 180 (o (ttns-tensors operator)) 181 (p (ttns-tensors ttns)) 182 (a (make-array l))) 183 (dotimes (i l 184 (make-ttns 185 :width (ttns-width ttns) 186 :height (ttns-height ttns) 187 :lattice (ttns-lattice ttns) 188 :tensors a)) 189 (setf (svref a i) 190 (tensor-fuse 191 (tensor-permute 192 (tensor-contract (svref o i) 5 (svref p i) 0) 193 '(0 1 3 5 7 2 4 6 8)) 194 '(1 3 5 7) 195 '(:normal :reverse :normal :reverse :normal)))))) 196 197 (defun ttns-create-particle (ttns x y dimension) 198 (flet ((f (i j o) 199 (let ((s (site-subscript (site-by-id ttns (list i j)))) 200 (p (ttns-tensors ttns))) 201 (setf (svref p s) 202 (tensor-fuse 203 (tensor-permute 204 (tensor-contract o 5 (svref p s) 0) 205 '(0 1 3 5 7 2 4 6 8)) 206 '(1 3 5 7) 207 '(:normal :reverse :normal :reverse :normal)))))) 208 (loop 209 with o = (xx-connecting-site-operator dimension) 210 for i below x 211 do (f i 0 o)) 212 (loop 213 with o = (yy-connecting-site-operator dimension) 214 for i from 1 below y 215 do (f x i o)) 216 (if (zerop y) 217 (f x y (x-creation-site-operator dimension)) 218 (progn 219 (f x 0 (xy-connecting-site-operator dimension)) 220 (f x y (y-creation-site-operator dimension)))))) 221 222 (defun pn-operator (width height x y dimension) 223 (let ((o (identity-operator width height dimension))) 224 (setf (tensor-at-site o (list x y)) 225 (pn-site-operator dimension)) 226 o)) 227 228 (defun mps-boundary-tensor () 229 (let ((m (list (make-segment :numbers '(0) :dimension 1)))) 230 (functional-tensor 231 (list m m m) 232 #'(lambda (s) (declare (ignore s)) 1)))) 233 234 (defun mps-overlap-boundary-tensor () 235 (let ((m (list (make-segment :numbers '(0) :dimension 1)))) 236 (functional-tensor 237 (list m m) 238 #'(lambda (s) (declare (ignore s)) 1)))) 239 240 (defun terminate-dummy-index (tensor index) 241 (tensor-contract 242 tensor 243 index 244 (functional-tensor 245 (list (list (make-segment :numbers '(0) :dimension 1))) 246 #'(lambda (s) (declare (ignore s)) 1)) 247 0)) 248 249 (defun terminate-dummy-indices (tensor indices) 250 (reduce #'terminate-dummy-index indices :initial-value tensor)) 251 252 (defun ttns-overlap-boundary-tensor (tensor1 tensor2) 253 (functional-tensor 254 (list (second (tensor-indices tensor2)) 255 (conjugate-index (second (tensor-indices tensor1)))) 256 #'(lambda (s) 257 (if (eql (subscript-subscript (first s)) 258 (subscript-subscript (second s))) 259 1 260 0)))) 261 262 (defun ttns-boundary-tensor (tensor1 tensor2) 263 (functional-tensor 264 (list (second (tensor-indices tensor2)) 265 (list (make-segment :numbers '(0) :dimension 1)) 266 (conjugate-index (second (tensor-indices tensor1)))) 267 #'(lambda (s) 268 (if (eql (subscript-subscript (first s)) 269 (subscript-subscript (third s))) 270 1 271 0)))) 272 273 (defun ttns-overlap (ttns1 ttns2) 274 (loop 275 with w = (ttns-width ttns1) 276 with h = (ttns-height ttns1) 277 with l = (ttns-overlap-boundary-tensor (tensor-at-site ttns1 '(0 0)) (tensor-at-site ttns2 '(0 0))) 278 for i below w 279 for u = (loop 280 for j from h above 0 281 for v = (mps-overlap-boundary-tensor) 282 then (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i j)) '(1 1))) 283 (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i j)) '(1 1)))) 284 (tensor-contract 285 (tensor-conjugate t2) 286 '(0 2) 287 (tensor-contract t1 2 v 1) 288 '(0 2))) 289 finally (return v)) 290 do 291 (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i 0)) '(3))) 292 (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i 0)) '(3)))) 293 (setf l (tensor-contract 294 (tensor-conjugate t2) 295 '(1 3 0) 296 (tensor-contract l 1 (tensor-contract u 1 t1 3) 2) 297 '(0 1 2)))) 298 finally (return (tensor-scalar (tensor-contract l '(0 1) (mps-overlap-boundary-tensor) '(0 1)))))) 299 300 (defun ttns-expectation (ttns1 ttns2 ttno) 301 (loop 302 with w = (ttns-width ttns1) 303 with h = (ttns-height ttns1) 304 with l = (ttns-boundary-tensor (tensor-at-site ttns1 '(0 0)) (tensor-at-site ttns2 '(0 0))) 305 for i below w 306 for u = (loop 307 for j from h above 0 308 for v = (mps-boundary-tensor) 309 then (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i j)) '(1 1))) 310 (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i j)) '(1 1))) 311 (to (terminate-dummy-indices (tensor-at-site ttno (list i j)) '(1 1)))) 312 (tensor-contract 313 (tensor-conjugate t2) 314 '(0 2) 315 (tensor-contract to '(2 3) (tensor-contract t1 2 v 2) '(3 0)) 316 '(0 3))) 317 finally (return v)) 318 do 319 (let ((t1 (terminate-dummy-indices (tensor-at-site ttns1 (list i 0)) '(3))) 320 (t2 (terminate-dummy-indices (tensor-at-site ttns2 (list i 0)) '(3))) 321 (to (terminate-dummy-indices (tensor-at-site ttno (list i 0)) '(3)))) 322 (setf l (tensor-contract 323 (tensor-conjugate t2) 324 '(0 1 3) 325 (tensor-contract 326 to 327 '(1 3 4) 328 (tensor-contract l 2 (tensor-contract u 2 t1 3) 3) 329 '(1 3 4)) 330 '(0 2 3)))) 331 finally (return (tensor-scalar (tensor-contract l '(0 1 2) (mps-boundary-tensor) '(0 1 2)))))) 332 333 (defun normalize-tensor-by-site (ttns site bond) 334 (let ((u (site-by-bond ttns site bond))) 335 (case bond 336 (:x- (multiple-value-bind (l s r) 337 (tensor-decompose (tensor-by-site ttns site) '(0 2 3 4) nil :left) 338 (declare (ignore s)) 339 (setf (tensor-by-site ttns site) 340 (tensor-permute l '(0 2 3 4 1))) 341 (setf (tensor-by-site ttns u) 342 (tensor-permute 343 (tensor-contract r 1 (tensor-by-site ttns u) 2) 344 '(2 0 1 3 4))))) 345 (:x+ (multiple-value-bind (l s r) 346 (tensor-decompose (tensor-by-site ttns site) '(0 1 3 4) nil :left) 347 (declare (ignore s)) 348 (setf (tensor-by-site ttns site) 349 (tensor-permute l '(0 1 3 4 2))) 350 (setf (tensor-by-site ttns u) 351 (tensor-permute 352 (tensor-contract r 1 (tensor-by-site ttns u) 1) 353 '(1 0 2 3 4))))) 354 (:y- (multiple-value-bind (l s r) 355 (tensor-decompose (tensor-by-site ttns site) '(0 1 2 4) nil :left) 356 (declare (ignore s)) 357 (setf (tensor-by-site ttns site) 358 (tensor-permute l '(0 1 2 4 3))) 359 (setf (tensor-by-site ttns u) 360 (tensor-permute 361 (tensor-contract r 1 (tensor-by-site ttns u) 4) 362 '(4 0 1 2 3))))) 363 (:y+ (multiple-value-bind (l s r) 364 (tensor-decompose (tensor-by-site ttns site) '(0 1 2 3) nil :left) 365 (declare (ignore s)) 366 (setf (tensor-by-site ttns site) l) 367 (setf (tensor-by-site ttns u) 368 (tensor-permute 369 (tensor-contract r 1 (tensor-by-site ttns u) 3) 370 '(3 0 1 2 4)))))))) 371 372 (defun compress-ttns (ttns dimension) 373 (let ((w (ttns-width ttns)) 374 (h (ttns-height ttns))) 375 ;; vertical top to bottom 376 (loop for i below w do 377 (loop for j from (1- h) above 0 do 378 (let* ((u (site-by-id ttns (list i j))) 379 (v (site-by-bond ttns u :y-))) 380 (multiple-value-bind (l s r) 381 (tensor-decompose (tensor-by-site ttns u) '(3) nil :right) 382 (declare (ignore s)) 383 (setf (tensor-by-site ttns u) 384 (tensor-permute r '(3 0 1 2 4))) 385 (setf (tensor-by-site ttns v) 386 (tensor-contract (tensor-by-site ttns v) 4 l 0)))))) 387 ;; horizontal left to right 388 (loop for i below (1- w) do 389 (let* ((u (site-by-id ttns (list i 0))) 390 (v (site-by-bond ttns u :x+))) 391 (multiple-value-bind (l s r) 392 (tensor-decompose (tensor-by-site ttns u) '(0 1 3 4) nil :left) 393 (declare (ignore s)) 394 (setf (tensor-by-site ttns u) 395 (tensor-permute l '(0 1 3 4 2))) 396 (setf (tensor-by-site ttns v) 397 (tensor-permute 398 (tensor-contract r 1 (tensor-by-site ttns v) 1) 399 '(1 0 2 3 4)))))) 400 ;; horizontal right to left 401 (loop for i from (1- w) downto 0 do 402 ;; vertical bottom to top 403 (loop for j from 0 below (1- h) do 404 (let* ((u (site-by-id ttns (list i j))) 405 (v (site-by-bond ttns u :y+))) 406 (multiple-value-bind (l s r) 407 (tensor-decompose (tensor-by-site ttns u) '(0 1 2 3) nil :left) 408 (declare (ignore s)) 409 (setf (tensor-by-site ttns u) l) 410 (setf (tensor-by-site ttns v) 411 (tensor-permute 412 (tensor-contract r 1 (tensor-by-site ttns v) 3) 413 '(3 0 1 2 4)))))) 414 ;; vertical top to bottom with truncation 415 (loop for j from (1- h) above 0 do 416 (let* ((u (site-by-id ttns (list i j))) 417 (v (site-by-bond ttns u :y-))) 418 (multiple-value-bind (l s r) 419 (tensor-decompose (tensor-by-site ttns u) '(3) dimension :right) 420 (declare (ignore s)) 421 (setf (tensor-by-site ttns u) 422 (tensor-permute r '(3 0 1 2 4))) 423 (setf (tensor-by-site ttns v) 424 (tensor-contract (tensor-by-site ttns v) 4 l 0))))) 425 ;; right to left with truncation 426 unless (zerop i) do 427 (let* ((u (site-by-id ttns (list i 0))) 428 (v (site-by-bond ttns u :x-))) 429 (multiple-value-bind (l s r) 430 (tensor-decompose (tensor-by-site ttns u) '(1) dimension :right) 431 (declare (ignore s)) 432 (setf (tensor-by-site ttns u) 433 (tensor-permute r '(1 0 2 3 4))) 434 (setf (tensor-by-site ttns v) 435 (tensor-permute 436 (tensor-contract (tensor-by-site ttns v) 2 l 0) 437 '(0 1 3 4 2)))))) 438 (setf (tensor-by-site ttns (site-by-id ttns '(0 0))) 439 (tensor-normalize (tensor-by-site ttns (site-by-id ttns '(0 0))))))) 440 441 (defun print-dimensions (ttns) 442 (let ((w (ttns-width ttns)) 443 (h (ttns-height ttns))) 444 (loop 445 for j from (1- h) downto 0 446 do 447 (loop 448 for i below w 449 do (format t " ~3D " (fifth 450 (tensor-dimensions 451 (tensor-at-site 452 ttns 453 (list i j)))))) 454 (terpri) 455 (loop 456 for i below w 457 for d = (tensor-dimensions (tensor-at-site ttns (list i j))) 458 do (format t "~3D * ~3D " (second d) (third d))) 459 (terpri) 460 (loop 461 for i below w 462 do (format t " ~3D " (fourth 463 (tensor-dimensions 464 (tensor-at-site 465 ttns 466 (list i j)))))) 467 (terpri)))) 468 469 (defun split-index-tensor (tensor index) 470 (let ((i (nth index (tensor-indices tensor)))) 471 (functional-tensor 472 (list (conjugate-index i) 473 i 474 (list (make-segment :numbers '(0) :dimension 1))) 475 #'(lambda (s) 476 (if (eql (subscript-subscript (first s)) 477 (subscript-subscript (second s))) 478 1 479 0))))) 480 481 (defun shift-right-bond-up (ttns site maxdim) 482 (let* ((u (tensor-by-site ttns site)) 483 (v (tensor-contract u 2 (split-index-tensor u 2) 0))) 484 (multiple-value-bind (l s r) (tensor-decompose v '(3 4) maxdim :right) 485 (declare (ignore s)) 486 (let ((w (site-by-bond ttns site :y+))) 487 (setf (tensor-by-site ttns w) 488 (tensor-fuse 489 (tensor-permute 490 (tensor-contract (tensor-by-site ttns w) 3 l 0) 491 '(0 1 2 5 3 4)) 492 '(1 2 4 5 6) 493 '(:normal :normal :normal :normal :normal)))) 494 (setf (tensor-by-site ttns site) 495 (tensor-permute 496 r 497 '(4 0 1 3 2)))))) 498 499 (defun shift-left-bond-up (ttns site maxdim) 500 (let* ((u (tensor-by-site ttns site)) 501 (v (tensor-contract u 1 (split-index-tensor u 1) 0))) 502 (multiple-value-bind (l s r) (tensor-decompose v '(3 4) maxdim :right) 503 (declare (ignore s)) 504 (let ((w (site-by-bond ttns site :y+))) 505 (setf (tensor-by-site ttns w) 506 (tensor-fuse 507 (tensor-permute 508 (tensor-contract (tensor-by-site ttns w) 3 l 0) 509 '(0 1 3 5 2 4)) 510 '(1 3 4 5 6) 511 '(:normal :reverse :normal :normal :normal)))) 512 (setf (tensor-by-site ttns site) 513 (tensor-permute 514 r 515 '(4 0 2 3 1)))))) 516 517 (defun shift-right-bond-down (ttns site maxdim) 518 (let* ((u (tensor-by-site ttns site)) 519 (v (tensor-contract u 2 (split-index-tensor u 2) 0))) 520 (multiple-value-bind (l s r) (tensor-decompose v '(0 1 5 3) maxdim :left) 521 (declare (ignore s)) 522 (setf (tensor-by-site ttns site) 523 (tensor-permute 524 l 525 '(0 1 2 4 3))) 526 (let ((w (site-by-bond ttns site :y-))) 527 (setf (tensor-by-site ttns w) 528 (tensor-fuse 529 (tensor-permute 530 (tensor-contract r 1 (tensor-by-site ttns w) 4) 531 '(5 2 0 1 3 4)) 532 '(1 2 3 5 6) 533 '(:normal :normal :normal :normal :normal))))))) 534 535 (defun shift-left-bond-down (ttns site maxdim) 536 (let* ((u (tensor-by-site ttns site)) 537 (v (tensor-contract u 1 (split-index-tensor u 1) 0))) 538 (multiple-value-bind (l s r) (tensor-decompose v '(0 5 1 3) maxdim :left) 539 (declare (ignore s)) 540 (setf (tensor-by-site ttns site) 541 (tensor-permute 542 l 543 '(0 1 2 4 3))) 544 (let ((w (site-by-bond ttns site :y-))) 545 (setf (tensor-by-site ttns w) 546 (tensor-fuse 547 (tensor-permute 548 (tensor-contract r 1 (tensor-by-site ttns w) 4) 549 '(5 1 0 2 3 4)) 550 '(1 2 4 5 6) 551 '(:normal :reverse :normal :normal :normal))))))) 552 553 ;; time evolution 554 555 (defun apply-single-site-operator-x (operator tensor) 556 (tensor-contract operator 1 tensor 0)) 557 558 (defun apply-double-site-operator-x (operator tensor1 tensor2 maxdim normalization) 559 (multiple-value-bind (l1 s1 r1) (tensor-decompose tensor1 '(1 3 4) nil :left) 560 (declare (ignore s1)) 561 (multiple-value-bind (l2 s2 r2) (tensor-decompose tensor2 '(1 0) nil :right) 562 (declare (ignore s2)) 563 (let ((u (tensor-contract 564 operator 565 '(1 3) 566 (tensor-contract r1 2 l2 0) 567 '(1 2)))) 568 (multiple-value-bind (l s r) 569 (tensor-decompose u '(2 0) maxdim normalization) 570 (declare (ignore s)) 571 (let ((v1 (tensor-permute 572 (tensor-contract l1 3 l 0) 573 '(1 3 4 0 2)))) 574 (let ((v2 (tensor-permute 575 (tensor-contract r 2 r2 0) 576 '(1 0 2 3 4)))) 577 (values v1 v2)))))))) 578 579 (defun apply-double-site-operator-y (operator tensor1 tensor2 maxdim normalization) 580 (multiple-value-bind (l1 s1 r1) (tensor-decompose tensor1 '(1 2 3) nil :left) 581 (declare (ignore s1)) 582 (multiple-value-bind (l2 s2 r2) (tensor-decompose tensor2 '(3 0) nil :right) 583 (declare (ignore s2)) 584 (let ((u (tensor-contract 585 operator 586 '(1 3) 587 (tensor-contract r1 2 l2 0) 588 '(1 2)))) 589 (multiple-value-bind (l s r) 590 (tensor-decompose u '(2 0) maxdim normalization) 591 (declare (ignore s)) 592 (values 593 (tensor-permute 594 (tensor-contract l1 3 l 0) 595 '(1 2 3 0 4)) 596 (tensor-permute 597 (tensor-contract r 2 r2 0) 598 '(3 0 1 2 4)))))))) 599 600 (defun normalize-line (ttns site direction) 601 (loop 602 for u = site then v 603 for v = (site-by-bond ttns u direction) 604 while v 605 do 606 (normalize-tensor-by-site ttns u direction))) 607 608 (defun apply-single-column-tebd-operators (ssps ttns x) 609 (loop 610 for u = (site-by-id ttns (list x 0)) then v 611 while u 612 for v = (site-by-bond ttns u :y+) 613 do 614 (setf (tensor-by-site ttns u) 615 (apply-single-site-operator-x (apply #'aref ssps (site-id u)) (tensor-by-site ttns u))))) 616 617 (defun apply-double-column-tebd-operators (dsps ttns x maxdim normalization) 618 (loop 619 for u = (site-by-id ttns (list x 0)) then w 620 while u 621 for v = (site-by-bond ttns u :x+) 622 for w = (site-by-bond ttns u :y+) 623 do 624 (multiple-value-bind (l r) 625 (apply-double-site-operator-x 626 (apply #'aref dsps (site-id u)) 627 (tensor-by-site ttns u) 628 (tensor-by-site ttns v) 629 maxdim 630 (if w 631 (case normalization 632 (:left :right) 633 (:right :left)) 634 normalization)) 635 (setf (tensor-by-site ttns u) l) 636 (setf (tensor-by-site ttns v) r)) 637 when w 638 do 639 (case normalization 640 (:left 641 (shift-right-bond-up ttns u maxdim) 642 (multiple-value-bind (l s r) 643 (tensor-decompose (tensor-by-site ttns w) '(0 1 3 4) nil :left) 644 (declare (ignore s)) 645 (setf (tensor-by-site ttns w) 646 (tensor-permute l '(0 1 3 4 2))) 647 (setf (tensor-by-site ttns v) 648 (tensor-permute 649 (tensor-contract r 1 (tensor-by-site ttns v) 1) 650 '(1 0 2 3 4)))) 651 (shift-left-bond-up ttns v maxdim)) 652 (:right 653 (shift-left-bond-up ttns v maxdim) 654 (let ((z (site-by-bond ttns w :x+))) 655 (multiple-value-bind (l s r) 656 (tensor-decompose (tensor-by-site ttns z) '(0 2 3 4) nil :left) 657 (declare (ignore s)) 658 (setf (tensor-by-site ttns z) 659 (tensor-permute l '(0 2 3 4 1))) 660 (setf (tensor-by-site ttns u) 661 (tensor-permute 662 (tensor-contract r 1 (tensor-by-site ttns u) 2) 663 '(2 0 1 3 4))))) 664 (shift-right-bond-up ttns u maxdim))))) 665 666 (defun apply-double-column-tebd-operators-down (dsps ttns x maxdim normalization) 667 (loop 668 with h = (ttns-height ttns) 669 for u = (site-by-id ttns (list x (1- h))) then w 670 while u 671 for v = (site-by-bond ttns u :x+) 672 for w = (site-by-bond ttns u :y-) 673 do 674 (multiple-value-bind (l r) 675 (apply-double-site-operator-x 676 (apply #'aref dsps (site-id u)) 677 (tensor-by-site ttns u) 678 (tensor-by-site ttns v) 679 maxdim 680 (if w 681 (case normalization 682 (:left :right) 683 (:right :left)) 684 normalization)) 685 (setf (tensor-by-site ttns u) l) 686 (setf (tensor-by-site ttns v) r)) 687 when w 688 do 689 (case normalization 690 (:left 691 (shift-right-bond-down ttns u maxdim) 692 (multiple-value-bind (l s r) 693 (tensor-decompose (tensor-by-site ttns w) '(0 1 3 4) nil :left) 694 (declare (ignore s)) 695 (setf (tensor-by-site ttns w) 696 (tensor-permute l '(0 1 3 4 2))) 697 (setf (tensor-by-site ttns v) 698 (tensor-permute 699 (tensor-contract r 1 (tensor-by-site ttns v) 1) 700 '(1 0 2 3 4)))) 701 (shift-left-bond-down ttns v maxdim)) 702 (:right 703 (shift-left-bond-down ttns v maxdim) 704 (let ((z (site-by-bond ttns w :x+))) 705 (multiple-value-bind (l s r) 706 (tensor-decompose (tensor-by-site ttns z) '(0 2 3 4) nil :left) 707 (declare (ignore s)) 708 (setf (tensor-by-site ttns z) 709 (tensor-permute l '(0 2 3 4 1))) 710 (setf (tensor-by-site ttns u) 711 (tensor-permute 712 (tensor-contract r 1 (tensor-by-site ttns u) 2) 713 '(2 0 1 3 4))))) 714 (shift-right-bond-down ttns u maxdim))))) 715 716 (defun apply-even-tebd-operators-x (ssps dsps ttns maxdim) 717 (let ((w (ttns-width ttns)) 718 (h (ttns-height ttns))) 719 (loop 720 for i below (1- w) by 2 721 do 722 (apply-double-column-tebd-operators dsps ttns i maxdim :left) 723 (unless (eql i (- w 2)) 724 (normalize-line ttns (site-by-id ttns (list (1+ i) (1- h))) :y-) 725 (normalize-tensor-by-site ttns (site-by-id ttns (list (1+ i) 0)) :x+))) 726 (when (oddp w) 727 (apply-single-column-tebd-operators ssps ttns (1- w)) 728 (normalize-line ttns (site-by-id ttns (list (1- w) (1- h))) :y-)))) 729 730 (defun apply-odd-tebd-operators-x (ssps dsps ttns maxdim) 731 (let ((w (ttns-width ttns)) 732 (h (ttns-height ttns))) 733 (when (evenp w) 734 (apply-single-column-tebd-operators ssps ttns (1- w)) 735 (normalize-line ttns (site-by-id ttns (list (1- w) 0)) :y+)) 736 (loop 737 for i from (if (evenp w) (- w 3) (- w 2)) downto 1 by 2 738 do 739 (unless (eql i (- w 2)) 740 (normalize-tensor-by-site ttns (site-by-id ttns (list (+ i 2) (1- h))) :x-) 741 (normalize-line ttns (site-by-id ttns (list (1+ i) (1- h))) :y-)) 742 (apply-double-column-tebd-operators dsps ttns i maxdim :right) 743 (normalize-tensor-by-site ttns (site-by-id ttns (list i (1- h))) :x-)))) 744 745 (defun apply-even-tebd-operators-x-down (ssps dsps ttns maxdim) 746 (let ((w (ttns-width ttns)) 747 (h (ttns-height ttns))) 748 (when (oddp w) 749 (apply-single-column-tebd-operators ssps ttns (1- w)) 750 (normalize-line ttns (site-by-id ttns (list (1- w) (1- h))) :y-)) 751 (loop 752 for i from (if (evenp w) (- w 2) (- w 3)) downto 0 by 2 753 do 754 (unless (eql i (- w 2)) 755 (normalize-tensor-by-site ttns (site-by-id ttns (list (+ i 2) 0)) :x-) 756 (normalize-line ttns (site-by-id ttns (list (1+ i) 0)) :y+)) 757 (apply-double-column-tebd-operators-down dsps ttns i maxdim :right) 758 (unless (zerop i) 759 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x-))))) 760 761 (defun apply-odd-tebd-operators-x-down (ssps dsps ttns maxdim) 762 (let ((w (ttns-width ttns)) 763 (h (ttns-height ttns))) 764 (normalize-tensor-by-site ttns (site-by-id ttns (list 0 (1- h))) :x+) 765 (loop 766 for i from 1 below (1- w) by 2 767 do 768 (apply-double-column-tebd-operators-down dsps ttns i maxdim :left) 769 (unless (eql i (- w 2)) 770 (normalize-line ttns (site-by-id ttns (list (1+ i) 0)) :y+) 771 (normalize-tensor-by-site ttns (site-by-id ttns (list (1+ i) (1- h))) :x+))) 772 (when (evenp w) 773 (apply-single-column-tebd-operators ssps ttns (1- w)) 774 (normalize-line ttns (site-by-id ttns (list (1- w) 0)) :y+)))) 775 776 (defun apply-even-tebd-operators-y (double ttns maxdim) 777 (loop 778 with w = (ttns-width ttns) 779 with h = (ttns-height ttns) 780 for i below w 781 do 782 (loop 783 for j below (1- h) by 2 784 for u = (site-by-id ttns (list i j)) 785 for v = (site-by-bond ttns u :y+) 786 do 787 (multiple-value-bind (l r) 788 (apply-double-site-operator-y 789 double 790 (tensor-by-site ttns u) 791 (tensor-by-site ttns v) 792 maxdim 793 :left) 794 (setf (tensor-by-site ttns u) l) 795 (setf (tensor-by-site ttns v) r)) 796 (unless (>= j (- h 3)) 797 (normalize-tensor-by-site ttns v :y+))) 798 (loop 799 for j from (1- h) above 0 800 do 801 (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-)) 802 (unless (eql i (1- w)) 803 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x+)))) 804 805 (defun apply-odd-tebd-operators-y (double ttns maxdim) 806 (loop 807 with w = (ttns-width ttns) 808 with h = (ttns-height ttns) 809 for i from (1- w) downto 0 810 do 811 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :y+) 812 (loop 813 for j from 1 below (1- h) by 2 814 for u = (site-by-id ttns (list i j)) 815 for v = (site-by-bond ttns u :y+) 816 do 817 (multiple-value-bind (l r) 818 (apply-double-site-operator-y 819 double 820 (tensor-by-site ttns u) 821 (tensor-by-site ttns v) 822 maxdim 823 :left) 824 (setf (tensor-by-site ttns u) l) 825 (setf (tensor-by-site ttns v) r)) 826 (unless (>= j (- h 3)) 827 (normalize-tensor-by-site ttns v :y+))) 828 (loop 829 for j from (1- h) above 0 830 do 831 (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-)) 832 (unless (zerop i) 833 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x-)))) 834 835 (defun apply-even-tebd-operators-y-down (double ttns maxdim) 836 (loop 837 with w = (ttns-width ttns) 838 with h = (ttns-height ttns) 839 for i from (1- w) downto 0 840 do 841 (loop 842 for j below (1- h) by 2 843 for u = (site-by-id ttns (list i j)) 844 for v = (site-by-bond ttns u :y+) 845 do 846 (multiple-value-bind (l r) 847 (apply-double-site-operator-y 848 double 849 (tensor-by-site ttns u) 850 (tensor-by-site ttns v) 851 maxdim 852 :left) 853 (setf (tensor-by-site ttns u) l) 854 (setf (tensor-by-site ttns v) r)) 855 (unless (>= j (- h 3)) 856 (normalize-tensor-by-site ttns v :y+))) 857 (loop 858 for j from (1- h) above 0 859 do 860 (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-)) 861 (unless (zerop i) 862 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x-)) 863 )) 864 865 (defun apply-odd-tebd-operators-y-down (double ttns maxdim) 866 (loop 867 with w = (ttns-width ttns) 868 with h = (ttns-height ttns) 869 for i below w 870 do 871 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :y+) 872 (loop 873 for j from 1 below (1- h) by 2 874 for u = (site-by-id ttns (list i j)) 875 for v = (site-by-bond ttns u :y+) 876 do 877 (multiple-value-bind (l r) 878 (apply-double-site-operator-y 879 double 880 (tensor-by-site ttns u) 881 (tensor-by-site ttns v) 882 maxdim 883 :left) 884 (setf (tensor-by-site ttns u) l) 885 (setf (tensor-by-site ttns v) r)) 886 (unless (>= j (- h 3)) 887 (normalize-tensor-by-site ttns v :y+))) 888 (loop 889 for j from (1- h) above 0 890 do 891 (normalize-tensor-by-site ttns (site-by-id ttns (list i j)) :y-)) 892 (unless (eql i (1- w)) 893 (normalize-tensor-by-site ttns (site-by-id ttns (list i 0)) :x+)))) 894 895 (defun make-ssps (width height dimension interaction potential step) 896 (loop 897 with a = (make-array (list width height)) 898 for i below width 899 do 900 (loop 901 for j below height 902 for h = (bose-hubbard-site-hamiltonian dimension 1 interaction (aref potential i j)) 903 do 904 (setf (aref a i j) (make-single-site-propagator h step))) 905 finally (return a))) 906 907 (defun make-dsps (width height dimension interaction potential step) 908 (loop 909 with a = (make-array (list (1- width) height)) 910 for i below (1- width) 911 do 912 (loop 913 for j below height 914 for f = (bose-hubbard-site-hamiltonian dimension 1 interaction (aref potential i j)) 915 for g = (bose-hubbard-site-hamiltonian dimension 1 interaction (aref potential (1+ i) j)) 916 do 917 (setf (aref a i j) (make-double-site-propagator f g step))) 918 finally (return a))) 919 920 (defun ttns-tebd-evolve (ttns ssps dsps idsp maxdim) 921 (apply-even-tebd-operators-y idsp ttns maxdim) 922 (apply-odd-tebd-operators-y idsp ttns maxdim) 923 (apply-even-tebd-operators-x ssps dsps ttns maxdim) 924 (apply-odd-tebd-operators-x ssps dsps ttns maxdim) 925 (apply-odd-tebd-operators-x-down ssps dsps ttns maxdim) 926 (apply-even-tebd-operators-x-down ssps dsps ttns maxdim) 927 (apply-odd-tebd-operators-y-down idsp ttns maxdim) 928 (apply-even-tebd-operators-y-down idsp ttns maxdim)) 929 930 (defun all-pn-operators (width height dimension) 931 (loop 932 with a = (make-array (list width height)) 933 for i below width 934 do 935 (loop 936 for j below height 937 do 938 (setf (aref a i j) 939 (pn-operator width height i j dimension))) 940 finally 941 (return a))) 942 943 (defun all-pn (ttns norm pn) 944 (loop 945 with w = (array-dimension pn 0) 946 with h = (array-dimension pn 1) 947 with a = (make-array (list w h)) 948 for i below w 949 do 950 (loop 951 for j below h 952 do 953 (setf (aref a i j) 954 (realpart (/ (ttns-expectation 955 ttns 956 ttns 957 (aref pn i j)) 958 norm)))) 959 finally 960 (return a))) 961 962 (defun print-pn (pn) 963 (loop 964 with w = (array-dimension pn 0) 965 with h = (array-dimension pn 1) 966 for j from (1- h) downto 0 967 do 968 (loop 969 for i below w 970 do 971 (format t 972 ;; "~,3f~:[ ~;~%~]" 973 "~F~:[ ~;~%~]" 974 (aref pn i j) 975 (eql i (1- w))))))