Precomputed derivatives III

Author: lkdvos

src/algorithms/derivatives/mpo_derivatives.jl

@@ -119,87 +119,112 @@ struct PrecomputedDerivative{
     allocator::A
 end
 
-const PrecomputedACDerivative{T, S} = PrecomputedDerivative{T, S, 3, 2, 2, 1}
-const PrecomputedAC2Derivative{T, S} = PrecomputedDerivative{T, S, 3, 2, 3, 2}
+const PrecomputedACDerivative{T, S} = PrecomputedDerivative{T, S, 2, 1, 2, 1}
 
 function prepare_operator!!(
         H::MPO_AC_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPSTensor},
         x::MPSTensor,
         backend::AbstractBackend, allocator
     )
-    @plansor backend = backend allocator = allocator begin
-        GL_O[-1 -2 -3; -4 -5] := H.leftenv[-1 1; -4] * H.operators[1][1 -2; -5 -3]
-    end
-    leftenv = TensorMap(GL_O)
-    rightenv = TensorMap(H.rightenv)
-    BraidingStyle(sectortype(rightenv)) === NoBraiding() ||
-        (rightenv = braid(rightenv, ((2, 1), (3,)), (1, 2, 3)))
+    F_left = fuser(scalartype(x), codomain(x)...)
+    x′ = F_left * x
+
+    leftenv = left_precontract_derivative(H.leftenv, H.operators[1], F_left, backend, allocator)
+    rightenv = right_precontract_derivative(H.rightenv, backend, allocator)
+
+    return PrecomputedDerivative(leftenv, rightenv, backend, allocator), x′
+end
 
-    return PrecomputedDerivative(leftenv, rightenv, backend, allocator), x
+function unprepare_operator!!(y::MPSBondTensor, ::PrecomputedACDerivative, x::MPSTensor)
+    F_left = fuser(scalartype(x), codomain(x)...)
+    return F_left' * y
 end
 
 function prepare_operator!!(
         H::MPO_AC2_Hamiltonian{<:MPSTensor, <:MPOTensor, <:MPOTensor, <:MPSTensor},
         x::MPOTensor,
         backend::AbstractBackend, allocator
     )
-    @plansor backend = backend allocator = allocator begin
-        GL_O[-1 -2 -3; -4 -5] := H.leftenv[-1 1; -4] * H.operators[1][1 -2; -5 -3]
-        O_GR[-1 -2 -3; -4 -5] := H.operators[2][-3 -5; -2 1] * H.rightenv[-1 1; -4]
-    end
-    leftenv = TensorMap(GL_O)
-    rightenv = TensorMap(O_GR)
-    BraidingStyle(sectortype(rightenv)) === NoBraiding() ||
-        (rightenv = braid(rightenv, ((3, 1, 2), (4, 5)), (1, 2, 3, 4, 5)))
-    return PrecomputedDerivative(leftenv, rightenv, backend, allocator), x
+    F_left = fuser(scalartype(x), codomain(x)...)
+    F_right = fuser(scalartype(x), domain(x)...)
+    x′ = F_left * x * F_right'
+
+    leftenv = left_precontract_derivative(H.leftenv, H.operators[1], F_left, backend, allocator)
+    rightenv = right_precontract_derivative(H.rightenv, H.operators[2], F_right, backend, allocator)
+
+    return PrecomputedDerivative(leftenv, rightenv, backend, allocator), x′
+end
+
+function unprepare_operator!!(y::MPSBondTensor, ::PrecomputedDerivative, x::MPOTensor)
+    F_left = fuser(scalartype(x), codomain(x)...)
+    F_right = fuser(scalartype(x), domain(x)...)
+    return F_left' * y * F_right
 end
 
-function (H::PrecomputedDerivative)(x::MPSTensor)
+function (H::PrecomputedDerivative)(x::MPSBondTensor)
     bstyle = BraidingStyle(sectortype(x))
     return _precontracted_ac_derivative(bstyle, x, H.leftenv, H.rightenv, H.backend, H.allocator)
 end
 
-function _precontracted_ac_derivative(::NoBraiding, x, leftenv, rightenv, backend, allocator)
-    return @planar backend = backend allocator = allocator begin
-        y[-1 -2; -3] ≔ leftenv[-1 -2 4; 1 2] * x[1 2; 3] * rightenv[3 4; -3]
-    end
-end
 function _precontracted_ac_derivative(::Bosonic, x, leftenv, rightenv, backend, allocator)
     return @tensor backend = backend allocator = allocator begin
-        y[-1 -2; -3] ≔ leftenv[-1 -2 4; 1 2] * x[1 2; 3] * rightenv[4 3; -3]
+        y[-1; -2] ≔ leftenv[-1 3; 1] * x[1; 2] * rightenv[3 2; -2]
     end
 end
 function _precontracted_ac_derivative(::BraidingStyle, x, leftenv, rightenv, backend, allocator)
