Add routines for CI optimization

src/PROPERTIES/properties_ci.irp.f

@@ -39,10 +39,10 @@ BEGIN_PROVIDER [ double precision, ci_h_psidet, (size_ci_h_psidet) ]
       T += det_alpha_grad_lapl(4,l,i)*det_beta_value (j)
     enddo
     do l=elec_beta_num+1,elec_num
-      T += det_beta_grad_lapl (4,l,j)*det_alpha_value(i) 
+      T += det_beta_grad_lapl (4,l,j)*det_alpha_value(i)
     enddo
     ci_h_psidet(k) = -0.5d0*T + E_pot * det_alpha_value(i)*det_beta_value (j)
-    ci_h_psidet(k) *= psidet_inv 
+    ci_h_psidet(k) *= psidet_inv
   enddo
 
   ci_h_psidet_min = min(ci_h_psidet_min,minval(ci_h_psidet))
@@ -54,10 +54,9 @@ END_PROVIDER
 BEGIN_PROVIDER [ double precision, ci_overlap_matrix, (size_ci_overlap_matrix) ]
   implicit none
   BEGIN_DOC
-  ! < det(i) | det(j) >
+  ! < det(i) |H| det(j) >
-  ! < det(i) | det(j) >
   !
-  ! Dimensions : det_num
+  ! Dimensions : det_num*det_num
   END_DOC
 
   integer                        :: i, j, k, l, m, n
@@ -70,7 +69,7 @@ BEGIN_PROVIDER [ double precision, ci_overlap_matrix, (size_ci_overlap_matrix) ]
     do l=1,det_num
       m = det_coef_matrix_rows(l)
       n = det_coef_matrix_columns(l)
-      ci_overlap_matrix(l) = det_alpha_value(m)*det_beta_value(n) * f
+      ci_overlap_matrix( det_num*(k-1) + l) = det_alpha_value(m)*det_beta_value(n) * f
     enddo
   enddo
 
