ZaletelOps

Fortran module ZaletelOps: August 2017 (dj, updated)

Constructs the MPO for time evolution with the Zaletel algorithm based on the MPO of the Hamiltonian.

Authors

• 13. 12. Wall

• 4. Jaschke

Details

The following procedures are important for call from other modules.

Procedures	include.f90	mpi.f90
zaletel_2nd
zaletel_4th

ZaletelOps_f90.extract_op_from_zaletel()

fortran-subroutine - April 2016 (dj) Extract from the big exponential the smaller operator.

Arguments

OpTYPE(tensorc), inout

On exit, this is the matrix in the MPO to be applied to the local Hilbert space. Needs to be allocated before calling this subroutine.

ZalTYPE(tensorc), in

Is the operator in the Zaletel notation on the complete operator exponential space.

b1INTEGER, in

Pick element corresponding to b1.

b2INTEGER, in

Pick element corresponding to b2.

locdINTEGER, in

Local dimension of the Hilbert space.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine extract_op_from_zaletel(Op, Zal, b1, b2, locd)
      type(tensorc), intent(inout) :: Op
      type(tensorc), intent(in) :: Zal
      integer, intent(in) :: b1, b2, locd

      ! Local variables
      ! ---------------

      ! indices
      integer :: i1, i2, j1, j2

      ! temporary array
      complex(KIND=rKind), dimension(:, :), allocatable :: tmp

      allocate(tmp(Zal%dl(1), Zal%dl(2)))
      tmp = reshape(Zal%elem(:Zal%dim), [Zal%dl(1), Zal%dl(2)])

      i1 = b1 * 2 * locd + b2 * locd + 1
      i2 = b1 * 2 * locd + b2 * locd + locd
      j1 = 1
      j2 = locd
      Op%elem = reshape(tmp(i1:i2, j1:j2), [product(Op%dl)])

      deallocate(tmp)

  end subroutine extract_op_from_zaletel

ZaletelOps_f90.build_hcb_idb()

fortran-subroutine - April 2016 (dj) Build the matrices Hcb and Idb, where Hcb is the creation operator for hard-core bosons and Idb is the identity matrix.

Arguments

HcbTYPE(tensor), out

Creation operator .

IdbTYPE(tensor), out

Creation operator .

Details

(template defined in ZaletelOps_template.f90)

Source Code

show / hide f90 code

  subroutine build_hcb_idb(Hcb, Idb)
      type(tensor), intent(out) :: Hcb, Idb

      call create(Hcb, [2, 2], init='0')
      Hcb%elem(2) = done

      call create(Idb, [2, 2], init='0')
      Idb%elem(1) = done
      Idb%elem(4) = done

  end subroutine build_hcb_idb

ZaletelOps_f90.get_zaletel_expm_mpo()

fortran-subroutine - April 2016 (dj) Build the exponential describing = langle 0 bar{0} | c_{a} c_{bar{a}}` .

Arguments

MatTYPE(matrixc), out

Contains on exit the argument in the exponential described above.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

AATYPE(tensor), in

Matrix belonging to the block A within the Zaletel nomenclature.

BBTYPE(tensor), in

Matrix belonging to the block B within the Zaletel nomenclature.

CCTYPE(tensor), in

Matrix belonging to the block C within the Zaletel nomenclature.

DDTYPE(tensor), in

Matrix belonging to the block D within the Zaletel nomenclature.

HcbTYPE(tensor), in

Auxiliary matrix (representing c^{dagger}, always?).

IdbTYPE(tensor), in

Auxiliary matrix (representing Id, always?).

locdINTEGER, in

Dimension of the local Hilbert space.

Source Code

show / hide f90 code

  subroutine get_zaletel_expm_mpo(Mat, dt, AA, BB, CC, DD, Hcb, Idb, locd)
      type(tensorc), intent(inout) :: Mat
      complex(kind=rkind), intent(in) :: dt
      type(tensor), intent(in) :: AA, BB, CC, DD
      type(tensor), intent(in) :: Hcb, Idb
      integer, intent(in) :: locd

      ! Local variables
      ! ---------------

      ! looping over the left operator c^{\dagger} respectively Id
      integer :: l1, l2, ll

      ! looping over the right operator \bar{c}^{\dagger} respectively Id
      integer :: r1, r2, rr

      ! subspace to be set in the matrix Arg%elem (i <-> row, j <-> col)
      integer :: i1, i2, j1, j2

      ! matrix in the exponential
      complex(KIND=rKind), dimension(:, :), allocatable :: marg

      allocate(marg(Mat%dl(1), Mat%dl(2)))
      marg = 0.0_rKind

      ! Take care of A with c^{\dagger} A \bar{c}^{\dagger}
      ! ---------------------------------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Hcb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Hcb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Hcb%elem(ll) &
                      * reshape(AA%elem, [AA%dl(1), AA%dl(2)]) &
                      * Hcb%elem(rr))
          end do
          end do
      end do
      end do

      ! Take care of B with c^{\dagger} B Id
      ! ------------------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Hcb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Idb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Hcb%elem(ll) &
                      * reshape(BB%elem, [BB%dl(1), BB%dl(2)]) &
                      * Idb%elem(rr) * sqrt(dt))
          end do
          end do
      end do
      end do

      ! Take care of C with Id C \bar{c}^{\dagger}
      ! ------------------------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Idb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Hcb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Idb%elem(ll) &
                      * reshape(CC%elem, [CC%dl(1), CC%dl(2)]) &
                      * Hcb%elem(rr) * sqrt(dt))
          end do
          end do
      end do
      end do

      ! Take care of D with Id D Id
      ! ---------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Idb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Idb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Idb%elem(ll) &
                      * reshape(DD%elem, [DD%dl(1), DD%dl(2)]) &
                      * Idb%elem(rr) * dt)
          end do
          end do
      end do
      end do

      call exp_kit_complex(Mat%elem, done, marg, Mat%dl(1))
      deallocate(marg)

  end subroutine get_zaletel_expm_mpo

ZaletelOps_f90.get_zaletel_expm_mpoc()

fortran-subroutine - April 2016 (dj) Build the exponential describing = langle 0 bar{0} | c_{a} c_{bar{a}}` .

Arguments

MatTYPE(matrixc), out

Contains on exit the argument in the exponential described above.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

AATYPE(tensorc), in

Matrix belonging to the block A within the Zaletel nomenclature.

BBTYPE(tensorc), in

Matrix belonging to the block B within the Zaletel nomenclature.

CCTYPE(tensorc), in

Matrix belonging to the block C within the Zaletel nomenclature.

DDTYPE(tensorc), in

Matrix belonging to the block D within the Zaletel nomenclature.

HcbTYPE(tensor), in

Auxiliary matrix (representing c^{dagger}, always?).

IdbTYPE(tensor), in

Auxiliary matrix (representing Id, always?).

locdINTEGER, in

Dimension of the local Hilbert space.

Source Code

show / hide f90 code

  subroutine get_zaletel_expm_mpoc(Mat, dt, AA, BB, CC, DD, Hcb, Idb, locd)
      type(tensorc), intent(inout) :: Mat
      complex(kind=rkind), intent(in) :: dt
      type(tensorc), intent(in) :: AA, BB, CC, DD
      type(tensor), intent(in) :: Hcb, Idb
      integer, intent(in) :: locd

      ! Local variables
      ! ---------------

      ! looping over the left operator c^{\dagger} respectively Id
      integer :: l1, l2, ll

      ! looping over the right operator \bar{c}^{\dagger} respectively Id
      integer :: r1, r2, rr

      ! subspace to be set in the matrix Arg%elem (i <-> row, j <-> col)
      integer :: i1, i2, j1, j2

      ! matrix in the exponential
      complex(KIND=rKind), dimension(:, :), allocatable :: marg

      allocate(marg(Mat%dl(1), Mat%dl(2)))
      marg = 0.0_rKind

      ! Take care of A with c^{\dagger} A \bar{c}^{\dagger}
      ! ---------------------------------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Hcb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Hcb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Hcb%elem(ll) &
                      * reshape(AA%elem, [AA%dl(1), AA%dl(2)]) &
                      * Hcb%elem(rr))
          end do
          end do
      end do
      end do

      ! Take care of B with c^{\dagger} B Id
      ! ------------------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Hcb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Idb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Hcb%elem(ll) &
                      * reshape(BB%elem, [BB%dl(1), BB%dl(2)]) &
                      * Idb%elem(rr) * sqrt(dt))
          end do
          end do
      end do
      end do

      ! Take care of C with Id C \bar{c}^{\dagger}
      ! ------------------------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Idb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Hcb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Idb%elem(ll) &
                      * reshape(CC%elem, [CC%dl(1), CC%dl(2)]) &
                      * Hcb%elem(rr) * sqrt(dt))
          end do
          end do
      end do
      end do

      ! Take care of D with Id D Id
      ! ---------------------------

      ll = 0
      do l2 = 1, 2
      do l1 = 1, 2
          ll = ll + 1
          if(Idb%elem(ll) == 0.0_rKind) cycle

          rr = 0
          do r2 = 1, 2
          do r1 = 1, 2
              rr = rr + 1
              if(Idb%elem(rr) == 0.0_rKind) cycle

              ! Get the block of the matrix
              i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + 1
              i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + locd
              j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + 1
              j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + locd

              marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                   + (Idb%elem(ll) &
                      * reshape(DD%elem, [DD%dl(1), DD%dl(2)]) &
                      * Idb%elem(rr) * dt)
          end do
          end do
      end do
      end do

      call exp_kit_complex(Mat%elem, done, marg, Mat%dl(1))
      deallocate(marg)

  end subroutine get_zaletel_expm_mpoc

ZaletelOps_f90.get_zaletel_expm_qmpo()

fortran-subroutine - April 2016 (dj) Build the exponential describing = langle 0 bar{0} | c_{a} c_{bar{a}}` .

Arguments

MatTYPE(matrixc), out

Contains on exit the argument in the exponential described above.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

AATYPE(qtensor), in

Matrix belonging to the block A within the Zaletel nomenclature.

BBTYPE(qtensor), in

Matrix belonging to the block B within the Zaletel nomenclature.

CCTYPE(qtensor), in

Matrix belonging to the block C within the Zaletel nomenclature.

DDTYPE(qtensor), in

Matrix belonging to the block D within the Zaletel nomenclature.

HcbTYPE(tensor), in

Auxiliary matrix (representing c^{dagger}, always?).

IdbTYPE(tensor), in

Auxiliary matrix (representing Id, always?).

locdINTEGER, in

Dimension of the local Hilbert space.

DmapTYPE(vector_int), in

Saving the dimensions of a block according to the hash.

StrmapTYPE(vector_int), in

Saving the sum of the previous dimension of smaller hashes in order to find the starting index writing a non-blocked matrix.

QidTYPE(qtensor), in

Identity matrix in q-notation.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_expm_qmpo(Mat, dt, AA, BB, CC, DD, Hcb, Idb, &
                                       locd, Dmap, Strmap, Qid)
      type(tensorc), intent(inout) :: Mat
      complex(kind=rkind), intent(in) :: dt
      type(qtensor), intent(in) :: AA, BB, CC, DD, Qid
      type(tensor), intent(in) :: Hcb, Idb
      integer, intent(in) :: locd
      type(vector_int) :: dmap, strmap

      ! Local variables
      ! ---------------

      ! looping over the left operator c^{\dagger} respectively Id
      integer :: l1, l2, ll

      ! looping over the right operator \bar{c}^{\dagger} respectively Id
      integer :: r1, r2, rr

      ! subspace to be set in the matrix Arg%elem (i <-> row, j <-> col)
      integer :: i1, i2, j1, j2

      ! indices with matching hash
      integer :: q1, q2

      ! looping over blocks in q-matrix
      integer :: kk

      ! save temporary hashes
      real(KIND=rKind) :: hsh1, hsh2

      ! matrix in the exponential
      complex(KIND=rKind), dimension(:, :), allocatable :: marg

      allocate(marg(Mat%dl(1), Mat%dl(1)))
      marg = 0.0_rKind

      ! Take care of A with c^{\dagger} A \bar{c}^{\dagger}
      ! ---------------------------------------------------

      do kk = 1, AA%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(AA, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(AA, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Hcb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Hcb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Hcb%elem(ll) &
                          * reshape(AA%Data(kk)%Tens%elem, &
                                    [AA%Data(kk)%Tens%dl(1), &
                                     AA%Data(kk)%Tens%dl(2)]) &
                          * Hcb%elem(rr))
              end do
              end do
          end do
          end do
      end do

      ! Take care of B with c^{\dagger} B Id
      ! ------------------------------------

      do kk = 1, BB%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(BB, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(BB, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Hcb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Idb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Hcb%elem(ll) &
                          * reshape(BB%Data(kk)%Tens%elem, &
                                    [BB%Data(kk)%Tens%dl(1), &
                                     BB%Data(kk)%Tens%dl(2)]) &
                          * Idb%elem(rr) * sqrt(dt))
              end do
              end do
          end do
          end do
      end do

      ! Take care of C with Id C \bar{c}^{\dagger}
      ! ------------------------------------------

      do kk = 1, CC%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(CC, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(CC, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Idb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Hcb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Idb%elem(ll) &
                          * reshape(CC%Data(kk)%Tens%elem, &
                                    [CC%Data(kk)%Tens%dl(1), &
                                     CC%Data(kk)%Tens%dl(2)]) &
                          * Hcb%elem(rr) * sqrt(dt))
              end do
              end do
          end do
          end do
      end do

      ! Take care of D with Id D Id
      ! ---------------------------

      do kk = 1, DD%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(DD, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(DD, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Idb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Idb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Idb%elem(ll) &
                          * reshape(DD%Data(kk)%Tens%elem, &
                                    [DD%Data(kk)%Tens%dl(1), &
                                     DD%Data(kk)%Tens%dl(2)]) &
                          * Idb%elem(rr) * (dt))
              end do
              end do
          end do
          end do
      end do

      call exp_kit_complex(Mat%elem, done, marg, Mat%dl(1))
      deallocate(marg)

  end subroutine get_zaletel_expm_qmpo

ZaletelOps_f90.get_zaletel_expm_qmpoc()

fortran-subroutine - April 2016 (dj) Build the exponential describing = langle 0 bar{0} | c_{a} c_{bar{a}}` .

