bug in mrcc with tiny matrix for the dressing

plugins/MRCC/H_apply.irp.f

@@ -2,18 +2,20 @@ use bitmasks
 BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import *
 
-s = H_apply("mrcc_simple")
+s = H_apply("mrcc")
-s.data["parameters"] = ", delta_ij_sd_, Ndet_sd"
+s.data["parameters"] = ", delta_ij_, delta_ii_,Ndet_cas, Ndet_non_cas"
 s.data["declarations"] += """
-    integer, intent(in) :: Ndet_sd
+    integer, intent(in) :: Ndet_cas,Ndet_non_cas
-    double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
+    double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*)
+    double precision, intent(in) :: delta_ii_(Ndet_cas,*)
 """
-s.data["keys_work"] = "call mrcc_dress_simple(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)"
+s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_cas,Ndet_non_cas,i_generator,key_idx,keys_out,N_int,iproc)"
-s.data["params_post"] += ", delta_ij_sd_, Ndet_sd"
+s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas"
-s.data["params_main"] += "delta_ij_sd_, Ndet_sd"
+s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas"
 s.data["decls_main"] += """
-    integer, intent(in) :: Ndet_sd
+    integer, intent(in) :: Ndet_cas,Ndet_non_cas
-    double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
+    double precision, intent(in) :: delta_ij_(Ndet_cas,Ndet_non_cas,*)
+    double precision, intent(in) :: delta_ii_(Ndet_cas,*)
 """
 s.data["finalization"] = ""
 s.data["copy_buffer"] = ""
@@ -22,11 +24,5 @@ s.data["size_max"] = "3072"
 print s
 
 
-
-
-s.data["subroutine"] = "H_apply_mrcc"
-s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)"
-print s
-
 END_SHELL
 

plugins/MRCC/davidson.irp.f

@@ -35,11 +35,16 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit,i
   allocate(H_jj(sze))
   
   !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP  SHARED(sze,H_jj,dets_in,Nint,istate,delta_ij)           &
+      !$OMP  SHARED(sze,H_jj,N_det_cas,dets_in,Nint,istate,delta_ij,delta_ii,idx_cas)           &
       !$OMP  PRIVATE(i)
   !$OMP DO SCHEDULE(guided)
   do i=1,sze
-    H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + delta_ij(i,i,istate)
+    H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) 
+  enddo
+  !$OMP END DO 
+  !$OMP DO SCHEDULE(guided)
+  do i=1,N_det_cas
+    H_jj(idx_cas(i)) +=  delta_ii(i,istate)
   enddo
   !$OMP END DO 
   !$OMP END PARALLEL
@@ -370,16 +375,16 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
   integer, allocatable           :: idx(:)
   double precision               :: hij
   double precision, allocatable  :: vt(:)
-  integer                        :: i,j,k,l, jj
+  integer                        :: i,j,k,l, jj,ii
   integer                        :: i0, j0
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
   ASSERT (n>0)
-  PROVIDE ref_bitmask_energy delta_ij
+  PROVIDE ref_bitmask_energy delta_ij delta_ii
   integer, parameter             :: block_size = 157
   !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP PRIVATE(i,hij,j,k,idx,jj,vt)                             &
+      !$OMP PRIVATE(i,hij,j,k,idx,jj,ii,vt)                             &
-      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
+      !$OMP SHARED(n_det_cas,n_det_non_cas,idx_cas,idx_non_cas,n,H_jj,u_0,keys_tmp,Nint,v_0,istate,delta_ij)
   !$OMP DO SCHEDULE(static)
   do i=1,n  
     v_0(i) = H_jj(i) * u_0(i)
@@ -389,20 +394,30 @@ subroutine H_u_0_mrcc(v_0,u_0,H_jj,n,keys_tmp,Nint,istate)
   Vt = 0.d0
   !$OMP DO SCHEDULE(guided)
   do i=1,n
-!   idx(0) = i
+    idx(0) = i
-!   call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
+    call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx)
-!   do jj=1,idx(0)
+    do jj=1,idx(0)
-!     j = idx(jj)
+      j = idx(jj)
 !     if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then
-    do j = 1, i-1
         call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij)
-        hij = hij + delta_ij(j,i,istate)
+        hij = hij 
         vt (i) = vt (i) + hij*u_0(j)
         vt (j) = vt (j) + hij*u_0(i)
 !     endif
     enddo
   enddo
   !$OMP END DO
+
+  !$OMP DO SCHEDULE(guided)
+  do ii=1,n_det_cas
+    i = idx_cas(ii)
+    do jj = 1, n_det_non_cas
+        j = idx_non_cas(jj)
+        vt (i) = vt (i) + delta_ij(j,i,istate)*u_0(j)
+        vt (j) = vt (j) + delta_ij(j,i,istate)*u_0(i)
+    enddo
+  enddo
+  !$OMP END DO
   !$OMP CRITICAL
   do i=1,n
     v_0(i) = v_0(i) + vt(i)

plugins/MRCC/mrcc_dress.irp.f

@@ -12,13 +12,14 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_cas_lock, (psi_det_size) ]
 
 END_PROVIDER
 
-subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,iproc)
+subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_cas, Ndet_non_cas,i_generator,n_selected,det_buffer,Nint,iproc)
  use bitmasks
  implicit none
 
   integer, intent(in)            :: i_generator,n_selected, Nint, iproc
