I am trying to write a code for an iterative solver of a set of non linear equations. Doing this needs computation of a lot of Tensor Contractions which I had been doing using TensorOperations.jl . My code is running fine but I also want to try to do the same using other Tensor frameworks. ITensor.jl is what I want to try out now to accomplish the same jobs. Before a small function to calculate a contraction and add it to a given workspace tensor looked like below
using ITensors,TensorOperations
using ITensors: permute!
include("ccd-helper.jl")
nv,nocc = 2,4;
function addres(R2u::Array{Float64,4},nv,nocc,T2::Array{Float64,4})
println("TensorOperations is being used")
g_voov = deserialize("g_voov.jlbin")
@tensor Trm2[a_1,a_2,i_1,i_2] := - g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_3,i_2]
@tensor R2u[a_1,a_2,i_1,i_2] += Trm2[a_1,a_2,i_1,i_2]
@tensor Trm5[a_1,a_2,i_1,i_2] := 2*g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_2,i_3]
@tensor R2u[a_1,a_2,i_1,i_2] += Trm5[a_1,a_2,i_1,i_2]
return R2u
end
Now, I cannot understand how to accomplish the same using ITensor.jl as elegantly as possible. The problem seems to stem from the fact that “Indices” are treated differently here and are just not dummy variables as they were when using TensorOperations.jl . I tried to implement the same in ITensors.jl as shown below
function addres(R2u::ITensor,nv,nocc,T2::ITensor)
println("ITensor is being used")
a_1 = Index(nv,"a_1")
a_2 = Index(nv,"a_2")
i_1 = Index(nocc,"i_1")
i_2 = Index(nocc,"i_2")
i_3 = Index(nocc,"i_3")
a_3 = Index(nv,"a_3")
g_voov = ITensor(deserialize("g_voov.jlbin"),a_1,i_3,i_1,a_3)
T2 = replaceinds(T2, inds(T2) => (a_2, a_3, i_3, i_2))
R2u = replaceinds(R2u, inds(R2u) => (a_1, a_2, i_1, i_2))
Trm2 = ITensor(a_1,a_2,i_1,i_2)
Trm5 = ITensor(a_1,a_2,i_1,i_2)
Trm2 = g_voov * T2
Trm2 = permute(Trm2,a_1,a_2,i_1,i_2)
R2u += Trm2
Trm5 = 2*g_voov * T2
Trm5 = permute(Trm5,a_1,a_2,i_1,i_2)
R2u += Trm5
return R2u
println("All Good")
end
But when trying to compare the outputs of the functions using the code below
initialize_cc_variables()
T2 = initialize_t2_only();
a = Index(nv,"a");
b = Index(nv,"b");
i = Index(nocc,"i");
j = Index(nocc,"j");
iT2 = ITensor(initialize_t2_only(),a,b,i,j);
R2u = zeros(Float64,nv,nv,nocc,nocc);
iR2u = ITensor(R2u,a,b,i,j);
R2u = addres(R2u,nv,nocc,T2)
iR2u = addres(iR2u,nv,nocc,iT2)
@show iR2u
display(R2u)
They seem to be giving different results. Is there any way to accomplish this code conversion without a significant amount of work ? (If you want to run the code exactly as shown above you can use the attached file which has been included at the start of the code, that file contains all the helper functions. The file t1.jl contains the current faulty code.)
ccd-helper.jl (26.6 KB)
t1.jl (1.7 KB)