tebdol

Simulation of ultracold atoms in optical lattices
git clone https://miroslavurbanek.com/git/tebdol.git
Log | Files | Refs | README

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)