-  integer, intent(in) :: Ndet
+  integer, intent(in) :: Ndet_cas, Ndet_non_cas
-  double precision, intent(inout) :: delta_ij_(Ndet,Ndet,*)
+  double precision, intent(inout) :: delta_ij_(Ndet_cas,Ndet_non_cas,*)
+  double precision, intent(inout) :: delta_ii_(Ndet_cas,*)
 
   integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
   integer                        :: i,j,k,l
@@ -42,7 +43,7 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro
 
   call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
 
-  allocate (dIa_hla(N_states,Ndet))
+  allocate (dIa_hla(N_states,Ndet_non_cas))
 
   ! |I>
 
@@ -136,12 +137,11 @@ subroutine mrcc_dress(delta_ij_,Ndet,i_generator,n_selected,det_buffer,Nint,ipro
        do l_sd=1,idx_alpha(0)
          k_sd = idx_alpha(l_sd)
          do i_state=1,N_states
-           delta_ij_(idx_non_cas(k_sd),idx_cas(i_I),i_state) += dIa_hla(i_state,k_sd)
+           delta_ij_(i_I,k_sd,i_state) += dIa_hla(i_state,k_sd)
-           delta_ij_(idx_cas(i_I),idx_non_cas(k_sd),i_state) += dIa_hla(i_state,k_sd)
            if(dabs(psi_cas_coef(i_I,i_state)).ge.5.d-5)then
-            delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state)
+            delta_ii_(i_I,i_state) -= dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_cas_coef(k_sd,i_state)
            else
-            delta_ij_(idx_cas(i_I),idx_cas(i_I),i_state)  = 0.d0
+            delta_ii_(i_I,i_state)  = 0.d0
            endif
          enddo
        enddo

plugins/MRCC/mrcc_utils.irp.f

@@ -12,13 +12,11 @@ do i=1,N_det_non_cas
      size(psi_cas_coef,1), n_states, ihpsi)
    call i_h_j(psi_non_cas(1,1,i),psi_non_cas(1,1,i),N_int,hii)
    do k=1,N_states
-
    lambda_pert(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii)
-   lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k)
     if (dabs(ihpsi(k)).le.1.d-3) then
-      lambda_mrcc(k,i) = 1.d0 / (psi_cas_energy_diagonalized(k)-hii)
+      lambda_mrcc(k,i) = lambda_pert(k,i)
-      icount_manu = icount_manu+1
+    else
-      cycle
+      lambda_mrcc(k,i) = psi_non_cas_coef(i,k)/ihpsi(k)
     endif
    enddo
  enddo
@@ -28,66 +26,49 @@ END_PROVIDER
 
 
 
+!BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ]
+!implicit none
+!BEGIN_DOC
+!! Dressing matrix in SD basis
+!END_DOC
+!delta_ij_non_cas = 0.d0
+!call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas)
+!END_PROVIDER
 
-BEGIN_PROVIDER [ character*(32), dressing_type ]
+ BEGIN_PROVIDER [ double precision, delta_ij, (N_det_cas,N_det_non_cas,N_states) ]
- implicit none
+&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_cas,N_states) ]
- BEGIN_DOC
- ! [ Simple | MRCC ]
- END_DOC
- dressing_type = "MRCC"
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, delta_ij_non_cas, (N_det_non_cas, N_det_non_cas,N_states) ]
- implicit none
- BEGIN_DOC
- ! Dressing matrix in SD basis
- END_DOC
- delta_ij_non_cas = 0.d0
- call H_apply_mrcc_simple(delta_ij_non_cas,N_det_non_cas)
-END_PROVIDER
-
-BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
  implicit none
  BEGIN_DOC
  ! Dressing matrix in N_det basis
  END_DOC
  integer :: i,j,m
  delta_ij = 0.d0
- if (dressing_type == "MRCC") then
+ delta_ii = 0.d0
-   call H_apply_mrcc(delta_ij,N_det)
+ call H_apply_mrcc(delta_ij,delta_ii,N_det_cas,N_det_non_cas)
- else if (dressing_type == "Simple") then
-   do m=1,N_states
-    do j=1,N_det_non_cas
-     do i=1,N_det_non_cas
-       delta_ij(idx_non_cas(i),idx_non_cas(j),m) = delta_ij_non_cas(i,j,m) 
-     enddo
-    enddo
-   enddo
- endif
- do i = 1, N_det
-  do j = 1, N_det
-   do m = 1, N_states
-    if(isnan(delta_ij(j,i,m)))then
-    delta_ij(j,i,m) = 0.d0
-    endif
-   enddo
-  enddo
- enddo
-
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]
+BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
  implicit none
  BEGIN_DOC
  ! Dressed H with Delta_ij
  END_DOC
- integer                        :: i, j
+ integer                        :: i, j,istate,ii,jj
+ do istate = 1,N_states
    do j=1,N_det
      do i=1,N_det
-     h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1)
+       h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) 
+     enddo
+   enddo
+   do ii = 1, N_det_cas
+     i =idx_cas(ii)
+     h_matrix_dressed(i,i,istate) += delta_ii(ii,istate)
+    do jj = 1, N_det_non_cas
+     j =idx_cas(jj)
+     h_matrix_dressed(i,j,istate) += delta_ij(ii,jj,istate)
+     h_matrix_dressed(j,i,istate) += delta_ij(ii,jj,istate)
+    enddo
    enddo 
  enddo
-
 END_PROVIDER