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