trotzky.lisp (4293B)
1 ;; equilibration in an optical lattice 2 3 ;; parameters 4 5 (defparameter *state* #p"inputs/trotzky-l-21-n-7/") 6 7 (defparameter *j* 4.1d-1) 8 (defparameter *k* 4.1d-3) 9 10 (defparameter *maxdim* '(nil nil 100)) 11 12 ;; libraries 13 14 (load "conf.lisp") 15 (require :tebdol) 16 17 (use-package '(:tensor :mpi :part :mps :tebd :bhm)) 18 (setf *print-level* 2) 19 20 ;; hamiltonian 21 22 (defun make-hamiltonian (dimension length j k) 23 (loop 24 with c = (/ (1- length) 2) 25 with m = (make-array length) 26 for i from (1+ (- c)) below 0 27 do 28 (let ((o (bose-hubbard-site-hamiltonian dimension j 1 (* 1/2 k i i)))) 29 (setf (svref m (+ c i)) o 30 (svref m (- c i)) o)) 31 finally 32 (setf (svref m c) (bose-hubbard-site-hamiltonian dimension j 1 (* 1/2 k 0 0))) 33 (let ((o (bose-hubbard-site-hamiltonian dimension j 1 (* 1/2 k c c)))) 34 (setf (svref m 0) (tensor-contract (left-open-boundary-tensor o) 1 o 0)) 35 (setf (svref m (1- length)) (tensor-contract o 3 (right-open-boundary-tensor o) 0))) 36 (return m))) 37 38 ;; odd particle number operator 39 ;; 40 ;; 0 0 41 ;; +------+ 42 ;; 0 | 1 0 | 43 ;; 0 | n 1 | 44 ;; +------+ 45 ;; 46 ;; 0 0 47 ;; +------+ 48 ;; 0 | 1 0 | 49 ;; 0 | 0 1 | 50 ;; +------+ 51 52 (defun site-operator (dimension type) 53 (let ((ops (list 54 (cons '1 #'(lambda (n) (declare (ignore n)) 1)) 55 (cons 'n #'(lambda (n) n)))) 56 (tbl (flet ((f (x y z) 57 (cons (list (make-subscript 58 :numbers (list (first x)) 59 :subscript (second x)) 60 (make-subscript 61 :numbers (list (first y)) 62 :subscript (second y))) 63 z))) 64 (case type 65 (:identity 66 (list 67 (f '(0 0) '(0 0) '1) 68 (f '(0 1) '(0 1) '1))) 69 (:particle-number 70 (list 71 (f '(0 0) '(0 0) '1) 72 (f '(0 1) '(0 0) 'n) 73 (f '(0 1) '(0 1) '1))))))) 74 (let ((l (list (make-segment :numbers '(0) :dimension 2)))) 75 (functional-tensor 76 (list l (ket-physical-index dimension) (bra-physical-index dimension) l) 77 #'(lambda (s) 78 (let ((x (find 79 (list (first s) (fourth s)) 80 tbl 81 :test #'equalp 82 :key #'car))) 83 (if x 84 (funcall (cdr (find (cdr x) ops :key #'car)) 85 (- (first (subscript-numbers (third s))))) 86 0))))))) 87 88 (defun make-particle-number-mpo (dimension length function) 89 (loop 90 with c = (/ (1- length) 2) 91 with m = (make-array length) 92 with e = (site-operator dimension :identity) 93 with p = (site-operator dimension :particle-number) 94 for i from (1+ (- c)) to (1- c) 95 do 96 (setf (svref m (+ i c)) 97 (if (funcall function i) p e)) 98 finally 99 (let ((o (if (funcall function (- c)) p e))) 100 (setf (svref m 0) (tensor-contract (left-open-boundary-tensor o) 1 o 0))) 101 (let ((o (if (funcall function c) p e))) 102 (setf (svref m (1- length)) (tensor-contract o 3 (right-open-boundary-tensor o) 0))) 103 (return m))) 104 105 ;; initialization 106 107 (mpi-init) 108 109 (defparameter *rank* (mpi-comm-rank *mpi-comm-world*)) 110 (defparameter *size* (mpi-comm-size *mpi-comm-world*)) 111 112 (defparameter *length* (state-parameter *state* :model-length)) 113 (defparameter *particles* (state-parameter *state* :particle-number)) 114 (defparameter *dimension* (state-parameter *state* :site-dimension)) 115 116 (defparameter *partition* (uniform-partition *size* *length*)) 117 118 (defparameter *mps* (load-mps *rank* *partition* *state*)) 119 (defparameter *mph* (make-hamiltonian *dimension* *length* *j* *k*)) 120 (defparameter *mpn* (make-particle-number-mpo *dimension* *length* #'oddp)) 121 122 (defparameter *count* 100) 123 (defparameter *step* (/ (* 5 (/ pi 2 *j*)) *count*)) 124 (defparameter *propagator* (make-mph-propagator *mph* *step*)) 125 126 (when (zerop *rank*) 127 (format 128 t 129 "# processes = ~a; length = ~a; particles = ~a; maxdim = ~a; j = ~f; k = ~f;~%" 130 *size* *length* *particles* *maxdim* *j* *k*) 131 (format t "# time nodd norm~%")) 132 133 (mps-prepare *rank* *length* *mps* *maxdim*) 134 135 ;; main loop 136 137 (loop for i from 0 to *count* do 138 (let ((e (realpart (mps-expectation-value *size* *rank* *length* *mpn* *mps*))) 139 (n (realpart (mps-overlap *size* *rank* *length* *mps* *mps*)))) 140 (if (= *rank* (1- *size*)) 141 (format t "~,2f ~f ~f~%" 142 (/ (* i *step*) (/ pi 2 *j*)) 143 (/ e n *particles*) 144 (sqrt n)))) 145 (when (< i *count*) 146 (tebd-evolve *rank* *length* *propagator* *mps* *maxdim*))) 147 148 (mpi-finalize)