# How to solve the non linear Coupled Cluster equations using Inexact Newton with DIIS?

**URL:** <https://discourse.julialang.org/t/how-to-solve-the-non-linear-coupled-cluster-equations-using-inexact-newton-with-diis/115291>\
**Category:** Chemistry\
**Tags:** nonlinear, tensors\
**Created:** [June 6, 2024, 6:31pm UTC](https://discourse.julialang.org/t/how-to-solve-the-non-linear-coupled-cluster-equations-using-inexact-newton-with-diis/115291 "2024-06-06T18:31:19Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [June 6, 2024, 6:31pm UTC](https://discourse.julialang.org/t/how-to-solve-the-non-linear-coupled-cluster-equations-using-inexact-newton-with-diis/115291/1 "2024-06-06T18:31:19Z")

</div>

I have been trying to solve the Coupled Cluster Amplitude Equations which are a set of non-linear equations using Julia. The aim is to find cluster amplitudes tensor `T2` such that the residual tensor (defined below) is zero (or has norm less than a certain tolerance value) :

```julia
function calculate_residual(g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,T2,fvv,foo,R2u,R2)
    @tensor begin
        R2u[a_1,a_2,i_1,i_2]= 0.5* g_vvoo[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] - g_vovo[a_2,i_3,a_3,i_2] * T2[a_1,a_3,i_1,i_3]- g_vovo[a_2,i_3,a_3,i_1] * T2[a_1,a_3,i_3,i_2]+ 2*g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_2,i_3]+ 0.5*g_oooo[i_3,i_4,i_1,i_2] * T2[a_1,a_2,i_3,i_4]+ fvv[a_2,a_3] * T2[a_1,a_3,i_1,i_2]+ 0.5*g_vvvv[a_1,a_2,a_3,a_4] * T2[a_3,a_4,i_1,i_2]- foo[i_3,i_2] * T2[a_1,a_2,i_1,i_3]+ 2*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_3,i_1,i_3]) * T2[a_2,a_4,i_2,i_4]- 2 *(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_3,i_1,i_3]) * T2[a_2,a_4,i_4,i_2]+ 0.5*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_3,i_3,i_1]) * T2[a_2,a_4,i_4,i_2]- 2*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_3,a_4,i_3,i_2]) * T2[a_1,a_2,i_1,i_4]- 1*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_4,i_1,i_3]) * T2[a_2,a_3,i_2,i_4]+ (g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_4,i_1,i_3]) * T2[a_2,a_3,i_4,i_2]+ (g_oovv[i_3,i_4,a_3,a_4] * T2[a_2,a_4,i_4,i_3]) * T2[a_1,a_3,i_1,i_2]+ (g_oovv[i_3,i_4,a_3,a_4] * T2[a_3,a_4,i_2,i_3]) * T2[a_1,a_2,i_1,i_4]- 2*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_2,a_3,i_4,i_3]) * T2[a_1,a_4,i_1,i_2]+ 0.5*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_4,i_3,i_2]) * T2[a_2,a_3,i_4,i_1]+ 0.5* (g_oovv[i_3,i_4,a_3,a_4] * T2[a_3,a_4,i_1,i_2]) * T2[a_1,a_2,i_3,i_4]
     end
     @tensor begin
        R2[a,b,i,j] := R2u[a,b,i,j] + R2u[b,a,j,i]
     end
     return R2
end #Only T2 (a rank 4 tensor) is variable, everything else
# are parameters

```

The problems I am facing can be categorized into two categories  
**One - Problems when trying to do it using a Package**  
I cannot find a package in Julia that can solve this set of equation using inexact newton (and also has the ability to implement DIIS) (where I can give my own Jacobian – as my function maps a rank 4 tensor to a rank 4 tensor, the Jacobian would be a rank 8 tensor!). I tried NLSolve.jl , it seems to incorrectly state that it’s method `anderson` is same as DIIS or Pulay mixing.

**Two - Problems in coding my own iterative solver with DIIS**  
After not getting an appropriate package, I tried to solve the equations by hand using the following initial solver :

> **Hand-written Iterative Solver**
>
> ```julia
> function ccd_by_hand(maxitr)
> R2u,R2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,K,fvv,foo,erhf = initialize_cc_variables();
> R_iter = zeros(typeof(R2), size(R2)) 
> T2_old = initialize_t2_only();
> normtol=1.0e-8
> # display(T2_old) # Using MP2 amplitudes
> e_new = calculate_ECCD(g_oovv,T2_old)
> e_old = copy(0.0)
> earr = []
> etol = 1.0e-6
> p_max = 6
> p_min = 2
> push!(earr,e_new[1]+erhf)
> println("Starting CCD Iterations with Convergence Criteria as ||R₂||<$normtol , ΔE<$etol and max iterations=$maxitr")
> println("-----------------------------------------------------------------------------------------------")
> for i in 1:maxitr
> @printf("Starting Iteration: %d Current Total Energy: %.8f\n", i, e_new[1]+erhf)
> println("Energy before updating = $(e_old[1]+erhf)")
> R2 = calculate_residual(g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,T2_old,fvv,foo,R2u,R2)
> conv,r2norm = check_convergence(R2,normtol,e_old,e_new,etol)
> if conv
> @printf("CCD Converged in %d iterations, no further updates required,current ||R2||=%.10f and ΔE = %.12f\n\n-------------------------\nE_Total = %.8f\n", i,r2norm,abs(e_old[1]-e_new[1]),e_new[1]+erhf)
> break
> end
> R_iter = calculate_R_iter(R_iter,R2,fvv,foo)
> @printf("CCD Not Converged in %d iterations, current ||R2|| = %.10f and ΔE = %.12f \n", i,r2norm,abs(e_old[1]-e_new[1]))
> push!(earr,e_new[1]+erhf)
> T2_old = update_T2(T2_old,R2,fvv,foo)
> e_old = e_new
> e_new = calculate_ECCD(g_oovv,T2_old)
> println("Energy after updating = $(e_new[1]+erhf)")
>         
> end
>     
> end
> 
> ```

The code does work and the solution converges, but as many might know, the Coupled Cluster Equations are very slow to converge and usually need some kind of convergence accelerators. I tried to implement DIIS using the residuals scaled by the elements of `foo` and `fvv` operators as the error vectors in the code below

> **Code with DIIS**
>
> ```julia
> using Printf,Plots
> include("ccd-helper.jl")
> function ccd_by_hand(maxitr)
> R2u,R2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,K,fvv,foo,erhf = initialize_cc_variables();
> R_iter = zeros(Float64,size(R2)) 
> T2_old = initialize_t2_only();
> normtol=1.0e-8
> # display(T2_old) # Using MP2 amplitudes
> e_new = calculate_ECCD(g_oovv,T2_old)
> e_old = copy(0.0)
> earr = []
> etol = 1.0e-6
> p_max = 6 # Maximum number of previous iterations to store for DIIS
> p_min = 2 # DIIS will start after this many iterations
> R_iter_storage = Array{typeof(R2)}(undef, 0)
> push!(earr,e_new[1]+erhf)
> println("Starting CCD Iterations with Convergence Criteria as ||R₂||<$normtol , ΔE<$etol and max iterations=$maxitr")
> println("-----------------------------------------------------------------------------------------------")
> for i in 1:maxitr
> @printf("\n-------------------------\nStarting Iteration: %d Current Total Energy: %.8f\n", i, e_new[1]+erhf)
> println("Energy before updating = $(e_old[1]+erhf)")
> R2 = calculate_residual(g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,T2_old,fvv,foo,R2u,R2)
> conv,r2norm = check_convergence(R2,normtol,e_old,e_new,etol)
> if conv
> @printf("CCD Converged in %d iterations, no further updates required,current ||R2||=%.10f and ΔE = %.12f\n\n-------------------------\nE_Total = %.8f\n", i,r2norm,abs(e_old[1]-e_new[1]),e_new[1]+erhf)
> break
> end
> R_iter = calculate_R_iter(R_iter,R2,fvv,foo)
> # display(R_iter)
> @printf("CCD Not Converged in %d iterations, current ||R2|| = %.10f and ΔE = %.12f \n", i,r2norm,abs(e_old[1]-e_new[1]))
> 
> if i >= p_min #DIIS 
> println("DIIS is being implemented in iteration $i")
> if i >= p_min + p_max - 1 # Adjusted condition to start popping correctly
> popfirst!(R_iter_storage)
> end
> push!(R_iter_storage, copy(R_iter))
> # R_iter = DIIS_by_hand(R_iter_storage, length(R_iter_storage))
> R_iter = DIIS_by_hand_constrained(R_iter_storage, length(R_iter_storage))
> println("NORM OF R_ITER = $(norm(R_iter))")
> elseif i == p_min - 1
> push!(R_iter_storage, copy(R_iter))
> end
> println(length(R_iter_storage))
> e_old = e_new
> e_new = calculate_ECCD(g_oovv,T2_old)
> T2_old = update_T2_using_R_iter(T2_old,R_iter)
> println("Energy after updating = $(e_new[1]+erhf)\n-------------------------\n")
> end
> return R_iter_storage
>     
> end
> maxitr = 200;
> R_iter_storage=ccd_by_hand(maxitr);
> 
> ```

However, on implementing this DIIS, the convergence becomes even slower now rather than accelerating and I cannot understand why. I tried to use the same implementation of DIIS as done in Fermi.jl . Below I have also given all the necessary subroutines I had defined in the above functions :

> **Definitions of user defined functions used in solver**
>
> ```julia
> function initialize_t2_only()
> h,g,hnuc,norb,nelec = main("water_dump.fcidump");
> T2 = zeros(size(g))
> f = zeros(size(h)) # fock matrix
> nocc = round(Int,nelec/2);
> nv = norb - nocc;
> erhf = calculate_ERHF(h,g,hnuc,nocc)
> ######################## Fock Matrix and K - Tensor Initialization ##############################
> for p in 1:norb , q in 1:norb
> s_pq=0
> for i in 1:nocc
> s_pq += 2*g[p,q,i,i]-g[p,i,i,q]
> end
> f[p,q] = h[p,q] + s_pq
> end
> K = zeros(nv,nv,nocc,nocc)
> # Note that K Tensor is indexed as K[a,b,i,j] where a,b are virtual orbitals and i,j are occupied orbitals
> # It has a lower dimensionality than the g tensor, as we dont have transition elements where 
> # the transition is happening from one virtual to another virtual or one occupied to another occupied.
> # In later versions of the code, we will use the g tensor directly to calculate the T2 amplitudes without
> # uselessly allocating an extra tensor.
> for a in 1:nv, b in 1:nv, i in 1:nocc, j in 1:nocc
> K[a,b,i,j] = g[nocc+a,i,nocc+b,j]
> end
> 
> ##################################################################################
> ################## Amplitude Initialization ######################################
> T2 = - K 
> for a in 1:nv, b in 1:nv, i in 1:nocc, j in 1:nocc
> T2[a,b,i,j] = T2[a,b,i,j]/(f[nocc+a,nocc+a] + f[nocc+b,nocc+b] - f[i,i] - f[j,j])
> # if abs(T2[a,b,i,j]) < 10e-8
> # T2[a,b,i,j] = 0.0
> # end
> end
> 
> ##################################################################################
> return T2
> end
> 
> function calculate_ECCD(g_oovv,T2)
> @tensoropt begin
> ECCSD[:] := 2*g_oovv[i_1,i_2,a_1,a_2] * T2[a_1,a_2,i_1,i_2] - g_oovv[i_1,i_2,a_1,a_2] * T2[a_1,a_2,i_2,i_1]
> end
> return ECCSD
> end
> 
> function calculate_residual(g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,T2,fvv,foo,R2u,R2)
> @tensor begin
> R2u[a_1,a_2,i_1,i_2]= 0.5* g_vvoo[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] - g_vovo[a_2,i_3,a_3,i_2] * T2[a_1,a_3,i_1,i_3]- g_vovo[a_2,i_3,a_3,i_1] * T2[a_1,a_3,i_3,i_2]+ 2*g_voov[a_1,i_3,i_1,a_3] * T2[a_2,a_3,i_2,i_3]+ 0.5*g_oooo[i_3,i_4,i_1,i_2] * T2[a_1,a_2,i_3,i_4]+ fvv[a_2,a_3] * T2[a_1,a_3,i_1,i_2]+ 0.5*g_vvvv[a_1,a_2,a_3,a_4] * T2[a_3,a_4,i_1,i_2]- foo[i_3,i_2] * T2[a_1,a_2,i_1,i_3]+ 2*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_3,i_1,i_3]) * T2[a_2,a_4,i_2,i_4]- 2 *(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_3,i_1,i_3]) * T2[a_2,a_4,i_4,i_2]+ 0.5*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_3,i_3,i_1]) * T2[a_2,a_4,i_4,i_2]- 2*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_3,a_4,i_3,i_2]) * T2[a_1,a_2,i_1,i_4]- 1*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_4,i_1,i_3]) * T2[a_2,a_3,i_2,i_4]+ (g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_4,i_1,i_3]) * T2[a_2,a_3,i_4,i_2]+ (g_oovv[i_3,i_4,a_3,a_4] * T2[a_2,a_4,i_4,i_3]) * T2[a_1,a_3,i_1,i_2]+ (g_oovv[i_3,i_4,a_3,a_4] * T2[a_3,a_4,i_2,i_3]) * T2[a_1,a_2,i_1,i_4]- 2*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_2,a_3,i_4,i_3]) * T2[a_1,a_4,i_1,i_2]+ 0.5*(g_oovv[i_3,i_4,a_3,a_4] * T2[a_1,a_4,i_3,i_2]) * T2[a_2,a_3,i_4,i_1]+ 0.5* (g_oovv[i_3,i_4,a_3,a_4] * T2[a_3,a_4,i_1,i_2]) * T2[a_1,a_2,i_3,i_4]
> end
> @tensor begin
> R2[a,b,i,j] := R2u[a,b,i,j] + R2u[b,a,j,i]
> end
> return R2
> end
> 
> function check_convergence(R2,normtol,e_old,e_new,etol)
> r2norm = sqrt(abs(dot(R2,R2)))
> if(r2norm < normtol && abs(e_new[1]-e_old[1])<etol)
> return true,r2norm
> else
> return false,r2norm
> end
>         
> end
> 
> function calculate_R_iter(R_iter,R2,fvv,foo)
> nv = size(fvv)[1]
> nocc = size(foo)[1]
> shiftp = 0.20
> for a in 1:nv, b in 1:nv, i in 1:nocc, j in 1:nocc
> R_iter[a,b,i,j] = R2[a,b,i,j]/(fvv[a,a] + fvv[b,b]-foo[i,i] - foo[j,j]+shiftp)
> end
> return R_iter
> end
> 
> function DIIS_by_hand_constrained(R_iter_storage,p)
> B = zeros(p+1,p+1)
> Bpp = dot(R_iter_storage[p],R_iter_storage[p])
> for i in 1:p
> for j in 1:i
> B[i,j] = dot(R_iter_storage[i],R_iter_storage[j])/Bpp
> B[j,i] = B[i,j]
> end
> end
> B[p+1, 1:p] .= -1
> B[1:p, p+1] .= -1
> B[p+1, p+1] = 0
> Id = zeros(p+1)
> Id[p+1]= -1
> B = Symmetric(B)
> display(B)
> # prob = LinearProblem(B,Id)
> # sol = solve(prob,LinearSolve.KrylovJL_CG())
> # C = sol.u
> # C = cholesky(B)\Id
> C = B\Id
> pop!(C) # As C[p+1] is the Lagrange multiplier
> s = zeros(size(R_iter_storage[1]))
> for i in 1:p
> s = s + C[i].*R_iter_storage[i]
> end
> println("The sum of coefficients of the DIIS are: ",sum(C))
> return s
> end
> 
> function update_T2_using_R_iter(T2,R_iter)
> nv = size(T2)[1]
> nocc = size(T2)[3]
> for a in 1:nv, b in 1:nv, i in 1:nocc, j in 1:nocc
> T2[a,b,i,j] = T2[a,b,i,j] - R_iter[a,b,i,j]
> end
> return T2
> end
> 
> Any suggestions that can help me progress in any of the two ways i'm exploring are really appreciated. 
> 
> ```

---

<div class="post-metadata">

**Author:** ![Uranium238](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/uranium238/32/209706_2.png) [@Uranium238](https://discourse.julialang.org/u/Uranium238)\
**Post date:** [June 16, 2024, 12:02pm UTC](https://discourse.julialang.org/t/how-to-solve-the-non-linear-coupled-cluster-equations-using-inexact-newton-with-diis/115291/2 "2024-06-16T12:02:09Z")

</div>

Due to the lack of any responses, I am using this reply to close this topic. There seems to be no packages available in Julia that can solve a non linear equation using inexact newton + DIIS. So the only way to implement it as of now is to code it by hand which at this point I have already done.