From a041df6641258b72899dc22729c54b14af208e04 Mon Sep 17 00:00:00 2001
From: arnab82 <arnabbachhar39@gmail.com>
Date: Tue, 24 Oct 2023 12:28:12 -0400
Subject: [PATCH] trust region for step added in newton

---
 src/finite_difference.jl |  3 ++-
 src/incore_cmf.jl        | 15 +++++++++------
 2 files changed, 11 insertions(+), 7 deletions(-)

diff --git a/src/finite_difference.jl b/src/finite_difference.jl
index 10b6f6c..233a56f 100644
--- a/src/finite_difference.jl
+++ b/src/finite_difference.jl
@@ -747,7 +747,7 @@ end
                                     stepsize = 1e-6)
 Compute orbital hessian with finite difference for cMF using orbital gradient 
 """
-function orbital_hessian_numerical(ints, clusters, kappa, fspace, ansatze,d::RDM1; verbose = 0,step_size = 1e-5,zero_intra_rots = true,maxiter_ci = 100, maxiter_d1 = 100, tol_d1 = 1e-6, tol_ci = 1e-8,sequential = false)
+function orbital_hessian_numerical(ints, clusters, kappa, fspace, ansatze,d::RDM1; verbose = 0,step_size = 5e-5,zero_intra_rots = true,maxiter_ci = 100, maxiter_d1 = 100, tol_d1 = 1e-6, tol_ci = 1e-8,sequential = false)
     n = length(kappa)
     #error("here.....")
     hessian = zeros(n, n)
@@ -782,6 +782,7 @@ function orbital_hessian_numerical(ints, clusters, kappa, fspace, ansatze,d::RDM
     norb = n_orb(ints)
     #central difference 
     for i in 1:n
+        println(i)
         x_plus_i = deepcopy(kappa)
         x_minus_i = deepcopy(kappa)
         x_plus_i[i] += step_size
diff --git a/src/incore_cmf.jl b/src/incore_cmf.jl
index cf83415..4914b6f 100644
--- a/src/incore_cmf.jl
+++ b/src/incore_cmf.jl
@@ -1635,7 +1635,8 @@ function cmf_oo_newton( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fsp
                     verbose         = 0, 
                     use_pyscf=true,
                     zero_intra_rots =true,
-                    sequential      = false
+                    sequential      = false,
+                    step_trust_region=false
     ) where T
     #={{{=#
     println(" Solve OO-CMF with newton")
@@ -1703,6 +1704,11 @@ function cmf_oo_newton( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fsp
         else
             step_i=-(pinv(h_i)*(g_i))#;atol=1e-8  
         end
+        if step_trust_region==true
+            if norm(step_i)> 0.95
+                step_i=step_i*0.95/norm(step_i)
+            end
+        end
         e = ei
         U = U*Ui    
         converged = norm(g_i) < tol_oo 
@@ -1711,7 +1717,7 @@ function cmf_oo_newton( ints_in::InCoreInts{T}, clusters::Vector{MOCluster}, fsp
             # display(packed_hessian)
             break
         else
-            @printf(" Step: %4i E: %16.12f G: %12.2e  \n", i, ei, norm(g_i))
+            @printf(" Step: %4i E: %16.12f G: %12.2e step_size: %12.2e \n", i, ei, norm(g_i), norm(step_i))
         end
     end
     return  e,U,d1  
@@ -1730,7 +1736,7 @@ function convert_pairs(original_list, rearranged_pairs)
 
     return reverted_pairs
 end
-function projection_vector(ansatze::Vector{<:Ansatz}, clusters, norb;gradient=false)
+function projection_vector(ansatze::Vector{<:Ansatz}, clusters, norb)
     n_dim = norb * (norb - 1) ÷ 2 #={{{=#
     
     tmp_mat = Matrix(1I, n_dim, n_dim)
@@ -1745,9 +1751,6 @@ function projection_vector(ansatze::Vector{<:Ansatz}, clusters, norb;gradient=fa
         append!(invar, tmp_global)
         # println(invar)
     end
-    if gradient == true
-        return invar
-    end
 
     fci = ActiveSpaceSolvers.FCIAnsatz(norb, 0, 0) #dummie FCI anstaz to generate all pairs
     full_list = ActiveSpaceSolvers.invariant_orbital_rotations(fci)