-    @planar backend = backend allocator = allocator begin
-        tmp[-1 -2 -3; -4] := leftenv[-1 -2 -3; 1 2] * x[1 2; -4]
+    return @planar backend = backend allocator = allocator begin
+        y[-1; -2] := leftenv[-1 2; 1] * x[1; 3] * τ'[3 2; 4 5] * rightenv[4 5; -2]
     end
-    tmp2 = braid(tmp, ((1, 2), (3, 4)), (1, 2, 4, 3))
+end
+function _precontracted_ac_derivative(::NoBraiding, x, leftenv, rightenv, backend, allocator)
     return @planar backend = backend allocator = allocator begin
-        y[-1 -2; -3] ≔ tmp2[-1 -2; 1 2] * rightenv[1 2; -3]
+        y[-1; -2] ≔ leftenv[-1 3; 1] * x[1; 2] * rightenv[2 3; -2]
     end
 end
 
-function (H::PrecomputedAC2Derivative)(x::MPOTensor)
-    bstyle = BraidingStyle(sectortype(x))
-    return _precontracted_ac2_derivative(bstyle, x, H.leftenv, H.rightenv, H.backend, H.allocator)
+left_precontract_derivative(arg, args...) = _left_precontract_derivative(BraidingStyle(sectortype(arg)), arg, args...)
+function _left_precontract_derivative(::BraidingStyle, leftenv, operator, F, backend, allocator)
+    @planar backend = backend allocator = allocator begin
+        GL_O[-1 -2 -3; -4 -5] := leftenv[-1 1; -4] * operator[1 -2; -5 -3]
+    end
+    return @planar backend = backend allocator = allocator begin
+        leftenv[-1 -2; -3] := F[-1; 3 4] * TensorMap(GL_O)[3 4 -2; 1 2] * F'[1 2; -3]
+    end
 end
 
-function _precontracted_ac2_derivative(::NoBraiding, x, leftenv, rightenv, backend, allocator)
+right_precontract_derivative(arg, args...) = _right_precontract_derivative(BraidingStyle(sectortype(arg)), arg, args...)
+_right_precontract_derivative(::NoBraiding, rightenv, backend, allocator) = TensorMap(rightenv)
+function _right_precontract_derivative(::Bosonic, rightenv, backend, allocator)
+    return @tensor backend = backend allocator = allocator begin
+        rightenv[-2 -1; -3] := TensorMap(rightenv)[-1 -2; -3]
+    end
+end
+function _right_precontract_derivative(::BraidingStyle, rightenv, backend, allocator)
     return @planar backend = backend allocator = allocator begin
-        y[-1 -2; -3 -4] ≔ leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * rightenv[3 4 5; -3 -4]
+        rightenv[-1 -2; -3] := τ[-1 -2; 1 2] * TensorMap(rightenv)[1 2; -3]
     end
 end
-function _precontracted_ac2_derivative(::Bosonic, x, leftenv, rightenv, backend, allocator)
+function _right_precontract_derivative(::NoBraiding, rightenv, operator, F, backend, allocator)
+    @planar backend = backend allocator = allocator begin
+        O_GR[-1 -2 -3; -4 -5] := operator[-3 -5; -2 1] * rightenv[-1 1; -4]
+    end
+    return @planar backend = backend allocator = allocator begin
+        rightenv[-1 -2; -3] := F[-1; 3 4] * TensorMap(O_GR)[3 4 -2; 1 2] * F'[1 2; -3]
+    end
+end
+function _right_precontract_derivative(::Bosonic, rightenv, operator, F, backend, allocator)
+    @tensor backend = backend allocator = allocator begin
+        O_GR[-1 -2 -3; -4 -5] := operator[-3 -5; -2 1] * rightenv[-1 1; -4]
+    end
     return @tensor backend = backend allocator = allocator begin
-        y[-1 -2; -3 -4] ≔ leftenv[-1 -2 5; 1 2] * x[1 2; 3 4] * rightenv[5 3 4; -3 -4]
+        rightenv[-2 -1; -3] := F[-1; 3 4] * TensorMap(O_GR)[3 4 -2; 1 2] * F'[1 2; -3]
     end
 end
-function _precontracted_ac2_derivative(::BraidingStyle, x, leftenv, rightenv, backend, allocator)
+function _right_precontract_derivative(::BraidingStyle, rightenv, operator, F, backend, allocator)
     @planar backend = backend allocator = allocator begin
-        tmp[-1 -2 -3; -4 -5] := leftenv[-1 -2 -3; 1 2] * x[1 2; -4 -5]
+        O_GR[-1 -2 -3; -4 -5] := operator[-3 -5; -2 1] * rightenv[-1 1; -4]
     end
-    tmp2 = braid(tmp, ((1, 2), (3, 4, 5)), (1, 2, 5, 3, 4))
     return @planar backend = backend allocator = allocator begin
-        y[-1 -2; -3 -4] := tmp2[-1 -2; 1 2 3] * rightenv[1 2 3; -3 -4]
+        rightenv[-1 -2; -3] := τ[-1 -2; 5 6] * F[5; 3 4] * TensorMap(O_GR)[3 4 6; 1 2] * F'[1 2; -3]
     end
 end