Arguments

MatTYPE(matrixc), out

Contains on exit the argument in the exponential described above.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

AATYPE(qtensorc), in

Matrix belonging to the block A within the Zaletel nomenclature.

BBTYPE(qtensorc), in

Matrix belonging to the block B within the Zaletel nomenclature.

CCTYPE(qtensorc), in

Matrix belonging to the block C within the Zaletel nomenclature.

DDTYPE(qtensorc), in

Matrix belonging to the block D within the Zaletel nomenclature.

HcbTYPE(tensor), in

Auxiliary matrix (representing c^{dagger}, always?).

IdbTYPE(tensor), in

Auxiliary matrix (representing Id, always?).

locdINTEGER, in

Dimension of the local Hilbert space.

DmapTYPE(vector_int), in

Saving the dimensions of a block according to the hash.

StrmapTYPE(vector_int), in

Saving the sum of the previous dimension of smaller hashes in order to find the starting index writing a non-blocked matrix.

QidTYPE(qtensorc), in

Identity matrix in q-notation.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_expm_qmpoc(Mat, dt, AA, BB, CC, DD, Hcb, Idb, &
                                       locd, Dmap, Strmap, Qid)
      type(tensorc), intent(inout) :: Mat
      complex(kind=rkind), intent(in) :: dt
      type(qtensorc), intent(in) :: AA, BB, CC, DD, Qid
      type(tensor), intent(in) :: Hcb, Idb
      integer, intent(in) :: locd
      type(vector_int) :: dmap, strmap

      ! Local variables
      ! ---------------

      ! looping over the left operator c^{\dagger} respectively Id
      integer :: l1, l2, ll

      ! looping over the right operator \bar{c}^{\dagger} respectively Id
      integer :: r1, r2, rr

      ! subspace to be set in the matrix Arg%elem (i <-> row, j <-> col)
      integer :: i1, i2, j1, j2

      ! indices with matching hash
      integer :: q1, q2

      ! looping over blocks in q-matrix
      integer :: kk

      ! save temporary hashes
      real(KIND=rKind) :: hsh1, hsh2

      ! matrix in the exponential
      complex(KIND=rKind), dimension(:, :), allocatable :: marg

      allocate(marg(Mat%dl(1), Mat%dl(1)))
      marg = 0.0_rKind

      ! Take care of A with c^{\dagger} A \bar{c}^{\dagger}
      ! ---------------------------------------------------

      do kk = 1, AA%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(AA, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(AA, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Hcb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Hcb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Hcb%elem(ll) &
                          * reshape(AA%Data(kk)%Tens%elem, &
                                    [AA%Data(kk)%Tens%dl(1), &
                                     AA%Data(kk)%Tens%dl(2)]) &
                          * Hcb%elem(rr))
              end do
              end do
          end do
          end do
      end do

      ! Take care of B with c^{\dagger} B Id
      ! ------------------------------------

      do kk = 1, BB%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(BB, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(BB, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Hcb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Idb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Hcb%elem(ll) &
                          * reshape(BB%Data(kk)%Tens%elem, &
                                    [BB%Data(kk)%Tens%dl(1), &
                                     BB%Data(kk)%Tens%dl(2)]) &
                          * Idb%elem(rr) * sqrt(dt))
              end do
              end do
          end do
          end do
      end do

      ! Take care of C with Id C \bar{c}^{\dagger}
      ! ------------------------------------------

      do kk = 1, CC%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(CC, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(CC, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Idb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Hcb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Idb%elem(ll) &
                          * reshape(CC%Data(kk)%Tens%elem, &
                                    [CC%Data(kk)%Tens%dl(1), &
                                     CC%Data(kk)%Tens%dl(2)]) &
                          * Hcb%elem(rr) * sqrt(dt))
              end do
              end do
          end do
          end do
      end do

      ! Take care of D with Id D Id
      ! ---------------------------

      do kk = 1, DD%nb
          ! Find corresponding block in identity matrix
          hsh1 = get_hash(DD, [1], kk)
          q1 = FindTagIndex(hsh1, Qid%hash(:Qid%nb))

          hsh2 = get_hash(DD, [2], kk)
          q2 = FindTagIndex(hsh2, Qid%hash(:Qid%nb))

          ll = 0
          do l2 = 1, 2
          do l1 = 1, 2
              ll = ll + 1
              if(Idb%elem(ll) == 0.0_rKind) cycle

              rr = 0
              do r2 = 1, 2
              do r1 = 1, 2
                  rr = rr + 1
                  if(Idb%elem(rr) == 0.0_rKind) cycle

                  ! Get the block of the matrix
                  i1 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       1 + Strmap%elem(q1)
                  i2 = (l1 - 1) * 2 * locd + (r1 - 1) * locd + &
                       Dmap%elem(q1) + Strmap%elem(q1)
                  j1 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       1 + Strmap%elem(q2)
                  j2 = (l2 - 1) * 2 * locd + (r2 - 1) * locd + &
                       Dmap%elem(q2) + Strmap%elem(q2)

                  marg(i1:i2, j1:j2) = marg(i1:i2, j1:j2) &
                       + (Idb%elem(ll) &
                          * reshape(DD%Data(kk)%Tens%elem, &
                                    [DD%Data(kk)%Tens%dl(1), &
                                     DD%Data(kk)%Tens%dl(2)]) &
                          * Idb%elem(rr) * (dt))
              end do
              end do
          end do
          end do
      end do

      call exp_kit_complex(Mat%elem, done, marg, Mat%dl(1))
      deallocate(marg)

  end subroutine get_zaletel_expm_qmpoc

ZaletelOps_f90.get_zaletel_1st_site_tensor()

fortran-subroutine - April 2016 (dj) Build the entries of the Zaletel MPO for the first (left) site in the system.

Arguments

ZalTYPE(sr_matrix_tensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of first (left) site.

HamTYPE(sr_matrix_tensor), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the first site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_1st_site_tensor(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_tensorc), intent(out) :: Zal
      type(sr_matrix_tensor), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ===============

      ! for looping / indexing
      integer :: k1, k2

      ! left and right chi of the Zaletel MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      ! Different parts in MPO
      type(tensorc) :: DD
      type(tensorc), pointer :: AA(:), BB(:), CC(:)

      ! Temporay matrix
      type(tensorc) :: Tmp

      ! zaletel exponential
      type(tensorc) :: ZalU

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_1st_site_tensor'

      call build_hcb_idb(Hcb, Idb)

      ! Get bond dimensions
      chil = 1
      chir = Ham%cbd - 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      k1 = 1
      allocate(Zal%row(chil))

      nel = Ham%row(Ham%rbd)%numel - 1
      Zal%Row(k1)%numel = nel
      allocate(Zal%row(k1)%ind(nel), Zal%row(k1)%Op(nel))
      Zal%row(k1)%ind = Ham%row(Ham%rbd)%ind(1:nel)

      locd = Ham%Row(1)%Op(1)%dl(1)

      ! Copy/set terms from Ham
      ! -----------------------

      k2 = 1
      call copy(DD, Ham%row(k1)%Op(k2))

      allocate(AA(chir - 1), BB(chir - 1), CC(chir - 1))
      do k2 = 1, nel - 1
          ! Create empty AA, BB
          call create(AA(k2), DD%dl, init='0')
          call create(BB(k2), DD%dl, init='0')

          ! Copy CC
          call copy(CC(k2), Ham%Row(k1)%Op(k2 + 1))
      end do

      ! W_D term is just exp(D)
      ! -----------------------

      k2 = 1
      call create(Zal%row(k1)%Op(k2), Ham%row(k1)%Op(k2)%dl)
      call copy(Tmp, DD)
      if(mpo_is_hermitian) then
          call expmh(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      else
          call expm(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      end if
      call destroy(Tmp)

      call create(ZalU, [4 * locd, 4 * locd], init='0')

      do k2 = 1, nel - 1
          call get_zaletel_expm(ZalU, dt, AA(k2), BB(k2), CC(k2), DD, &
                               Hcb, Idb, locd)
          call create(Zal%Row(k1)%Op(k2 + 1), DD%dl)
          call extract_op_from_zaletel(Zal%Row(k1)%Op(k2 + 1), &
                                       ZalU, 0, 1, locd)

          call destroy(AA(k2))
          call destroy(BB(k2))
          call destroy(CC(k2))
      end do

      call destroy(DD)
      call destroy(ZalU)
      deallocate(AA, BB, CC)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_1st_site_tensor <<< <<<'

  end subroutine get_zaletel_1st_site_tensor

ZaletelOps_f90.get_zaletel_1st_site_tensorc()

fortran-subroutine - April 2016 (dj) Build the entries of the Zaletel MPO for the first (left) site in the system.

Arguments

ZalTYPE(sr_matrix_tensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of first (left) site.

HamTYPE(sr_matrix_tensorc), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the first site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_1st_site_tensorc(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_tensorc), intent(out) :: Zal
      type(sr_matrix_tensorc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ===============

      ! for looping / indexing
      integer :: k1, k2

      ! left and right chi of the Zaletel MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      ! Different parts in MPO
      type(tensorc) :: DD
      type(tensorc), pointer :: AA(:), BB(:), CC(:)

      ! Temporay matrix
      type(tensorc) :: Tmp

      ! zaletel exponential
      type(tensorc) :: ZalU

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_1st_site_tensorc'

      call build_hcb_idb(Hcb, Idb)

      ! Get bond dimensions
      chil = 1
      chir = Ham%cbd - 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      k1 = 1
      allocate(Zal%row(chil))

      nel = Ham%row(Ham%rbd)%numel - 1
      Zal%Row(k1)%numel = nel
      allocate(Zal%row(k1)%ind(nel), Zal%row(k1)%Op(nel))
      Zal%row(k1)%ind = Ham%row(Ham%rbd)%ind(1:nel)

      locd = Ham%Row(1)%Op(1)%dl(1)

      ! Copy/set terms from Ham
      ! -----------------------

      k2 = 1
      call copy(DD, Ham%row(k1)%Op(k2))

      allocate(AA(chir - 1), BB(chir - 1), CC(chir - 1))
      do k2 = 1, nel - 1
          ! Create empty AA, BB
          call create(AA(k2), DD%dl, init='0')
          call create(BB(k2), DD%dl, init='0')

          ! Copy CC
          call copy(CC(k2), Ham%Row(k1)%Op(k2 + 1))
      end do

      ! W_D term is just exp(D)
      ! -----------------------

      k2 = 1
      call create(Zal%row(k1)%Op(k2), Ham%row(k1)%Op(k2)%dl)
      call copy(Tmp, DD)
      if(mpo_is_hermitian) then
          call expmh(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      else
          call expm(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      end if
      call destroy(Tmp)

      call create(ZalU, [4 * locd, 4 * locd], init='0')

      do k2 = 1, nel - 1
          call get_zaletel_expm(ZalU, dt, AA(k2), BB(k2), CC(k2), DD, &
                               Hcb, Idb, locd)
          call create(Zal%Row(k1)%Op(k2 + 1), DD%dl)
          call extract_op_from_zaletel(Zal%Row(k1)%Op(k2 + 1), &
                                       ZalU, 0, 1, locd)

          call destroy(AA(k2))
          call destroy(BB(k2))
          call destroy(CC(k2))
      end do

      call destroy(DD)
      call destroy(ZalU)
      deallocate(AA, BB, CC)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_1st_site_tensorc <<< <<<'

  end subroutine get_zaletel_1st_site_tensorc

ZaletelOps_f90.get_zaletel_1st_site_qtensor()

fortran-subroutine - April 2016 (dj) Build the entries of the Zaletel MPO for the first (left) site in the system.

Arguments

ZalTYPE(sr_matrix_qtensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of first (left) site.

HamTYPE(sr_matrix_qtensor), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the first site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_1st_site_qtensor(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_qtensorc), intent(out) :: Zal
      type(sr_matrix_qtensor), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! bond dimensions of MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! for indexing / looping
      integer :: k1, k2, jj, idx

      ! number of elements in sparse row
      integer :: nel

      ! for mapping between matrices and q-matrices
      type(vector_int) :: Dmap, Strmap

      ! Different parts in MPO
      type(qtensorc) :: DD
      type(qtensorc), pointer :: AA(:), BB(:), CC(:)

      ! Zaletel exponential and propagator for local Hilbert space
      type(tensorc) :: ZalU, Unonq

      ! ID as q-matrix / temporary q-matrix
      type(qtensorc) :: Qid, Tmp

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_first_site_qtensor'

      call build_hcb_idb(Hcb, Idb)

      chil = 1
      chir = Ham%cbd - 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      k1 = 1
      allocate(Zal%row(chil))

      nel = Ham%Row(Ham%rbd)%numel - 1
      Zal%Row(k1)%numel = nel
      allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
      Zal%Row(k1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

      ! Copy/set terms from Ham
      ! -----------------------

      call copy(Qid, Ham%Row(Ham%rbd)%Op(Ham%cbd))

      k2 = 1
      call copy(DD, Ham%Row(k1)%Op(k2))

      allocate(AA(chir - 1), BB(chir - 1), CC(chir - 1))
      do k2 = 1, nel - 1
          ! Create empty AA, BB
          call create(AA(k2), Qid%nqs, 1)
          call create(BB(k2), Qid%nqs, 1)

          ! Copy C
          call copy(CC(k2), Ham%Row(k1)%Op(k2 + 1))
      end do

      ! W_D term is just exp(D) and D is blockdiagonal
      ! ----------------------------------------------

      ! Get the blocks, set them to zero in case local terms does not
      ! access all

      k2 = 1
      call copy(Zal%Row(k1)%Op(k2), Qid)
      call copy(Tmp, DD)

      ! Exponentiate
      do jj = 1, Tmp%nb
          idx = FindTagIndex(Tmp%hash(jj), &
                             Zal%Row(k1)%Op(k2)%hash(1:Zal%Row(k1)%Op(k2)%nb))
          if(mpo_is_hermitian) then
              call expmh(Zal%Row(k1)%Op(k2)%data(idx)%Tens, &
                        dt, Tmp%Data(jj)%Tens, 1)
          else
              call expm(Zal%Row(k1)%Op(k2)%data(idx)%Tens, &
                       dt, Tmp%Data(jj)%Tens, 1)
          end if
      end do

      call destroy(Tmp)
      call skim(Zal%Row(k1)%Op(k2), 10.0_rKind * numzero)

      ! Other terms W_X (cannot be created as q-matrices)
      ! ---------------

      ! Store identity matrix and get mapping from/to q-matrices
      call get_mapping_qtensorc(Qid, Dmap, Strmap, locd)

      call create(ZalU, [4 * locd, 4 * locd], init='0')
      call create(Unonq, [locd, locd], init='0')

      do k2 = 1, (nel - 1)
          call get_zaletel_expm(ZalU, dt, AA(k2), BB(k2), CC(k2), &
                               DD, Hcb, Idb, locd, Dmap, Strmap, &
                               Qid)
          call extract_op_from_zaletel(Unonq, ZalU, 0, 1, locd)
          call extract_qop_from_op(Zal%Row(k1)%Op(k2 + 1), Unonq, Dmap, &
                                   Strmap, locd, CC(k2), Qid)

          call destroy(AA(k2))
          call destroy(BB(k2))
          call destroy(CC(k2))
      end do

      deallocate(AA, BB, CC, Dmap%elem, Strmap%elem)
      call destroy(DD)
      call destroy(ZalU)
      call destroy(Unonq)
      call destroy(Qid)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_first_site_qtensor <<< <<<'

  end subroutine get_zaletel_1st_site_qtensor

ZaletelOps_f90.get_zaletel_1st_site_qtensorc()

fortran-subroutine - April 2016 (dj) Build the entries of the Zaletel MPO for the first (left) site in the system.

Arguments

ZalTYPE(sr_matrix_qtensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of first (left) site.

HamTYPE(sr_matrix_qtensorc), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the first site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_1st_site_qtensorc(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_qtensorc), intent(out) :: Zal
      type(sr_matrix_qtensorc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! bond dimensions of MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! for indexing / looping
      integer :: k1, k2, jj, idx

      ! number of elements in sparse row
      integer :: nel

      ! for mapping between matrices and q-matrices
      type(vector_int) :: Dmap, Strmap

      ! Different parts in MPO
      type(qtensorc) :: DD
      type(qtensorc), pointer :: AA(:), BB(:), CC(:)

      ! Zaletel exponential and propagator for local Hilbert space
      type(tensorc) :: ZalU, Unonq

      ! ID as q-matrix / temporary q-matrix
      type(qtensorc) :: Qid, Tmp

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_first_site_qtensorc'

      call build_hcb_idb(Hcb, Idb)

      chil = 1
      chir = Ham%cbd - 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      k1 = 1
      allocate(Zal%row(chil))

      nel = Ham%Row(Ham%rbd)%numel - 1
      Zal%Row(k1)%numel = nel
      allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
      Zal%Row(k1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

      ! Copy/set terms from Ham
      ! -----------------------

      call copy(Qid, Ham%Row(Ham%rbd)%Op(Ham%cbd))

      k2 = 1
      call copy(DD, Ham%Row(k1)%Op(k2))

      allocate(AA(chir - 1), BB(chir - 1), CC(chir - 1))
      do k2 = 1, nel - 1
          ! Create empty AA, BB
          call create(AA(k2), Qid%nqs, 1)
          call create(BB(k2), Qid%nqs, 1)

          ! Copy C
          call copy(CC(k2), Ham%Row(k1)%Op(k2 + 1))
      end do

      ! W_D term is just exp(D) and D is blockdiagonal
      ! ----------------------------------------------

      ! Get the blocks, set them to zero in case local terms does not
      ! access all

      k2 = 1
      call copy(Zal%Row(k1)%Op(k2), Qid)
      call copy(Tmp, DD)

      ! Exponentiate
      do jj = 1, Tmp%nb
          idx = FindTagIndex(Tmp%hash(jj), &
                             Zal%Row(k1)%Op(k2)%hash(1:Zal%Row(k1)%Op(k2)%nb))
          if(mpo_is_hermitian) then
              call expmh(Zal%Row(k1)%Op(k2)%data(idx)%Tens, &
                        dt, Tmp%Data(jj)%Tens, 1)
          else
              call expm(Zal%Row(k1)%Op(k2)%data(idx)%Tens, &
                       dt, Tmp%Data(jj)%Tens, 1)
          end if
      end do

      call destroy(Tmp)
      call skim(Zal%Row(k1)%Op(k2), 10.0_rKind * numzero)

      ! Other terms W_X (cannot be created as q-matrices)
      ! ---------------

      ! Store identity matrix and get mapping from/to q-matrices
      call get_mapping_qtensorc(Qid, Dmap, Strmap, locd)

      call create(ZalU, [4 * locd, 4 * locd], init='0')
      call create(Unonq, [locd, locd], init='0')

      do k2 = 1, (nel - 1)
          call get_zaletel_expm(ZalU, dt, AA(k2), BB(k2), CC(k2), &
                               DD, Hcb, Idb, locd, Dmap, Strmap, &
                               Qid)
          call extract_op_from_zaletel(Unonq, ZalU, 0, 1, locd)
          call extract_qop_from_op(Zal%Row(k1)%Op(k2 + 1), Unonq, Dmap, &
                                   Strmap, locd, CC(k2), Qid)

          call destroy(AA(k2))
          call destroy(BB(k2))
          call destroy(CC(k2))
      end do

      deallocate(AA, BB, CC, Dmap%elem, Strmap%elem)
      call destroy(DD)
      call destroy(ZalU)
      call destroy(Unonq)
      call destroy(Qid)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_first_site_qtensorc <<< <<<'

  end subroutine get_zaletel_1st_site_qtensorc

ZaletelOps_f90.extract_qop_from_op_qtensor()

fortran-subroutine - April 2016 (dj) Extract blocks according to quantum numbers from operator on complete local Hilbert space.

Arguments

QPropTYPE(qtensorc), out

On exit, the q matrix of Prop is generated according to the information given.

PropTYPE(tensorc), in

Matrix on without q-structure to be parsed.

DmapTYPE(vector_int), in

Vector containing the local dimensions of each block in the q-matrix (gained from identity).

StrmapTYPE(vector_int), in

Cumulative sum of the dimension of each block (gained from identity).

locdINTEGER, in

Dimension of the matrix Prop.

QRefTYPE(qtensor), in

If QRef has no entries (and therefore norm zero), QProp will be empty as well.

QidTYPE(qtensor), in

The identity matrix as q-matrix.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine extract_qop_from_op_qtensor(QProp, Prop, Dmap, Strmap, &
                                              locd, QRef, Qid)
      type(qtensorc), intent(out) :: QProp
      type(tensorc), intent(in) :: Prop
      type(vector_int), intent(in) :: Dmap, Strmap
      integer, intent(in) :: locd
      type(qtensor), intent(in) :: QRef, Qid

      ! Local variables
      ! ---------------

      ! for looping
      integer :: jj

      ! for indexing
      integer :: idx

      ! dimension of current block
      integer :: d1, d2

      ! total number of symmetries
      integer :: snqs

      ! temporary matrix for subblock
      type(tensorc) :: Tmp

      ! Prop as array
      complex(KIND=rKind), dimension(:, :), allocatable :: parr

      ! storing hash temporarily
      real(KIND=rKind) :: hsh

      !DBGprint *, '>>> >>> extract_qop_from_op_qtensor'

      call create(Qprop, Qid%nqs, Qid%nb)

      ! Fast return, leave propagator empty if reference has no entries
      if(Qref%nb == 0) return

      allocate(parr(Prop%dl(1), Prop%dl(2)))
      parr = reshape(Prop%elem(:Prop%dim), [Prop%dl(1), Prop%dl(2)])

      snqs = sum(Qid%nqs)

      do jj = 1, Qid%nb
          d2 = Dmap%elem(jj)

          hsh = prime_hash(sumq(Qid%Data(jj)%qq(snqs + 1:2 * snqs), &
                                QRef%Data(1)%qq(:snqs) - &
                                Qref%Data(1)%qq(snqs + 1:2 * snqs), &
                                Qid%nqs), [1], Qid%nqs)
      ! To-Do: Why last line? Check!

          idx = FindTagIndex(hsh, Qid%hash(1:Qid%nb))
          if(idx < 1) cycle

          d1 = Dmap%elem(idx)
          call create(Tmp, [d1, d2])
          Tmp%elem = reshape(&
               parr(Strmap%elem(idx) + 1:Strmap%elem(idx) + d1, &
                    Strmap%elem(jj) + 1:Strmap%elem(jj) + d2), [Tmp%dim])

          if(norm(Tmp) > 10.0_rKind * numzero) then
              Qprop%nb = Qprop%nb + 1
              allocate(Qprop%Data(QProp%nb)%qq(2 * snqs))
              QProp%Data(QProp%nb)%qq = [Qid%Data(idx)%qq(:snqs), &
                                         Qid%Data(jj)%qq(:snqs)]
              call set_hash(QProp, [2], QProp%nb)
              call create(QProp%Data(QProp%nb)%Tens, Tmp%dl)
              QProp%Data(QProp%nb)%Tens%elem = Tmp%elem
      ! To-do: Why does next line not work instead of two previous lines.
              !call copy(QProp%Data(QProp%nb)%Tens, Tmp)
          end if
          call destroy(Tmp)
      end do
      call skim(QProp, 10.0_rKind * numzero)

      deallocate(parr)

      !DBGprint *, 'extract_qop_from_op_qtensor <<< <<<'

  end subroutine extract_qop_from_op_qtensor

ZaletelOps_f90.extract_qop_from_op_qtensorc()

fortran-subroutine - April 2016 (dj) Extract blocks according to quantum numbers from operator on complete local Hilbert space.

Arguments

QPropTYPE(qtensorc), out

On exit, the q matrix of Prop is generated according to the information given.

PropTYPE(tensorc), in

Matrix on without q-structure to be parsed.

DmapTYPE(vector_int), in

Vector containing the local dimensions of each block in the q-matrix (gained from identity).

StrmapTYPE(vector_int), in

Cumulative sum of the dimension of each block (gained from identity).

locdINTEGER, in

Dimension of the matrix Prop.

QRefTYPE(qtensorc), in

If QRef has no entries (and therefore norm zero), QProp will be empty as well.

QidTYPE(qtensorc), in

The identity matrix as q-matrix.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine extract_qop_from_op_qtensorc(QProp, Prop, Dmap, Strmap, &
                                              locd, QRef, Qid)
      type(qtensorc), intent(out) :: QProp
      type(tensorc), intent(in) :: Prop
      type(vector_int), intent(in) :: Dmap, Strmap
      integer, intent(in) :: locd
      type(qtensorc), intent(in) :: QRef, Qid

      ! Local variables
      ! ---------------

      ! for looping
      integer :: jj

      ! for indexing
      integer :: idx

      ! dimension of current block
      integer :: d1, d2

      ! total number of symmetries
      integer :: snqs

      ! temporary matrix for subblock
      type(tensorc) :: Tmp

      ! Prop as array
      complex(KIND=rKind), dimension(:, :), allocatable :: parr

      ! storing hash temporarily
      real(KIND=rKind) :: hsh

      !DBGprint *, '>>> >>> extract_qop_from_op_qtensorc'

      call create(Qprop, Qid%nqs, Qid%nb)

      ! Fast return, leave propagator empty if reference has no entries
      if(Qref%nb == 0) return

      allocate(parr(Prop%dl(1), Prop%dl(2)))
      parr = reshape(Prop%elem(:Prop%dim), [Prop%dl(1), Prop%dl(2)])

      snqs = sum(Qid%nqs)

      do jj = 1, Qid%nb
          d2 = Dmap%elem(jj)

          hsh = prime_hash(sumq(Qid%Data(jj)%qq(snqs + 1:2 * snqs), &
                                QRef%Data(1)%qq(:snqs) - &
                                Qref%Data(1)%qq(snqs + 1:2 * snqs), &
                                Qid%nqs), [1], Qid%nqs)
      ! To-Do: Why last line? Check!

          idx = FindTagIndex(hsh, Qid%hash(1:Qid%nb))
          if(idx < 1) cycle

          d1 = Dmap%elem(idx)
          call create(Tmp, [d1, d2])
          Tmp%elem = reshape(&
               parr(Strmap%elem(idx) + 1:Strmap%elem(idx) + d1, &
                    Strmap%elem(jj) + 1:Strmap%elem(jj) + d2), [Tmp%dim])

          if(norm(Tmp) > 10.0_rKind * numzero) then
              Qprop%nb = Qprop%nb + 1
              allocate(Qprop%Data(QProp%nb)%qq(2 * snqs))
              QProp%Data(QProp%nb)%qq = [Qid%Data(idx)%qq(:snqs), &
                                         Qid%Data(jj)%qq(:snqs)]
              call set_hash(QProp, [2], QProp%nb)
              call create(QProp%Data(QProp%nb)%Tens, Tmp%dl)
              QProp%Data(QProp%nb)%Tens%elem = Tmp%elem
      ! To-do: Why does next line not work instead of two previous lines.
              !call copy(QProp%Data(QProp%nb)%Tens, Tmp)
          end if
          call destroy(Tmp)
      end do
      call skim(QProp, 10.0_rKind * numzero)

      deallocate(parr)

      !DBGprint *, 'extract_qop_from_op_qtensorc <<< <<<'

  end subroutine extract_qop_from_op_qtensorc

ZaletelOps_f90.get_mapping_qtensor()

fortran-subroutine - April 2016 (dj) Get local dimension and position of each block in a q-matrix.

Arguments

QidTYPE(qtensor), in

Identity matrix as q-matrix representation. It is diagonal and used to identity the dimension and position of each block. Order of blocks is taken as is.

DmapTYPE(vector_int), out

Contains the dimension of each block.

StrmapTYPE(vector_int), out

Contains the cumulated dimensions of the previous blocks. The ii-th block starts therefore at Strmap%elem(ii) + 1.

locdINTEGER, out

Adds up the dimension of each block gaining the dimension of the matrix without q-representation.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_mapping_qtensor(Qid, Dmap, Strmap, locd)
      type(qtensor), intent(in) :: Qid
      type(vector_int), intent(out) :: Dmap, Strmap
      integer, intent(out) :: locd

      ! Local variables
      ! ---------------

      integer :: kk

      !DBGprint *, '>>> >>> get_mapping_qtensor'

      allocate(Dmap%elem(Qid%nb), Strmap%elem(Qid%nb))
      Strmap%elem(1) = 0
      locd = 0

      do kk = 1, Qid%nb
          Dmap%elem(kk) = Qid%Data(kk)%Tens%dl(1)
          if(kk > 1) Strmap%elem(kk) = Strmap%elem(kk - 1) + Dmap%elem(kk - 1)
          locd = locd + Dmap%elem(kk)
      end do

      !DBGprint *, 'get_mapping_qtensor <<< <<<'

  end subroutine get_mapping_qtensor

ZaletelOps_f90.get_mapping_qtensorc()

fortran-subroutine - April 2016 (dj) Get local dimension and position of each block in a q-matrix.

Arguments

QidTYPE(qtensorc), in

Identity matrix as q-matrix representation. It is diagonal and used to identity the dimension and position of each block. Order of blocks is taken as is.

DmapTYPE(vector_int), out

Contains the dimension of each block.

StrmapTYPE(vector_int), out

Contains the cumulated dimensions of the previous blocks. The ii-th block starts therefore at Strmap%elem(ii) + 1.

locdINTEGER, out

Adds up the dimension of each block gaining the dimension of the matrix without q-representation.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_mapping_qtensorc(Qid, Dmap, Strmap, locd)
      type(qtensorc), intent(in) :: Qid
      type(vector_int), intent(out) :: Dmap, Strmap
      integer, intent(out) :: locd

      ! Local variables
      ! ---------------

      integer :: kk

      !DBGprint *, '>>> >>> get_mapping_qtensorc'

      allocate(Dmap%elem(Qid%nb), Strmap%elem(Qid%nb))
      Strmap%elem(1) = 0
      locd = 0

      do kk = 1, Qid%nb
          Dmap%elem(kk) = Qid%Data(kk)%Tens%dl(1)
          if(kk > 1) Strmap%elem(kk) = Strmap%elem(kk - 1) + Dmap%elem(kk - 1)
          locd = locd + Dmap%elem(kk)
      end do

      !DBGprint *, 'get_mapping_qtensorc <<< <<<'

  end subroutine get_mapping_qtensorc

ZaletelOps_f90.get_zaletel_last_site_tensor()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the last (right) site of the system for a qMPO.

Arguments

ZalTYPE(sr_matrix_tensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of last (right) site.

HamTYPE(sr_matrix_tensor), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the last site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_last_site_tensor(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_tensorc), intent(out) :: Zal
      type(sr_matrix_tensor), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ===============

      ! for looping / indexing
      integer :: k1, k2

      ! left and right chi of the Zaletel MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      ! Different parts in MPO
      type(tensorc) :: DD
      type(tensorc), pointer :: AA(:), BB(:), CC(:)

      ! Temporay matrix
      type(tensorc) :: Tmp

      ! zaletel exponential
      type(tensorc) :: ZalU

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_last_site_tensor'

      call build_hcb_idb(Hcb, Idb)

      ! Get bond dimensions
      chil = Ham%rbd - 1
      chir = 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      k1 = 1
      allocate(Zal%row(chil))

      nel = 1
      do k1 = 1, chil
          Zal%Row(k1)%numel = nel
          allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
          Zal%Row(k1)%ind = 1
      end do

      locd = Ham%Row(1)%Op(1)%dl(1)

      ! Copy/set terms from Ham
      ! -----------------------

      k1 = chil + 1
      k2 = 1
      call copy(DD, Ham%row(k1)%Op(k2))

      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))
      do k1 = 1, chil - 1
          ! Create empty AA, CC
          call create(AA(k1), DD%dl, init='0')
          call create(CC(k1), DD%dl, init='0')

          ! Copy BB
          call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
      end do

      ! W_D term is just exp(D)
      ! -----------------------

      k1 = 1
      k2 = 1
      call create(Zal%row(k1)%Op(k2), DD%dl)
      call copy(Tmp, DD)

      if(mpo_is_hermitian) then
          call expmh(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      else
          call expm(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      end if
      call destroy(Tmp)
      call create(ZalU, [4 * locd, 4 * locd], init='0')

      do k1 = 1, chil - 1
          call get_zaletel_expm(ZalU, dt, AA(k1), BB(k1), CC(k1), DD, Hcb, &
                               Idb, locd)
          call create(Zal%Row(k1 + 1)%Op(k2), DD%dl)
          call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(k2), &
                                       ZalU, 1, 0, locd)

          call destroy(AA(k1))
          call destroy(BB(k1))
          call destroy(CC(k1))
      end do

      call destroy(DD)
      call destroy(ZalU)
      deallocate(AA, BB, CC)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_last_site_tensor <<< <<<'

  end subroutine get_zaletel_last_site_tensor

ZaletelOps_f90.get_zaletel_last_site_tensorc()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the last (right) site of the system for a qMPO.

Arguments

ZalTYPE(sr_matrix_tensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of last (right) site.

HamTYPE(sr_matrix_tensorc), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the last site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_last_site_tensorc(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_tensorc), intent(out) :: Zal
      type(sr_matrix_tensorc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ===============

      ! for looping / indexing
      integer :: k1, k2

      ! left and right chi of the Zaletel MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      ! Different parts in MPO
      type(tensorc) :: DD
      type(tensorc), pointer :: AA(:), BB(:), CC(:)

      ! Temporay matrix
      type(tensorc) :: Tmp

      ! zaletel exponential
      type(tensorc) :: ZalU

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_last_site_tensorc'

      call build_hcb_idb(Hcb, Idb)

      ! Get bond dimensions
      chil = Ham%rbd - 1
      chir = 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      k1 = 1
      allocate(Zal%row(chil))

      nel = 1
      do k1 = 1, chil
          Zal%Row(k1)%numel = nel
          allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
          Zal%Row(k1)%ind = 1
      end do

      locd = Ham%Row(1)%Op(1)%dl(1)

      ! Copy/set terms from Ham
      ! -----------------------

      k1 = chil + 1
      k2 = 1
      call copy(DD, Ham%row(k1)%Op(k2))

      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))
      do k1 = 1, chil - 1
          ! Create empty AA, CC
          call create(AA(k1), DD%dl, init='0')
          call create(CC(k1), DD%dl, init='0')

          ! Copy BB
          call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
      end do

      ! W_D term is just exp(D)
      ! -----------------------

      k1 = 1
      k2 = 1
      call create(Zal%row(k1)%Op(k2), DD%dl)
      call copy(Tmp, DD)

      if(mpo_is_hermitian) then
          call expmh(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      else
          call expm(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      end if
      call destroy(Tmp)
      call create(ZalU, [4 * locd, 4 * locd], init='0')

      do k1 = 1, chil - 1
          call get_zaletel_expm(ZalU, dt, AA(k1), BB(k1), CC(k1), DD, Hcb, &
                               Idb, locd)
          call create(Zal%Row(k1 + 1)%Op(k2), DD%dl)
          call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(k2), &
                                       ZalU, 1, 0, locd)

          call destroy(AA(k1))
          call destroy(BB(k1))
          call destroy(CC(k1))
      end do

      call destroy(DD)
      call destroy(ZalU)
      deallocate(AA, BB, CC)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_last_site_tensorc <<< <<<'

  end subroutine get_zaletel_last_site_tensorc

ZaletelOps_f90.get_zaletel_last_site_qtensor()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the last (right) site of the system.

Arguments

ZalTYPE(sr_matrix_qtensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of last (right) site.

HamTYPE(sr_matrix_qtensor), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the last site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_last_site_qtensor(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_qtensorc), intent(out) :: Zal
      type(sr_matrix_qtensor), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! bond dimensions of MPO
      integer :: chil, chir

      ! local dimension (of complete Hilbert space)
      integer :: locd

      ! for indexing / looping
      integer :: k1, k2, jj, idx

      ! number of elements in sparse row
      integer :: nel

      ! for mapping between matrices and q-matrices
      type(vector_int) :: Dmap, Strmap

      ! Different parts in MPO
      type(qtensorc) :: DD
      type(qtensorc), pointer :: AA(:), BB(:), CC(:)

      ! Zaletel exponential and propagator for local Hilbert space
      type(tensorc) :: ZalU, Unonq

      ! ID as q-matrix / temporary q-matrix
      type(qtensorc) :: Qid, Tmp

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_last_site_qtensor'

      call build_hcb_idb(Hcb, Idb)

      chil = Ham%rbd - 1
      chir = 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      allocate(Zal%Row(chil))
      nel = 1
      do k1 = 1, chil
          Zal%Row(k1)%numel = nel
          allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
          Zal%Row(k1)%ind = 1
      end do

      ! Copy/set terms from Ham
      ! -----------------------

      call copy(Qid, Ham%Row(1)%Op(Ham%cbd))

      k1 = chil + 1
      k2 = 1
      call copy(DD, Ham%Row(k1)%Op(k2))

      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))
      do k1 = 1, chil - 1
          ! Create empty AA, CC
          call create(AA(k1), Qid%nqs, 1)
          call create(CC(k1), Qid%nqs, 1)

          ! Copy BB
          call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
      end do

      ! W_D term is just exp(D) and D is blockdiagonal
      ! ----------------------------------------------

      ! Get the blocks, set to 0 in case local terms does not access
      k1 = 1
      k2 = 1
      call copy(Zal%Row(k1)%Op(k2), Qid)
      call copy(Tmp, DD)

      ! Exponentiate
      do jj = 1, Tmp%nb
          idx = FindTagIndex(Tmp%hash(jj), &
                             Zal%Row(k1)%Op(k2)%hash(1:Zal%Row(k1)%Op(k2)%nb))
          if(mpo_is_hermitian) then
              call expmh(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                        dt, Tmp%Data(jj)%Tens, 1)
          else
              call expm(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                       dt, Tmp%Data(jj)%Tens, 1)
          end if
      end do

      call destroy(Tmp)
      call skim(Zal%Row(k1)%Op(k2), 10.0_rKind * numzero)

      ! Other terms W_X (cannot be created as q-matrices)
      ! ---------------

      ! Get mapping from/to q-matrices
      call get_mapping_qtensorc(Qid, Dmap, Strmap, locd)

      call create(ZalU, [4 * locd, 4 * locd], init='0')
      call create(Unonq, [locd, locd], init='0')

      do k1 = 1, chil - 1
          call get_zaletel_expm(ZalU, dt, AA(k1), BB(k1), CC(k1), &
                               DD, Hcb, Idb, locd, Dmap, Strmap, &
                               Qid)
          call extract_op_from_zaletel(Unonq, ZalU, 1, 0, locd)
          call extract_qop_from_op(Zal%Row(k1 + 1)%Op(k2), Unonq, Dmap, &
                                   Strmap, locd, BB(k1), Qid)

          call destroy(AA(k1))
          call destroy(BB(k1))
          call destroy(CC(k1))
      end do

      deallocate(AA, BB, CC, Dmap%elem, Strmap%elem)
      call destroy(DD)
      call destroy(ZalU)
      call destroy(Unonq)
      call destroy(Qid)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_last_site_qtensor <<< <<<'

  end subroutine get_zaletel_last_site_qtensor

ZaletelOps_f90.get_zaletel_last_site_qtensorc()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the last (right) site of the system.

Arguments

ZalTYPE(sr_matrix_qtensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of last (right) site.

HamTYPE(sr_matrix_qtensorc), in

MPO-matrix (rank-4 tensor) of the Hamiltonian of the last site.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_last_site_qtensorc(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_qtensorc), intent(out) :: Zal
      type(sr_matrix_qtensorc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! bond dimensions of MPO
      integer :: chil, chir

      ! local dimension (of complete Hilbert space)
      integer :: locd

      ! for indexing / looping
      integer :: k1, k2, jj, idx

      ! number of elements in sparse row
      integer :: nel

      ! for mapping between matrices and q-matrices
      type(vector_int) :: Dmap, Strmap

      ! Different parts in MPO
      type(qtensorc) :: DD
      type(qtensorc), pointer :: AA(:), BB(:), CC(:)

      ! Zaletel exponential and propagator for local Hilbert space
      type(tensorc) :: ZalU, Unonq

      ! ID as q-matrix / temporary q-matrix
      type(qtensorc) :: Qid, Tmp

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_last_site_qtensorc'

      call build_hcb_idb(Hcb, Idb)

      chil = Ham%rbd - 1
      chir = 1
      Zal%rbd = chil
      Zal%cbd = chir

      ! Allocate space
      allocate(Zal%Row(chil))
      nel = 1
      do k1 = 1, chil
          Zal%Row(k1)%numel = nel
          allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
          Zal%Row(k1)%ind = 1
      end do

      ! Copy/set terms from Ham
      ! -----------------------

      call copy(Qid, Ham%Row(1)%Op(Ham%cbd))

      k1 = chil + 1
      k2 = 1
      call copy(DD, Ham%Row(k1)%Op(k2))

      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))
      do k1 = 1, chil - 1
          ! Create empty AA, CC
          call create(AA(k1), Qid%nqs, 1)
          call create(CC(k1), Qid%nqs, 1)

          ! Copy BB
          call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
      end do

      ! W_D term is just exp(D) and D is blockdiagonal
      ! ----------------------------------------------

      ! Get the blocks, set to 0 in case local terms does not access
      k1 = 1
      k2 = 1
      call copy(Zal%Row(k1)%Op(k2), Qid)
      call copy(Tmp, DD)

      ! Exponentiate
      do jj = 1, Tmp%nb
          idx = FindTagIndex(Tmp%hash(jj), &
                             Zal%Row(k1)%Op(k2)%hash(1:Zal%Row(k1)%Op(k2)%nb))
          if(mpo_is_hermitian) then
              call expmh(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                        dt, Tmp%Data(jj)%Tens, 1)
          else
              call expm(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                       dt, Tmp%Data(jj)%Tens, 1)
          end if
      end do

      call destroy(Tmp)
      call skim(Zal%Row(k1)%Op(k2), 10.0_rKind * numzero)

      ! Other terms W_X (cannot be created as q-matrices)
      ! ---------------

      ! Get mapping from/to q-matrices
      call get_mapping_qtensorc(Qid, Dmap, Strmap, locd)

      call create(ZalU, [4 * locd, 4 * locd], init='0')
      call create(Unonq, [locd, locd], init='0')

      do k1 = 1, chil - 1
          call get_zaletel_expm(ZalU, dt, AA(k1), BB(k1), CC(k1), &
                               DD, Hcb, Idb, locd, Dmap, Strmap, &
                               Qid)
          call extract_op_from_zaletel(Unonq, ZalU, 1, 0, locd)
          call extract_qop_from_op(Zal%Row(k1 + 1)%Op(k2), Unonq, Dmap, &
                                   Strmap, locd, BB(k1), Qid)

          call destroy(AA(k1))
          call destroy(BB(k1))
          call destroy(CC(k1))
      end do

      deallocate(AA, BB, CC, Dmap%elem, Strmap%elem)
      call destroy(DD)
      call destroy(ZalU)
      call destroy(Unonq)
      call destroy(Qid)
      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_last_site_qtensorc <<< <<<'

  end subroutine get_zaletel_last_site_qtensorc

ZaletelOps_f90.get_zaletel_bulk_site_tensor()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the bulk sites of the system.

Arguments

ZalTYPE(sr_matrix_tensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of site.

HamTYPE(sr_matrix_tensor), in

MPO-matrix (rank-4 tensor) of the Hamiltonian.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_bulk_site_tensor(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_tensorc), intent(out) :: Zal
      type(sr_matrix_tensor), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! For looping / indexing
      integer :: k1, k2

      ! Bond dimension in Zaletel MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      type(tensorc) :: ZalU

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      ! Matrices of MPO in Zaletel notation
      type(tensorc) :: DD
      type(tensorc), pointer :: AA(:), BB(:), CC(:)

      ! Temporary matrices
      type(tensorc) :: Tmp, TmpA, TmpB

      !DBGprint *, '>>> >>> get_zaletel_bulk_site_tensor'

      call build_hcb_idb(Hcb, Idb)

      chil = Ham%rbd - 1
      chir = Ham%cbd - 1

      ! Copy element DD
      call copy(DD, Ham%Row(Ham%rbd)%Op(1))
      locd = DD%dl(1)

      ! Take care of AA, BB, CC
      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))

      do k1 = 1, (chil - 1)
          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              call copy(AA(k1), Ham%Row(k1 + 1)%Op(1))
              call create(BB(k1), DD%dl, init='0')
          else
              call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
              if(Ham%row(k1 + 1)%numel > 1) then
                  call copy(AA(k1), Ham%Row(k1+1)%Op(2))
              else
                  call create(AA(k1), DD%dl, init='0')
              end if
          end if
      end do

      !do k2 = 1, chil - 1!Ham%Row(Ham%rbd)%numel - 2
      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call copy(CC(k2), Ham%Row(Ham%rbd)%Op(k2 + 1))
      end do

      ! Allocate space for rows
      Zal%rbd = chil
      Zal%cbd = chir
      allocate(Zal%Row(chil))

      ! First row
      ! ---------

      k1 = 1
      nel = Ham%Row(Ham%rbd)%numel - 1
      Zal%Row(k1)%numel = nel
      allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
      Zal%Row(k1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

      ! D-term is element (1, 1)
      k2 = 1
      call copy(Tmp, DD)
      call create(Zal%Row(k1)%Op(k2), DD%dl, init='0')
      if(mpo_is_hermitian) then
          call expmh(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      else
          call expm(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      end if
      call destroy(Tmp)

      call create(ZalU, [4 * locd, 4 * locd], init='0')

      ! Remaining part of the first row (1, :)
      call create(Tmp, DD%dl, init='0')
      do k2 = 1, (nel - 1)
          call get_zaletel_expm(ZalU, dt, Tmp, BB(k2), CC(k2), DD, Hcb, &
                               Idb, locd)
          call create(Zal%Row(k1)%Op(k2 + 1), DD%dl)
          call extract_op_from_zaletel(Zal%Row(k1)%Op(k2 + 1), &
                                       ZalU, 0, 1, locd)
      end do
      call destroy(Tmp)

      ! Everything apart from first row
      ! -------------------------------

      rows: do k1 = 1, chil - 1
          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              nel = 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(k1 + 1)%ind(1)

              call copy(Tmp, Ham%Row(k1 + 1)%Op(1))
              call create(TmpA, DD%dl, init='0')
              call get_zaletel_expm(ZalU, dt, Tmp, TmpA, TmpA, DD, Hcb, &
                                   Idb, locd)
              call create(Zal%Row(k1 + 1)%Op(1), DD%dl)
              call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(1), ZalU, &
                                           1, 1, locd)

              call destroy(Tmp)
              call destroy(TmpA)
          else
              nel = Ham%Row(Ham%rbd)%numel - 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

              ! First column
              call create(Tmp, DD%dl, init='0')
              call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), Tmp, DD, Hcb, &
                                   Idb, locd)
              call create(Zal%Row(k1 + 1)%Op(1), DD%dl)
              call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(1), ZalU, &
                                           1, 0, locd)
              call destroy(Tmp)

              ! Remaining columns
              cols: do k2 = 1, (nel - 1)

                  if(k1 == k2) then
                      call copy(Tmp, AA(k1))
                  else
                      call create(Tmp, DD%dl, init='0')
                  end if

                  call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), CC(k2), DD, &
                                       Hcb, Idb, locd)

                  ! Calculate commutator between terms
      ! To-Do: commutator or anti-commutator
                  call contr(TmpA, BB(k1), CC(k2), [2], [1])
                  call contr(TmpB, CC(k2), BB(k1), [2], [1])
                  call gaxpy(TmpA, done, TmpB)
                  if(norm(TmpA) < numzero) then
                      call destroy(TmpA)
                      call copy(TmpA, Tmp)
                  end if
                  call create(Zal%Row(k1 + 1)%Op(k2 + 1), DD%dl)
                  call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(k2 + 1), &
                                               ZalU, 1, 1, locd)

                  call destroy(Tmp)
                  call destroy(TmpA)
                  call destroy(TmpB)
              end do cols
          end if
      end do rows

      call destroy(DD)

      do k1 = 1, chil - 1
          call destroy(AA(k1))
          call destroy(BB(k1))
      end do

      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call destroy(CC(k2))
      end do
      deallocate(AA, BB, CC)

      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_bulk_site_tensor <<< <<<'

  end subroutine get_zaletel_bulk_site_tensor

ZaletelOps_f90.get_zaletel_bulk_site_tensorc()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the bulk sites of the system.

Arguments

ZalTYPE(sr_matrix_tensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of site.

HamTYPE(sr_matrix_tensorc), in

MPO-matrix (rank-4 tensor) of the Hamiltonian.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_bulk_site_tensorc(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_tensorc), intent(out) :: Zal
      type(sr_matrix_tensorc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! For looping / indexing
      integer :: k1, k2

      ! Bond dimension in Zaletel MPO
      integer :: chil, chir

      ! local dimension (local Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      type(tensorc) :: ZalU

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      ! Matrices of MPO in Zaletel notation
      type(tensorc) :: DD
      type(tensorc), pointer :: AA(:), BB(:), CC(:)

      ! Temporary matrices
      type(tensorc) :: Tmp, TmpA, TmpB

      !DBGprint *, '>>> >>> get_zaletel_bulk_site_tensorc'

      call build_hcb_idb(Hcb, Idb)

      chil = Ham%rbd - 1
      chir = Ham%cbd - 1

      ! Copy element DD
      call copy(DD, Ham%Row(Ham%rbd)%Op(1))
      locd = DD%dl(1)

      ! Take care of AA, BB, CC
      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))

      do k1 = 1, (chil - 1)
          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              call copy(AA(k1), Ham%Row(k1 + 1)%Op(1))
              call create(BB(k1), DD%dl, init='0')
          else
              call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
              if(Ham%row(k1 + 1)%numel > 1) then
                  call copy(AA(k1), Ham%Row(k1+1)%Op(2))
              else
                  call create(AA(k1), DD%dl, init='0')
              end if
          end if
      end do

      !do k2 = 1, chil - 1!Ham%Row(Ham%rbd)%numel - 2
      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call copy(CC(k2), Ham%Row(Ham%rbd)%Op(k2 + 1))
      end do

      ! Allocate space for rows
      Zal%rbd = chil
      Zal%cbd = chir
      allocate(Zal%Row(chil))

      ! First row
      ! ---------

      k1 = 1
      nel = Ham%Row(Ham%rbd)%numel - 1
      Zal%Row(k1)%numel = nel
      allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
      Zal%Row(k1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

      ! D-term is element (1, 1)
      k2 = 1
      call copy(Tmp, DD)
      call create(Zal%Row(k1)%Op(k2), DD%dl, init='0')
      if(mpo_is_hermitian) then
          call expmh(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      else
          call expm(Zal%Row(k1)%Op(k2), dt, Tmp, 1)
      end if
      call destroy(Tmp)

      call create(ZalU, [4 * locd, 4 * locd], init='0')

      ! Remaining part of the first row (1, :)
      call create(Tmp, DD%dl, init='0')
      do k2 = 1, (nel - 1)
          call get_zaletel_expm(ZalU, dt, Tmp, BB(k2), CC(k2), DD, Hcb, &
                               Idb, locd)
          call create(Zal%Row(k1)%Op(k2 + 1), DD%dl)
          call extract_op_from_zaletel(Zal%Row(k1)%Op(k2 + 1), &
                                       ZalU, 0, 1, locd)
      end do
      call destroy(Tmp)

      ! Everything apart from first row
      ! -------------------------------

      rows: do k1 = 1, chil - 1
          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              nel = 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(k1 + 1)%ind(1)

              call copy(Tmp, Ham%Row(k1 + 1)%Op(1))
              call create(TmpA, DD%dl, init='0')
              call get_zaletel_expm(ZalU, dt, Tmp, TmpA, TmpA, DD, Hcb, &
                                   Idb, locd)
              call create(Zal%Row(k1 + 1)%Op(1), DD%dl)
              call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(1), ZalU, &
                                           1, 1, locd)

              call destroy(Tmp)
              call destroy(TmpA)
          else
              nel = Ham%Row(Ham%rbd)%numel - 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

              ! First column
              call create(Tmp, DD%dl, init='0')
              call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), Tmp, DD, Hcb, &
                                   Idb, locd)
              call create(Zal%Row(k1 + 1)%Op(1), DD%dl)
              call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(1), ZalU, &
                                           1, 0, locd)
              call destroy(Tmp)

              ! Remaining columns
              cols: do k2 = 1, (nel - 1)

                  if(k1 == k2) then
                      call copy(Tmp, AA(k1))
                  else
                      call create(Tmp, DD%dl, init='0')
                  end if

                  call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), CC(k2), DD, &
                                       Hcb, Idb, locd)

                  ! Calculate commutator between terms
      ! To-Do: commutator or anti-commutator
                  call contr(TmpA, BB(k1), CC(k2), [2], [1])
                  call contr(TmpB, CC(k2), BB(k1), [2], [1])
                  call gaxpy(TmpA, done, TmpB)
                  if(norm(TmpA) < numzero) then
                      call destroy(TmpA)
                      call copy(TmpA, Tmp)
                  end if
                  call create(Zal%Row(k1 + 1)%Op(k2 + 1), DD%dl)
                  call extract_op_from_zaletel(Zal%Row(k1 + 1)%Op(k2 + 1), &
                                               ZalU, 1, 1, locd)

                  call destroy(Tmp)
                  call destroy(TmpA)
                  call destroy(TmpB)
              end do cols
          end if
      end do rows

      call destroy(DD)

      do k1 = 1, chil - 1
          call destroy(AA(k1))
          call destroy(BB(k1))
      end do

      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call destroy(CC(k2))
      end do
      deallocate(AA, BB, CC)

      call destroy(Hcb)
      call destroy(Idb)

      !DBGprint *, 'get_zaletel_bulk_site_tensorc <<< <<<'

  end subroutine get_zaletel_bulk_site_tensorc

ZaletelOps_f90.get_zaletel_bulk_site_qtensor()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the bulk sites of the system for symmetry conserving qMPOs.

Arguments

ZalTYPE(sr_matrix_qtensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of site.

HamTYPE(sr_matrix_qtensor), in

MPO-matrix (rank-4 tensor) of the Hamiltonian.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_bulk_site_qtensor(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_qtensorc), intent(out) :: Zal
      type(sr_Matrix_qtensor), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! bond dimension in Zaletel MPO
      integer :: chil, chir

      ! local dimension (of complete Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      ! for indexing / looping
      integer :: k1, k2, jj, idx

      ! for mapping between matrices and q-matrices
      type(vector_int) :: Dmap, Strmap

      ! Different parts in MPO
      type(qtensorc) :: DD
      type(qtensorc), pointer :: AA(:), BB(:), CC(:)

      ! Zaletel exponential and propagator for local Hilbert space
      type(tensorc) :: ZalU, Unonq

      ! ID as q-matrix / temporary q-matrices
      type(qtensorc) :: Qid, Tmp, TmpA, TmpB

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_bulk_site_qtensor'

      call build_hcb_idb(Hcb, Idb)

      ! Copy identity and get mapping between q-matrices and matrices
      call copy(Qid, Ham%Row(1)%Op(1))
      call get_mapping_qtensorc(Qid, Dmap, Strmap, locd)
      call create(ZalU, [4 * locd, 4 * locd], init='0')
      call create(Unonq, [locd, locd], init='0')

      chil = Ham%rbd - 1
      chir = Ham%cbd - 1

      ! Copy element DD
      call copy(DD, Ham%Row(Ham%rbd)%Op(1))

      ! Take care of AA, BB, CC
      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))

      do k1 = 1, chil - 1
          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              call copy(AA(k1), Ham%Row(k1 + 1)%Op(1))
              call create(BB(k1), Qid%nqs, 1)
          else
              call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
              if(Ham%Row(k1 + 1)%numel > 1) then
                  call copy(AA(k1), Ham%Row(k1 + 1)%Op(2))
              else
                  call create(AA(k1), Qid%nqs, 1)
              end if
          end if
      end do

      !do k2 = 1, chil - 1
      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call copy(CC(k2), Ham%Row(Ham%rbd)%Op(k2 + 1))
      end do

      ! Allocate space for rows
      Zal%rbd = chil
      Zal%cbd = chir
      allocate(Zal%Row(chil))

      ! First row
      ! ---------

      k1 = 1
      nel = Ham%Row(Ham%rbd)%numel - 1
      Zal%row(k1)%numel = nel
      allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
      Zal%Row(k1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

      ! D-term is element (1, 1)
      k2 = 1
      call copy(Tmp, DD)
      call copy(Zal%Row(k1)%Op(k2), Qid)

      do jj = 1, Tmp%nb
          idx = FindTagIndex(Tmp%hash(jj), &
                             Zal%Row(k1)%Op(k2)%hash(1:Zal%Row(k1)%Op(k2)%nb))
          if(mpo_is_hermitian) then
              call expmh(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                        dt, Tmp%Data(jj)%Tens, 1)
          else
              call expm(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                       dt, Tmp%Data(jj)%Tens, 1)
          end if
      end do
      call destroy(Tmp)
      call skim(Zal%Row(k1)%Op(k2), 10.0_rKind * numzero)

      ! Remaining part of first row (1, :)
      call create(Tmp, Qid%nqs, 1)

      do k2 = 1, (nel - 1)
          call get_zaletel_expm(ZalU, dt, Tmp, BB(k2), CC(k2), DD, Hcb, Idb, &
                               locd, Dmap, Strmap, Qid)
          call extract_op_from_zaletel(Unonq, ZalU, 0, 1, locd)
          call extract_qop_from_op(Zal%Row(k1)%Op(k2 + 1), Unonq, Dmap, &
                                   Strmap, locd, CC(k2), Qid)
          call skim(Zal%Row(k1)%Op(k2 + 1), 10.0_rKind * numzero)
      end do

      call destroy(Tmp)

      ! Everything apart from the first row
      ! -----------------------------------

      rows: do k1 = 1, chil - 1

          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              nel = 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(k1 + 1)%ind(1)

              call copy(Tmp, Ham%Row(k1 + 1)%Op(1))
              call create(TmpA, Qid%nqs, 1)

              call get_zaletel_expm(ZalU, dt, Tmp, TmpA, TmpA, DD, Hcb, Idb, &
                                   locd, Dmap, Strmap, Qid)
              call extract_op_from_zaletel(Unonq, ZalU, 1, 1, locd)
              call extract_qop_from_op(Zal%Row(k1 + 1)%Op(1), Unonq, Dmap, &
                                       Strmap, locd, Qid, Qid)
              call skim(Zal%Row(k1 + 1)%Op(1), 10.0_rKind * numzero)

              call destroy(Tmp)
              call destroy(TmpA)
          else
              nel = Ham%Row(Ham%rbd)%numel - 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

              ! First column
              call create(Tmp, Qid%nqs, 1)
              call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), Tmp, DD, Hcb, Idb, &
                                   locd, Dmap, Strmap, Qid)
              call extract_op_from_zaletel(Unonq, ZalU, 1, 0, locd)
              call extract_qop_from_op(Zal%Row(k1 + 1)%Op(1), Unonq, &
                                       Dmap, Strmap, locd, BB(k1), Qid)
              call skim(Zal%Row(k1 + 1)%Op(1), 10.0_rKind * numzero)
              call destroy(Tmp)

              ! Remaining columns
              cols: do k2 = 1, (nel - 1)
                  if(k1 == k2) then
                      call copy(Tmp, AA(k1))
                  else
                      call create(Tmp, Qid%nqs, 1)
                  end if

                  call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), CC(k2), DD, &
                                       Hcb, Idb, locd, Dmap, Strmap, &
                                       Qid)

                  call extract_op_from_zaletel(Unonq, ZalU, 1, 1, locd)

                  ! Calculate commutators between terms
      ! To-do: commutator or anti-commutator
                  call contr(TmpA, BB(k1), CC(k2), [2], [1])
                  call contr(TmpB, CC(k2), BB(k1), [2], [1])
                  call gaxpy(TmpA, done, TmpB)

                  if(TmpA%nb == 0) then
                      call destroy(TmpA)
                      call copy(TmpA, Tmp)
                  end if

                  call extract_qop_from_op(Zal%Row(k1 + 1)%Op(k2 + 1), &
                                           Unonq, Dmap, Strmap, locd, &
                                           TmpA, Qid)
                  call skim(Zal%Row(k1 + 1)%Op(k2 + 1), 10.0_rKind * numzero)

                  call destroy(Tmp)
                  call destroy(TmpA)
                  call destroy(TmpB)
              end do cols
          end if
      end do rows

      call destroy(DD)

      do k1 = 1, chil - 1
          call destroy(AA(k1))
          call destroy(BB(k1))
      end do

      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call destroy(CC(k2))
      end do
      deallocate(AA, BB, CC)

      call destroy(Hcb)
      call destroy(Idb)
      call destroy(Qid)

      !DBGprint *, 'get_zaletel_bulk_site_qtensor <<< <<<'

  end subroutine get_zaletel_bulk_site_qtensor

ZaletelOps_f90.get_zaletel_bulk_site_qtensorc()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for the bulk sites of the system for symmetry conserving qMPOs.

Arguments

ZalTYPE(sr_matrix_qtensorc), out

On exit, MPO-matrix (rank-4 tensor) for time evolution of site.

HamTYPE(sr_matrix_qtensorc), in

MPO-matrix (rank-4 tensor) of the Hamiltonian.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_bulk_site_qtensorc(Zal, Ham, dt, mpo_is_hermitian)
      type(sr_matrix_qtensorc), intent(out) :: Zal
      type(sr_Matrix_qtensorc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! bond dimension in Zaletel MPO
      integer :: chil, chir

      ! local dimension (of complete Hilbert space)
      integer :: locd

      ! number of elementes
      integer :: nel

      ! for indexing / looping
      integer :: k1, k2, jj, idx

      ! for mapping between matrices and q-matrices
      type(vector_int) :: Dmap, Strmap

      ! Different parts in MPO
      type(qtensorc) :: DD
      type(qtensorc), pointer :: AA(:), BB(:), CC(:)

      ! Zaletel exponential and propagator for local Hilbert space
      type(tensorc) :: ZalU, Unonq

      ! ID as q-matrix / temporary q-matrices
      type(qtensorc) :: Qid, Tmp, TmpA, TmpB

      ! 2x2 operators
      type(tensor) :: Hcb, Idb

      !DBGprint *, '>>> >>> get_zaletel_bulk_site_qtensorc'

      call build_hcb_idb(Hcb, Idb)

      ! Copy identity and get mapping between q-matrices and matrices
      call copy(Qid, Ham%Row(1)%Op(1))
      call get_mapping_qtensorc(Qid, Dmap, Strmap, locd)
      call create(ZalU, [4 * locd, 4 * locd], init='0')
      call create(Unonq, [locd, locd], init='0')

      chil = Ham%rbd - 1
      chir = Ham%cbd - 1

      ! Copy element DD
      call copy(DD, Ham%Row(Ham%rbd)%Op(1))

      ! Take care of AA, BB, CC
      allocate(AA(chil - 1), BB(chil - 1), CC(chil - 1))

      do k1 = 1, chil - 1
          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              call copy(AA(k1), Ham%Row(k1 + 1)%Op(1))
              call create(BB(k1), Qid%nqs, 1)
          else
              call copy(BB(k1), Ham%Row(k1 + 1)%Op(1))
              if(Ham%Row(k1 + 1)%numel > 1) then
                  call copy(AA(k1), Ham%Row(k1 + 1)%Op(2))
              else
                  call create(AA(k1), Qid%nqs, 1)
              end if
          end if
      end do

      !do k2 = 1, chil - 1
      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call copy(CC(k2), Ham%Row(Ham%rbd)%Op(k2 + 1))
      end do

      ! Allocate space for rows
      Zal%rbd = chil
      Zal%cbd = chir
      allocate(Zal%Row(chil))

      ! First row
      ! ---------

      k1 = 1
      nel = Ham%Row(Ham%rbd)%numel - 1
      Zal%row(k1)%numel = nel
      allocate(Zal%Row(k1)%ind(nel), Zal%Row(k1)%Op(nel))
      Zal%Row(k1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

      ! D-term is element (1, 1)
      k2 = 1
      call copy(Tmp, DD)
      call copy(Zal%Row(k1)%Op(k2), Qid)

      do jj = 1, Tmp%nb
          idx = FindTagIndex(Tmp%hash(jj), &
                             Zal%Row(k1)%Op(k2)%hash(1:Zal%Row(k1)%Op(k2)%nb))
          if(mpo_is_hermitian) then
              call expmh(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                        dt, Tmp%Data(jj)%Tens, 1)
          else
              call expm(Zal%Row(k1)%Op(k2)%Data(idx)%Tens, &
                       dt, Tmp%Data(jj)%Tens, 1)
          end if
      end do
      call destroy(Tmp)
      call skim(Zal%Row(k1)%Op(k2), 10.0_rKind * numzero)

      ! Remaining part of first row (1, :)
      call create(Tmp, Qid%nqs, 1)

      do k2 = 1, (nel - 1)
          call get_zaletel_expm(ZalU, dt, Tmp, BB(k2), CC(k2), DD, Hcb, Idb, &
                               locd, Dmap, Strmap, Qid)
          call extract_op_from_zaletel(Unonq, ZalU, 0, 1, locd)
          call extract_qop_from_op(Zal%Row(k1)%Op(k2 + 1), Unonq, Dmap, &
                                   Strmap, locd, CC(k2), Qid)
          call skim(Zal%Row(k1)%Op(k2 + 1), 10.0_rKind * numzero)
      end do

      call destroy(Tmp)

      ! Everything apart from the first row
      ! -----------------------------------

      rows: do k1 = 1, chil - 1

          if(Ham%Row(k1 + 1)%ind(1) > 1) then
              nel = 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(k1 + 1)%ind(1)

              call copy(Tmp, Ham%Row(k1 + 1)%Op(1))
              call create(TmpA, Qid%nqs, 1)

              call get_zaletel_expm(ZalU, dt, Tmp, TmpA, TmpA, DD, Hcb, Idb, &
                                   locd, Dmap, Strmap, Qid)
              call extract_op_from_zaletel(Unonq, ZalU, 1, 1, locd)
              call extract_qop_from_op(Zal%Row(k1 + 1)%Op(1), Unonq, Dmap, &
                                       Strmap, locd, Qid, Qid)
              call skim(Zal%Row(k1 + 1)%Op(1), 10.0_rKind * numzero)

              call destroy(Tmp)
              call destroy(TmpA)
          else
              nel = Ham%Row(Ham%rbd)%numel - 1
              Zal%Row(k1 + 1)%numel = nel
              allocate(Zal%Row(k1 + 1)%ind(nel), Zal%Row(k1 + 1)%Op(nel))
              Zal%Row(k1 + 1)%ind = Ham%Row(Ham%rbd)%ind(1:nel)

              ! First column
              call create(Tmp, Qid%nqs, 1)
              call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), Tmp, DD, Hcb, Idb, &
                                   locd, Dmap, Strmap, Qid)
              call extract_op_from_zaletel(Unonq, ZalU, 1, 0, locd)
              call extract_qop_from_op(Zal%Row(k1 + 1)%Op(1), Unonq, &
                                       Dmap, Strmap, locd, BB(k1), Qid)
              call skim(Zal%Row(k1 + 1)%Op(1), 10.0_rKind * numzero)
              call destroy(Tmp)

              ! Remaining columns
              cols: do k2 = 1, (nel - 1)
                  if(k1 == k2) then
                      call copy(Tmp, AA(k1))
                  else
                      call create(Tmp, Qid%nqs, 1)
                  end if

                  call get_zaletel_expm(ZalU, dt, Tmp, BB(k1), CC(k2), DD, &
                                       Hcb, Idb, locd, Dmap, Strmap, &
                                       Qid)

                  call extract_op_from_zaletel(Unonq, ZalU, 1, 1, locd)

                  ! Calculate commutators between terms
      ! To-do: commutator or anti-commutator
                  call contr(TmpA, BB(k1), CC(k2), [2], [1])
                  call contr(TmpB, CC(k2), BB(k1), [2], [1])
                  call gaxpy(TmpA, done, TmpB)

                  if(TmpA%nb == 0) then
                      call destroy(TmpA)
                      call copy(TmpA, Tmp)
                  end if

                  call extract_qop_from_op(Zal%Row(k1 + 1)%Op(k2 + 1), &
                                           Unonq, Dmap, Strmap, locd, &
                                           TmpA, Qid)
                  call skim(Zal%Row(k1 + 1)%Op(k2 + 1), 10.0_rKind * numzero)

                  call destroy(Tmp)
                  call destroy(TmpA)
                  call destroy(TmpB)
              end do cols
          end if
      end do rows

      call destroy(DD)

      do k1 = 1, chil - 1
          call destroy(AA(k1))
          call destroy(BB(k1))
      end do

      do k2 = 1, Ham%Row(Ham%rbd)%numel - 2
          call destroy(CC(k2))
      end do
      deallocate(AA, BB, CC)

      call destroy(Hcb)
      call destroy(Idb)
      call destroy(Qid)

      !DBGprint *, 'get_zaletel_bulk_site_qtensorc <<< <<<'

  end subroutine get_zaletel_bulk_site_qtensorc

ZaletelOps_f90.get_zaletel_mpo()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for all sites.

Arguments

ZalTYPE(mpoc), out

On exit this is the Zaletel-MPO for the time evolution of a time step dt.

HamTYPE(mpo), in

The MPO of the Hamiltonian is needed to build the Zaletel-MPO.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_mpo(Zal, Ham, dt, mpo_is_hermitian)
      type(mpoc), intent(out) :: Zal
      type(mpo), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! short cut to system size
      integer :: ll

      ! Number of sites and allocate MPOs, copy ti flat
      ll = Ham%ll
      Zal%ll = ll
      allocate(Zal%Ws(ll))
      Zal%ti = Ham%ti

      if(Ham%ti) then
          ! Translation invariant
          call get_zaletel_1st_site(Zal%Wl, Ham%Wl, dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Wr, Ham%Wr, dt, &
                                     mpo_is_hermitian)
          call get_zaletel_bulk_site(Zal%Wb, Ham%Wb, dt, &
                                     mpo_is_hermitian)

          call set_timpo_pointers(Zal)
      else
          ! Translational variant
          call get_zaletel_1st_site(Zal%Ws(1), Ham%Ws(1), dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Ws(ll), Ham%Ws(ll), dt, &
                                     mpo_is_hermitian)

          do ii = 2, (ll - 1)
              call get_zaletel_bulk_site(Zal%Ws(ii), Ham%Ws(ii), dt, &
                                         mpo_is_hermitian)
          end do
      end if

  end subroutine get_zaletel_mpo

ZaletelOps_f90.get_zaletel_mpoc()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for all sites.

Arguments

ZalTYPE(mpoc), out

On exit this is the Zaletel-MPO for the time evolution of a time step dt.

HamTYPE(mpoc), in

The MPO of the Hamiltonian is needed to build the Zaletel-MPO.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_mpoc(Zal, Ham, dt, mpo_is_hermitian)
      type(mpoc), intent(out) :: Zal
      type(mpoc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! short cut to system size
      integer :: ll

      ! Number of sites and allocate MPOs, copy ti flat
      ll = Ham%ll
      Zal%ll = ll
      allocate(Zal%Ws(ll))
      Zal%ti = Ham%ti

      if(Ham%ti) then
          ! Translation invariant
          call get_zaletel_1st_site(Zal%Wl, Ham%Wl, dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Wr, Ham%Wr, dt, &
                                     mpo_is_hermitian)
          call get_zaletel_bulk_site(Zal%Wb, Ham%Wb, dt, &
                                     mpo_is_hermitian)

          call set_timpo_pointers(Zal)
      else
          ! Translational variant
          call get_zaletel_1st_site(Zal%Ws(1), Ham%Ws(1), dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Ws(ll), Ham%Ws(ll), dt, &
                                     mpo_is_hermitian)

          do ii = 2, (ll - 1)
              call get_zaletel_bulk_site(Zal%Ws(ii), Ham%Ws(ii), dt, &
                                         mpo_is_hermitian)
          end do
      end if

  end subroutine get_zaletel_mpoc

ZaletelOps_f90.get_zaletel_qmpo()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for all sites.

Arguments

ZalTYPE(qmpoc), out

On exit this is the Zaletel-MPO for the time evolution of a time step dt.

HamTYPE(qmpo), in

The MPO of the Hamiltonian is needed to build the Zaletel-MPO.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_qmpo(Zal, Ham, dt, mpo_is_hermitian)
      type(qmpoc), intent(out) :: Zal
      type(qmpo), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! short cut to system size
      integer :: ll

      ! Number of sites and allocate MPOs, copy ti flat
      ll = Ham%ll
      Zal%ll = ll
      allocate(Zal%Ws(ll))
      Zal%ti = Ham%ti

      if(Ham%ti) then
          ! Translation invariant
          call get_zaletel_1st_site(Zal%Wl, Ham%Wl, dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Wr, Ham%Wr, dt, &
                                     mpo_is_hermitian)
          call get_zaletel_bulk_site(Zal%Wb, Ham%Wb, dt, &
                                     mpo_is_hermitian)

          call set_timpo_pointers(Zal)
      else
          ! Translational variant
          call get_zaletel_1st_site(Zal%Ws(1), Ham%Ws(1), dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Ws(ll), Ham%Ws(ll), dt, &
                                     mpo_is_hermitian)

          do ii = 2, (ll - 1)
              call get_zaletel_bulk_site(Zal%Ws(ii), Ham%Ws(ii), dt, &
                                         mpo_is_hermitian)
          end do
      end if

  end subroutine get_zaletel_qmpo

ZaletelOps_f90.get_zaletel_qmpoc()

fortran-subroutine - April 2016 (dj) Build the Zaletel MPO for all sites.

Arguments

ZalTYPE(qmpoc), out

On exit this is the Zaletel-MPO for the time evolution of a time step dt.

HamTYPE(qmpoc), in

The MPO of the Hamiltonian is needed to build the Zaletel-MPO.

dtCOMPLEX, in

Time step to be included in the exponential already multiplied with (- eye) and scalar according to order.

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine get_zaletel_qmpoc(Zal, Ham, dt, mpo_is_hermitian)
      type(qmpoc), intent(out) :: Zal
      type(qmpoc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! short cut to system size
      integer :: ll

      ! Number of sites and allocate MPOs, copy ti flat
      ll = Ham%ll
      Zal%ll = ll
      allocate(Zal%Ws(ll))
      Zal%ti = Ham%ti

      if(Ham%ti) then
          ! Translation invariant
          call get_zaletel_1st_site(Zal%Wl, Ham%Wl, dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Wr, Ham%Wr, dt, &
                                     mpo_is_hermitian)
          call get_zaletel_bulk_site(Zal%Wb, Ham%Wb, dt, &
                                     mpo_is_hermitian)

          call set_timpo_pointers(Zal)
      else
          ! Translational variant
          call get_zaletel_1st_site(Zal%Ws(1), Ham%Ws(1), dt, &
                                    mpo_is_hermitian)
          call get_zaletel_last_site(Zal%Ws(ll), Ham%Ws(ll), dt, &
                                     mpo_is_hermitian)

          do ii = 2, (ll - 1)
              call get_zaletel_bulk_site(Zal%Ws(ii), Ham%Ws(ii), dt, &
                                         mpo_is_hermitian)
          end do
      end if

  end subroutine get_zaletel_qmpoc

ZaletelOps_f90.zaletel_2nd_mpo()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 2nd order time evolution.

Arguments

ZalsTYPE(mpoc)(2), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(mpo), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye).

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_2nd_mpo(Zals, Ham, dt, mpo_is_hermitian)
      type(mpoc), dimension(2), intent(out) :: Zals
      type(mpo), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 2
          sc = dt * sc2(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_2nd_mpo

ZaletelOps_f90.zaletel_2nd_mpoc()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 2nd order time evolution.

Arguments

ZalsTYPE(mpoc)(2), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(mpoc), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye).

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_2nd_mpoc(Zals, Ham, dt, mpo_is_hermitian)
      type(mpoc), dimension(2), intent(out) :: Zals
      type(mpoc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 2
          sc = dt * sc2(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_2nd_mpoc

ZaletelOps_f90.zaletel_2nd_qmpo()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 2nd order time evolution.

Arguments

ZalsTYPE(qmpoc)(2), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(qmpo), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye).

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_2nd_qmpo(Zals, Ham, dt, mpo_is_hermitian)
      type(qmpoc), dimension(2), intent(out) :: Zals
      type(qmpo), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 2
          sc = dt * sc2(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_2nd_qmpo

ZaletelOps_f90.zaletel_2nd_qmpoc()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 2nd order time evolution.

Arguments

ZalsTYPE(qmpoc)(2), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(qmpoc), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye).

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_2nd_qmpoc(Zals, Ham, dt, mpo_is_hermitian)
      type(qmpoc), dimension(2), intent(out) :: Zals
      type(qmpoc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 2
          sc = dt * sc2(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_2nd_qmpoc

ZaletelOps_f90.zaletel_4th_mpo()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 4th order time evolution.

Arguments

ZalsTYPE(mpoc)(4), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(mpo), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye)

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_4th_mpo(Zals, Ham, dt, mpo_is_hermitian)
      type(mpoc), dimension(4), intent(out) :: Zals
      type(mpo), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 4
          sc = dt * sc4(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_4th_mpo

ZaletelOps_f90.zaletel_4th_mpoc()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 4th order time evolution.

Arguments

ZalsTYPE(mpoc)(4), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(mpoc), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye)

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_4th_mpoc(Zals, Ham, dt, mpo_is_hermitian)
      type(mpoc), dimension(4), intent(out) :: Zals
      type(mpoc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 4
          sc = dt * sc4(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_4th_mpoc

ZaletelOps_f90.zaletel_4th_qmpo()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 4th order time evolution.

Arguments

ZalsTYPE(qmpoc)(4), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(qmpo), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye)

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_4th_qmpo(Zals, Ham, dt, mpo_is_hermitian)
      type(qmpoc), dimension(4), intent(out) :: Zals
      type(qmpo), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 4
          sc = dt * sc4(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_4th_qmpo

ZaletelOps_f90.zaletel_4th_qmpoc()

fortran-subroutine - April 2016 (dj) Generate the Zaletel MPOs for 4th order time evolution.

Arguments

ZalsTYPE(qmpoc)(4), out

The array of MPOs is filled corresponding to the second order.

HamTYPE(qmpoc), in

Derive Zaletel MPO from this Hamiltonian-MPO.

dtREAL, in

Time step of the time evolution scaled with (- eye)

mpo_is_hermitianLOGICAL, in

Necessary to distighuish betweew exponentials used for local term.

Details

(template defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine zaletel_4th_qmpoc(Zals, Ham, dt, mpo_is_hermitian)
      type(qmpoc), dimension(4), intent(out) :: Zals
      type(qmpoc), intent(in) :: Ham
      complex(kind=rkind), intent(in) :: dt
      logical, intent(in) :: mpo_is_hermitian

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! Scalar to be included in the exponential
      complex(KIND=rKind) :: sc

      do ii = 1, 4
          sc = dt * sc4(ii)
          call get_zaletel(Zals(ii), Ham, sc, mpo_is_hermitian)
      end do

  end subroutine zaletel_4th_qmpoc

ZaletelOps_f90.LRK2_mpo()

fortran-subroutine - August 2017 (dj, updated) Propagate Psi with the LRK algorithm.

Arguments

PsiTYPE(mpsc), inout

On entry current state, on exit MPS after evolution step.

deltatREAL, in

Time step already scaled with (- eye)

HamTYPE(mpo), inout

Hamiltonian MPO for evolution of this step.

CpTYPE(ConvParam), in

Maximal bond dimension allowed for Psi.

convergedLOGICAL, out

On exit flag if step of algorithm converged according to convergence parameters.

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

initmethodCHARACTER, in

Not referenced, init method by default ‘F’ for LRK2

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK2_mpo(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(mpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(mpo), intent(inout) :: Ham
      type(ConvParam), intent(in) :: Cp
      logical, intent(in) :: mpo_is_hermitian
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(mpsc) :: Psi_fit

      ! MPOs of the propagator
      type(mpoc), dimension(2) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK2_mpo:'//&
      !        'Zaletel / LRK2 does not have PBC.', 99, 'ZaletelOps_include.f90:2393', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK2_mpo : set_hash failed.', &
      !               'ZaletelOps_include.f90:2400', errst=errst)) return

      call zaletel_2nd(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 2
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', renorm, errst=errst)
          !if(prop_error('LRK2_mpo : fit_hpsi failed.', &
          !              errst=errst)) return
          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK2_mpo

ZaletelOps_f90.LRK2_mpoc()

fortran-subroutine - August 2017 (dj, updated) Propagate Psi with the LRK algorithm.

Arguments

PsiTYPE(mpsc), inout

On entry current state, on exit MPS after evolution step.

deltatREAL, in

Time step already scaled with (- eye)

HamTYPE(mpoc), inout

Hamiltonian MPO for evolution of this step.

CpTYPE(ConvParam), in

Maximal bond dimension allowed for Psi.

convergedLOGICAL, out

On exit flag if step of algorithm converged according to convergence parameters.

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

initmethodCHARACTER, in

Not referenced, init method by default ‘F’ for LRK2

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK2_mpoc(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(mpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(mpoc), intent(inout) :: Ham
      type(ConvParam), intent(in) :: Cp
      logical, intent(in) :: mpo_is_hermitian
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(mpsc) :: Psi_fit

      ! MPOs of the propagator
      type(mpoc), dimension(2) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK2_mpoc:'//&
      !        'Zaletel / LRK2 does not have PBC.', 99, 'ZaletelOps_include.f90:2393', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK2_mpoc : set_hash failed.', &
      !               'ZaletelOps_include.f90:2400', errst=errst)) return

      call zaletel_2nd(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 2
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', renorm, errst=errst)
          !if(prop_error('LRK2_mpoc : fit_hpsi failed.', &
          !              errst=errst)) return
          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK2_mpoc

ZaletelOps_f90.LRK2_qmpo()

fortran-subroutine - August 2017 (dj, updated) Propagate Psi with the LRK algorithm.

Arguments

PsiTYPE(qmpsc), inout

On entry current state, on exit MPS after evolution step.

deltatREAL, in

Time step already scaled with (- eye)

HamTYPE(qmpo), inout

Hamiltonian MPO for evolution of this step.

CpTYPE(ConvParam), in

Maximal bond dimension allowed for Psi.

convergedLOGICAL, out

On exit flag if step of algorithm converged according to convergence parameters.

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

initmethodCHARACTER, in

Not referenced, init method by default ‘F’ for LRK2

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK2_qmpo(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(qmpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(qmpo), intent(inout) :: Ham
      type(ConvParam), intent(in) :: Cp
      logical, intent(in) :: mpo_is_hermitian
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(qmpsc) :: Psi_fit

      ! MPOs of the propagator
      type(qmpoc), dimension(2) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK2_qmpo:'//&
      !        'Zaletel / LRK2 does not have PBC.', 99, 'ZaletelOps_include.f90:2393', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK2_qmpo : set_hash failed.', &
      !               'ZaletelOps_include.f90:2400', errst=errst)) return

      call zaletel_2nd(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 2
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', renorm, errst=errst)
          !if(prop_error('LRK2_qmpo : fit_hpsi failed.', &
          !              errst=errst)) return
          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK2_qmpo

ZaletelOps_f90.LRK2_qmpoc()

fortran-subroutine - August 2017 (dj, updated) Propagate Psi with the LRK algorithm.

Arguments

PsiTYPE(qmpsc), inout

On entry current state, on exit MPS after evolution step.

deltatREAL, in

Time step already scaled with (- eye)

HamTYPE(qmpoc), inout

Hamiltonian MPO for evolution of this step.

CpTYPE(ConvParam), in

Maximal bond dimension allowed for Psi.

convergedLOGICAL, out

On exit flag if step of algorithm converged according to convergence parameters.

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

initmethodCHARACTER, in

Not referenced, init method by default ‘F’ for LRK2

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK2_qmpoc(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(qmpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(qmpoc), intent(inout) :: Ham
      type(ConvParam), intent(in) :: Cp
      logical, intent(in) :: mpo_is_hermitian
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(qmpsc) :: Psi_fit

      ! MPOs of the propagator
      type(qmpoc), dimension(2) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK2_qmpoc:'//&
      !        'Zaletel / LRK2 does not have PBC.', 99, 'ZaletelOps_include.f90:2393', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK2_qmpoc : set_hash failed.', &
      !               'ZaletelOps_include.f90:2400', errst=errst)) return

      call zaletel_2nd(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 2
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', renorm, errst=errst)
          !if(prop_error('LRK2_qmpoc : fit_hpsi failed.', &
          !              errst=errst)) return
          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK2_qmpoc

ZaletelOps_f90.LRK4_mpo()

fortran-subroutine - August 2017 (dj, updated) Propagate psi with the TEBD algorithm.

Arguments

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK4_mpo(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(mpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(mpo), intent(inout) :: Ham
      logical, intent(in) :: mpo_is_hermitian
      type(ConvParam), intent(in) :: Cp
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(mpsc) :: Psi_fit

      ! MPOs of the propagator
      type(mpoc), dimension(4) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK4_mpo:'//&
      !        'Zaletel / LRK4 does not have PBC.', 99, 'ZaletelOps_include.f90:2496', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK4_mpo : set_hash failed.', &
      !               'ZaletelOps_include.f90:2503', errst=errst)) return

      call zaletel_4th(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 4
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', 'M', errst=errst)
          !if(prop_error('LRK4_mpo : fit_hpsi failed.', &
          !              errst=errst)) return

          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK4_mpo

ZaletelOps_f90.LRK4_mpoc()

fortran-subroutine - August 2017 (dj, updated) Propagate psi with the TEBD algorithm.

Arguments

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK4_mpoc(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(mpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(mpoc), intent(inout) :: Ham
      logical, intent(in) :: mpo_is_hermitian
      type(ConvParam), intent(in) :: Cp
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(mpsc) :: Psi_fit

      ! MPOs of the propagator
      type(mpoc), dimension(4) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK4_mpoc:'//&
      !        'Zaletel / LRK4 does not have PBC.', 99, 'ZaletelOps_include.f90:2496', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK4_mpoc : set_hash failed.', &
      !               'ZaletelOps_include.f90:2503', errst=errst)) return

      call zaletel_4th(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 4
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', 'M', errst=errst)
          !if(prop_error('LRK4_mpoc : fit_hpsi failed.', &
          !              errst=errst)) return

          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK4_mpoc

ZaletelOps_f90.LRK4_qmpo()

fortran-subroutine - August 2017 (dj, updated) Propagate psi with the TEBD algorithm.

Arguments

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK4_qmpo(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(qmpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(qmpo), intent(inout) :: Ham
      logical, intent(in) :: mpo_is_hermitian
      type(ConvParam), intent(in) :: Cp
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(qmpsc) :: Psi_fit

      ! MPOs of the propagator
      type(qmpoc), dimension(4) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK4_qmpo:'//&
      !        'Zaletel / LRK4 does not have PBC.', 99, 'ZaletelOps_include.f90:2496', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK4_qmpo : set_hash failed.', &
      !               'ZaletelOps_include.f90:2503', errst=errst)) return

      call zaletel_4th(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 4
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', 'M', errst=errst)
          !if(prop_error('LRK4_qmpo : fit_hpsi failed.', &
          !              errst=errst)) return

          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK4_qmpo

ZaletelOps_f90.LRK4_qmpoc()

fortran-subroutine - August 2017 (dj, updated) Propagate psi with the TEBD algorithm.

Arguments

cerrREAL, inout

The cumulated error. The error done during this subroutine is added to the incoming value.

renormLOGICAL, in

Flag if state vector should be renormalized to 1 (true). ‘N’ : do not normalize (default); ‘M’ : normalize for MPS 1 / sqrt(norm);

pbcLOGICAL, in

If PBC are used in any rule set. There is a check before, but the debugging mode should ensure for now that there are no calls with such a check.

Details

Currently uses complete reorthogonalization for stability. (defined in ZaletelOps_include.f90)

Source Code

show / hide f90 code

  subroutine LRK4_qmpoc(Psi, deltat, Ham, Cp, mpo_is_hermitian, &
                           converged, cerr, renorm, pbc, errst)
      type(qmpsc), intent(inout) :: Psi
      complex(KIND=rKind), intent(in) :: deltat
      type(qmpoc), intent(inout) :: Ham
      logical, intent(in) :: mpo_is_hermitian
      type(ConvParam), intent(in) :: Cp
      logical, intent(out) :: converged
      real(KIND=rKind), intent(inout) :: cerr
      character, intent(in) :: renorm
      logical, intent(in) :: pbc
      integer, intent(out), optional :: errst

      ! Local variables
      ! ---------------

      ! for looping
      integer :: ii

      ! error in one step of fitting U | Psi>
      real(KIND=rKind) :: err

      ! Result of the fit
      type(qmpsc) :: Psi_fit

      ! MPOs of the propagator
      type(qmpoc), dimension(4) :: ZalUs

      !if(present(errst)) errst = 0

      !if(pbc) then
      !    errst = raise_error('LRK4_qmpoc:'//&
      !        'Zaletel / LRK4 does not have PBC.', 99, 'ZaletelOps_include.f90:2496', &
      !        errst=errst)
      !    return
      !end if

      call set_hash(Ham, [2], errst=errst)
      !if(prop_error('LRK4_qmpoc : set_hash failed.', &
      !               'ZaletelOps_include.f90:2503', errst=errst)) return

      call zaletel_4th(ZalUs, Ham, deltat, mpo_is_hermitian)

      do ii = 1, 4
          call Fit_Hpsi(Psi_fit, ZalUs(ii), Psi, Cp, err, converged, &
                        'F', 'E', 'M', errst=errst)
          !if(prop_error('LRK4_qmpoc : fit_hpsi failed.', &
          !              errst=errst)) return

          cerr = cerr + err

          call destroy(Psi)
          call copy(Psi, Psi_fit)
          call destroy(Psi_fit)

          call destroy(ZalUs(ii))
      end do

  end subroutine LRK4_qmpoc