Improved MRCC

plugins/MRCC_CASSD/EZFIO.cfg

@@ -7,11 +7,11 @@ interface: ezfio
 type: Threshold
 doc: Threshold on the convergence of the MRCC energy
 interface: ezfio,provider,ocaml
-default: 1.e-7
+default: 1.e-5
 
 [n_it_mrcc_max]
 type: Strictly_positive_int
 doc: Maximum number of MRCC iterations
 interface: ezfio,provider,ocaml
-default: 20
+default: 10
 

plugins/MRCC_CASSD/mrcc_cassd.irp.f

@@ -1,13 +1,47 @@
 program mrcc
   implicit none
-  if (.not.read_wf) then
-    print *,  'read_wf has to be true.'
-    stop 1
-  endif
+  read_wf = .True.
+  SOFT_TOUCH read_wf
   call print_cas_coefs
-  call run_mrcc
+  call set_generators_bitmasks_as_holes_and_particles
+  call run
 end
 
+subroutine run
+  implicit none
+  
+  integer :: i
+
+  double precision :: E_new, E_old, delta_e
+  integer :: iteration
+  double precision :: E_past(4), lambda
+  E_new = 0.d0
+  delta_E = 1.d0
+  iteration = 0
+  call diagonalize_ci_dressed(lambda)
+  lambda = 1.d0
+  do while (delta_E > thresh_mrcc)
+    iteration += 1
+    print *,  '===========================' 
+    print *,  'MRCC Iteration', iteration
+    print *,  '===========================' 
+    print *,  ''
+    E_old = sum(ci_energy_dressed)
+    call write_double(6,ci_energy_dressed(1),"MRCC energy")
+    call diagonalize_ci_dressed(lambda)
+    E_new = sum(ci_energy_dressed)
+    delta_E = dabs(E_new - E_old)
+    call save_wavefunction
+    call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
+    if (iteration > n_it_mrcc_max) then
+      exit
+    endif
+  enddo
+  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
+
+end
+
+
 subroutine print_cas_coefs
   implicit none
 
@@ -18,7 +52,7 @@ subroutine print_cas_coefs
     print *,  psi_cas_coef(i,:)
     call debug_det(psi_cas(1,1,i),N_int)
   enddo
-
   call write_double(6,ci_energy(1),"Initial CI energy")
+
 end
 

plugins/MRCC_Utils/mrcc_general.irp.f

@@ -11,12 +11,13 @@ subroutine mrcc_iterations
 
   double precision :: E_new, E_old, delta_e
   integer :: iteration,i_oscillations
-  double precision :: E_past(4)
+  double precision :: E_past(4), lambda
   E_new = 0.d0
   delta_E = 1.d0
   iteration = 0
   j = 1
   i_oscillations = 0
+  lambda = 1.d0
   do while (delta_E > 1.d-7)
     iteration += 1
     print *,  '===========================' 
@@ -25,10 +26,15 @@ subroutine mrcc_iterations
     print *,  ''
     E_old = sum(ci_energy_dressed)
     call write_double(6,ci_energy_dressed(1),"MRCC energy")
-    call diagonalize_ci_dressed
+    call diagonalize_ci_dressed(lambda)
     E_new = sum(ci_energy_dressed)
     delta_E = dabs(E_new - E_old)
-
+!    if (E_new > E_old) then
+!       lambda = lambda * 0.7d0
+!    else
+!       lambda = min(1.d0, lambda * 1.1d0)
+!    endif
+!    print *,  'energy lambda ', lambda
     E_past(j) = E_new
     j +=1
     call save_wavefunction

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -4,29 +4,44 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
  ! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
  END_DOC
  integer :: i,k
- double precision               :: ihpsi(N_states),ihpsi_current(N_states)
- integer :: i_pert_count
- 
- i_pert_count = 0
- lambda_mrcc = 0.d0
- 
-  do i=1,N_det_non_ref
-    call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
-    do k=1,N_states
-      if (ihpsi_current(k) == 0.d0) then
-        ihpsi_current(k) = 1.d-32
+  double precision               :: ihpsi_current(N_states)
+  integer :: i_pert_count
+  double precision :: hii, lambda_pert
+  
+  i_pert_count = 0
+  lambda_mrcc = 0.d0
+  
+   do i=1,N_det_non_ref
+     call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref, &
+                  size(psi_ref_coef,1), N_states,ihpsi_current)
+     call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
+     do k=1,N_states
+       if (ihpsi_current(k) == 0.d0) then
+         ihpsi_current(k) = 1.d-32
+       endif
+       lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k) 
+      if ( dabs(psi_non_ref_coef(i,k)*ihpsi_current(k)) < 1.d-8 ) then
+           i_pert_count += 1
+             lambda_mrcc(k,i) = 0.d0
+!          lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
+!          if((ihpsi_current(k) * lambda_pert) < 0.5d0 * psi_non_ref_coef_restart(i,k) ) then
+!             lambda_mrcc(k,i) = 0.d0
+!          endif
       endif
-      if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then
-        i_pert_count +=1
-      else
-        lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
-      endif
-    enddo
- enddo
+       double precision, parameter :: x = 2.d0
+       if (lambda_mrcc(k,i) > x) then
+         lambda_mrcc(k,i) = x
+       else if (lambda_mrcc(k,i) < -x) then
+         lambda_mrcc(k,i) = -x
+       endif
+     enddo
+  enddo
+ 
+  print*,'N_det_non_ref = ',N_det_non_ref
+  print*,'Number of ignored determinants = ',i_pert_count
+  print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
+  print*,'lambda min/max = ',maxval(dabs(lambda_mrcc)), minval(dabs(lambda_mrcc))
 
- print*,'N_det_non_ref = ',N_det_non_ref
- print*,'Number of ignored determinants = ',i_pert_count
- print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
 END_PROVIDER
 
 
@@ -174,17 +189,19 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
 
 END_PROVIDER
 
-subroutine diagonalize_CI_dressed
+subroutine diagonalize_CI_dressed(lambda)
   implicit none
   BEGIN_DOC
 !  Replace the coefficients of the CI states by the coefficients of the 
 !  eigenstates of the CI matrix
   END_DOC
+  double precision, intent(in) :: lambda
   integer :: i,j
   do j=1,N_states_diag
     do i=1,N_det
-      psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
+      psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j) 
     enddo
+    call normalize(psi_coef(1,j), N_det)
   enddo
   SOFT_TOUCH psi_coef