@@ -79,3 +78,142 @@ BEGIN_PROVIDER [ double precision, ci_overlap_matrix, (size_ci_overlap_matrix) ]
   SOFT_TOUCH ci_overlap_matrix_min ci_overlap_matrix_max
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ double precision, ci_h_matrix, (size_ci_h_matrix) ]
+  implicit none
+  BEGIN_DOC
+  ! < det(i) |H| det(j) >
+  !
+  ! Dimensions : det_num*det_num
+  END_DOC
+
+  integer                        :: i, j, k, l, m, n, e
+  double precision :: f, g, h, T, V
+
+  do l=1,det_num
+      m = det_coef_matrix_rows(l)
+      n = det_coef_matrix_columns(l)
+      ! Lapl D 
+      g = 0.d0
+      do e=1,elec_alpha_num
+        g += det_alpha_grad_lapl(4,e,m) * det_beta_value (n)
+      enddo
+      do e=elec_beta_num+1,elec_num
+        g += det_alpha_value(m) * det_beta_grad_lapl(4,e,n)
+      enddo
+      T = g
+      ! D (Lapl J)/J
+      g = 0.d0
+      do e=1,elec_num
+        g += jast_lapl_jast_inv(e)
+      enddo
+      T += det_alpha_value(m) * det_beta_value(n) * g
+      ! 2 (grad D).(Grad J)/J
+      g = 0.d0
+      do e=1,elec_alpha_num
+        g += &
+          det_alpha_grad_lapl(1,e,m) * jast_grad_jast_inv_x(e) + &
+          det_alpha_grad_lapl(2,e,m) * jast_grad_jast_inv_y(e) + &
+          det_alpha_grad_lapl(3,e,m) * jast_grad_jast_inv_z(e)
+      enddo
+      h = 0.d0
+      do e=elec_alpha_num+1,elec_num
+        h += &
+          det_beta_grad_lapl(1,e,n) * jast_grad_jast_inv_x(e) + &
+          det_beta_grad_lapl(2,e,n) * jast_grad_jast_inv_y(e) + &
+          det_beta_grad_lapl(3,e,n) * jast_grad_jast_inv_z(e)
+      enddo
+      T += 2.d0*( g * det_beta_value(n) + h * det_alpha_value(m) )
+      g = det_alpha_value(m)*det_beta_value(n)
+      V = E_pot* g
+      do e=1,elec_alpha_num
+        V -= pseudo_non_local(e)* g
+        V += det_alpha_pseudo(e,m) * det_beta_value(n)
+      enddo
+      do e=elec_alpha_num+1,elec_num
+        V -= pseudo_non_local(e)* g
+        V += det_alpha_value(m) * det_beta_pseudo(e,n)
+      enddo
+      f = -0.5d0*T + V
+      f *= psidet_inv * psidet_inv
+      do k=1,det_num
+        i = det_coef_matrix_rows(k)
+        j = det_coef_matrix_columns(k)
+        ci_h_matrix( det_num*(l-1) + k) = f * &
+           det_alpha_value(i)*det_beta_value (j)
+      enddo
+  enddo
+
+  ci_h_matrix_min = min(ci_h_matrix_min,minval(ci_h_matrix))
+  ci_h_matrix_max = max(ci_h_matrix_max,maxval(ci_h_matrix))
+  SOFT_TOUCH ci_h_matrix_min ci_h_matrix_max
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ci_h_matrix_diag, (size_ci_h_matrix_diag) ]
+  implicit none
+  BEGIN_DOC
+  ! < det(i) |H| det(j) >
+  !
+  ! Dimensions : det_num
+  END_DOC
+
+  integer                        :: i, j, k, l, m, n, e
+  double precision :: f, g, h, T, V
+
+  do l=1,det_num
+      m = det_coef_matrix_rows(l)
+      n = det_coef_matrix_columns(l)
+      ! Lapl D 
+      g = 0.d0
+      do e=1,elec_alpha_num
+        g += det_alpha_grad_lapl(4,e,m) * det_beta_value (n)
+      enddo
+      do e=elec_beta_num+1,elec_num
+        g += det_alpha_value(m) * det_beta_grad_lapl(4,e,n)
+      enddo
+      T = g
+      ! D (Lapl J)/J
+      g = 0.d0
+      do e=1,elec_num
+        g += jast_lapl_jast_inv(e)
+      enddo
+      T += det_alpha_value(m) * det_beta_value(n) * g
+      ! 2 (grad D).(Grad J)/J
+      g = 0.d0
+      do e=1,elec_alpha_num
+        g += &
+          det_alpha_grad_lapl(1,e,m) * jast_grad_jast_inv_x(e) + &
+          det_alpha_grad_lapl(2,e,m) * jast_grad_jast_inv_y(e) + &
+          det_alpha_grad_lapl(3,e,m) * jast_grad_jast_inv_z(e)
+      enddo
+      h = 0.d0
+      do e=elec_alpha_num+1,elec_num
+        h += &
+          det_beta_grad_lapl(1,e,n) * jast_grad_jast_inv_x(e) + &
+          det_beta_grad_lapl(2,e,n) * jast_grad_jast_inv_y(e) + &
+          det_beta_grad_lapl(3,e,n) * jast_grad_jast_inv_z(e)
+      enddo
+      T += 2.d0*( g * det_beta_value(n) + h * det_alpha_value(m) )
+      g = det_alpha_value(m)*det_beta_value(n)
+      V = E_pot* g
+      do e=1,elec_alpha_num
+        V -= pseudo_non_local(e)* g
+        V += det_alpha_pseudo(e,m) * det_beta_value(n)
+      enddo
+      do e=elec_alpha_num+1,elec_num
+        V -= pseudo_non_local(e)* g
+        V += det_alpha_value(m) * det_beta_pseudo(e,n)
+      enddo
+      f = -0.5d0*T + V
+      f *= psidet_inv * psidet_inv
+      ci_h_matrix_diag(l) = f * &
+           det_alpha_value(m)*det_beta_value (n)
+  enddo
+
+  ci_h_matrix_diag_min = min(ci_h_matrix_diag_min,minval(ci_h_matrix_diag))
+  ci_h_matrix_diag_max = max(ci_h_matrix_diag_max,maxval(ci_h_matrix_diag))
+  SOFT_TOUCH ci_h_matrix_diag_min ci_h_matrix_diag_max
+END_PROVIDER
+
+