Basically, I am trying to solve the Coupled Cluster Equations using Julia and the package NLsolve which can be used to solve such nonlinear equations. My initial strategy was to define a function as follows that accepts only one matrix as a parameter

```
function f_old(T2)
R2u,R2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,K,fvv,foo,erhf = initialize_cc_variables();
R2 = calculate_residual(g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,T2,fvv,foo,R2u,R2);
R2
end
T2guess = initialize_t2_only()
sol=nlsolve(f_old, T2guess)
finalamps = sol.zero
R2u,R2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,K,fvv,foo,erhf = initialize_cc_variables()
ECCSD = calculate_ECCSD(g_oovv,finalamps)
ECCSD[1]+erhf
```

which gives the output

```
* Inf-norm of residuals: 0.000000
* Iterations: 15
* Convergence: true
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: true
* Function Calls (f): 16
* Jacobian Calls (df/dx): 16
```

The variables `g_voooo`

etc are very difficult to compute, and hence itâ€™s a bad idea to initialize them every time the function `f`

is called (they remain constant throughout every iteration). So I tried to follow this post by @ufechner7 , and as per the post, I modified my function and succeeding statements as follows

```
function f!(F,T2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,fvv,foo,R2u,R2)
R2 = calculate_residual(g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,T2,fvv,foo,R2u,R2);
F = R2
F
end
R2u,R2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,K,fvv,foo,erhf = initialize_cc_variables();
T2guess = initialize_t2_only()
sol = nlsolve((F,T2)->f!(F,T2,g_vvoo,g_voov,g_vovo,g_oovv,g_oooo,g_vvvv,fvv,foo,R2u,R2), T2guess)
finalamps = sol.zero
ECCSD = calculate_ECCSD(g_oovv,finalamps)
ECCSD[1]+erhf
```

But the later gives me the output

```
* Inf-norm of residuals: 0.067053
* Iterations: 1000
* Convergence: false
* |x - x'| < 0.0e+00: false
* |f(x)| < 1.0e-08: false
* Function Calls (f): 1
* Jacobian Calls (df/dx): 1
```

which seems to be incorrect because it seems like it converges in just one iteration? Also the the quantity `ECCSD[1]+erhf`

in the later case has a less negative value (a better solution should produce a more negative value of this quantity since this is basically the electronic energy). Any help is greatly appreciated !