Reduced error in MRCC

plugins/MRCC_Utils/davidson.irp.f

@@ -35,21 +35,20 @@ subroutine davidson_diag_mrcc(dets_in,u_in,energies,dim_in,sze,N_st,N_st_diag,Ni
   PROVIDE mo_bielec_integrals_in_map
   allocate(H_jj(sze))
   
+  H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint) 
   !$OMP PARALLEL DEFAULT(NONE)                                       &
       !$OMP  SHARED(sze,H_jj,N_det_ref,dets_in,Nint,istate,delta_ii,idx_ref)           &
       !$OMP  PRIVATE(i)
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO 
-  do i=1,sze
+  do i=2,sze
     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_ref
-    H_jj(idx_ref(i)) +=  delta_ii(istate,i)
-  enddo
-  !$OMP END DO 
   !$OMP END PARALLEL
 
+  do i=1,N_det_ref
+    H_jj(idx_ref(i)) +=  delta_ii(istate,i)
+  enddo
   call davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate)
   deallocate (H_jj)
 end
@@ -224,17 +223,6 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s
           W(i,k,iter+1) = 0.d0
         enddo
       enddo
-!      do k=1,N_st_diag
-!         do iter2=1,iter
-!          do l=1,N_st_diag
-!            do i=1,sze
-!              U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1)
-!              W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1)
-!            enddo
-!          enddo
-!        enddo
-!      enddo
-!
 !
       call dgemm('N','N', sze, N_st_diag, N_st_diag*iter,            &
           1.d0, U, size(U,1), y, size(y,1)*size(y,2), 0.d0, U(1,1,iter+1), size(U,1))
@@ -276,27 +264,11 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,N_s
       
       do k=1,N_st_diag
 
-!        do iter2=1,iter
-!          do l=1,N_st_diag
-!            c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze)
-!            do i=1,sze
-!              U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter2)
-!            enddo
-!          enddo
-!        enddo
-!
         call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1),  &
               U(1,k,iter+1),1,0.d0,c,1)
         call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
               c,1,1.d0,U(1,k,iter+1),1)
-!
+
-!        do l=1,k-1
-!          c(1) = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze)
-!          do i=1,sze
-!            U(i,k,iter+1) = U(i,k,iter+1) - c(1) * U(i,l,iter+1)
-!          enddo
-!        enddo
-!
         call dgemv('T',sze,k-1,1.d0,U(1,1,iter+1),size(U,1),   &
             U(1,k,iter+1),1,0.d0,c,1)
         call dgemv('N',sze,k-1,-1.d0,U(1,1,iter+1),size(U,1),        &
@@ -429,7 +401,7 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
   allocate(vt(sze_8,N_st))
   Vt = 0.d0
   
-  !$OMP DO SCHEDULE(dynamic)
+  !$OMP DO SCHEDULE(static,1)
   do sh=1,shortcut(0,1)
     do sh2=sh,shortcut(0,1)
       exa = 0
@@ -468,9 +440,9 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
       enddo
     enddo
   enddo
-  !$OMP END DO NOWAIT
+  !$OMP END DO 
   
-  !$OMP DO SCHEDULE(dynamic)
+  !$OMP DO SCHEDULE(static,1)
   do sh=1,shortcut(0,2)
     do i=shortcut(sh,2),shortcut(sh+1,2)-1
       org_i = sort_idx(i,2)
@@ -490,7 +462,7 @@ subroutine H_u_0_mrcc_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,istate_in,N_st,sze_8)
       end do
     end do
   enddo
-  !$OMP END DO NOWAIT
+  !$OMP END DO
   
   !$OMP DO 
   do ii=1,n_det_ref
@@ -562,22 +534,23 @@ subroutine davidson_diag_mrcc_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_dia
   PROVIDE mo_bielec_integrals_in_map 
   allocate(H_jj(sze), S2_jj(sze))
   
+  H_jj(1) = diag_h_mat_elem(dets_in(1,1,1),Nint)
+  call get_s2(dets_in(1,1,1),dets_in(1,1,1),Nint,S2_jj(1))
   !$OMP PARALLEL DEFAULT(NONE)                                       &
       !$OMP  SHARED(sze,H_jj,S2_jj, dets_in,Nint,N_det_ref,delta_ii, &
       !$OMP  idx_ref, istate)  &
       !$OMP  PRIVATE(i)
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO 
-  do i=1,sze
+  do i=2,sze
     H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint)
     call get_s2(dets_in(1,1,i),dets_in(1,1,i),Nint,S2_jj(i))
   enddo
   !$OMP END DO 
-  !$OMP DO SCHEDULE(guided)
+  !$OMP END PARALLEL
+
   do i=1,N_det_ref
     H_jj(idx_ref(i)) +=  delta_ii(istate,i)
   enddo
-  !$OMP END DO 
-  !$OMP END PARALLEL
 
   call davidson_diag_hjj_sjj_mrcc(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit,istate)
   deallocate (H_jj,S2_jj)
@@ -1051,7 +1024,7 @@ subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_i
   Vt = 0.d0
   St = 0.d0
 
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(static,1)
   do sh=1,shortcut(0,1)
     do sh2=sh,shortcut(0,1)
       exa = 0
@@ -1094,7 +1067,7 @@ subroutine H_S2_u_0_mrcc_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,istate_i
     enddo
   enddo
   !$OMP END DO
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(static,1)
   do sh=1,shortcut(0,2)
     do i=shortcut(sh,2),shortcut(sh+1,2)-1
       org_i = sort_idx(i,2)

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -771,10 +771,8 @@ END_PROVIDER
     factor = 1.d0
     resold = huge(1.d0)
 
-    do k=0,10*hh_nex
+    do k=0,hh_nex/4
       res = 0.d0
-      !$OMP PARALLEL default(shared) private(cx, i, a_col, a_coll) reduction(+:res)
-      !$OMP DO
       do a_coll = 1, n_exc_active
         a_col = active_pp_idx(a_coll)
         cx = 0.d0
@@ -785,21 +783,20 @@ END_PROVIDER
         res = res + (X_new(a_col) - X(a_col))*(X_new(a_col) - X(a_col))
         X(a_col) = X_new(a_col)
       end do
-      !$OMP END DO
-      !$OMP END PARALLEL
       
       if (res > resold) then
         factor = factor * 0.5d0
       endif
       resold = res
       
-      if(iand(k, 4095) == 0) then
+      if(iand(k, 127) == 0) then
         print *, "res ", k, res
       end if
       
       if(res < 1d-10) exit
     end do
     dIj_unique(1:size(X), s) = X(1:size(X))
+    print *, "res ", k, res
 
   enddo
 
@@ -831,21 +828,23 @@ END_PROVIDER
         
    do s=1,N_states
      norm = 0.d0
-     double precision               :: f
+     double precision               :: f, g, gmax
+     gmax = 1.d0*maxval(dabs(psi_non_ref_coef(:,s)))
      do i=1,N_det_non_ref
-       if (rho_mrcc(i,s) == 0.d0) then
-         rho_mrcc(i,s) = 1.d-32
-       endif
-
        if (lambda_type == 2) then
          f = 1.d0
        else
+        if (rho_mrcc(i,s) == 0.d0) then
+          cycle
+        endif
         ! f is such that f.\tilde{c_i} = c_i
         f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
 
         ! Avoid numerical instabilities
-        f = min(f,2.d0)
+!        g = 1.d0+dabs(gmax / psi_non_ref_coef(i,s) )
-        f = max(f,-2.d0)
+        g = 2.d0+100.d0*exp(-20.d0*dabs(psi_non_ref_coef(i,s)/gmax))
+        f = min(f, g)
+        f = max(f,-g)
       endif
 
        norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
@@ -1087,6 +1086,22 @@ end function
   end do
   hh_shortcut(hh_shortcut(0)+1) = s+1
   
+  if (hh_shortcut(0) > N_hh_exists) then
+    print *,  'Error in ', irp_here
+    print *,  'hh_shortcut(0) :', hh_shortcut(0)
+    print *,  'N_hh_exists : ', N_hh_exists
+    print *,  'Is your active space defined?'
+    stop
+  endif
+
+  if (hh_shortcut(hh_shortcut(0)+1)-1 > N_pp_exists) then
+    print *,  'Error 1 in ', irp_here
+    print *,  'hh_shortcut(hh_shortcut(0)+1)-1 :', hh_shortcut(hh_shortcut(0)+1)-1 
+    print *,  'N_pp_exists : ', N_pp_exists
+    print *,  'Is your active space defined?'
+    stop
+  endif
+
   do s=2,4,2
     do i=1,hh_shortcut(0)
       if(hh_exists(s, i) == 0) then
@@ -1097,6 +1112,7 @@ end function
       end if
     end do
     
+
     do i=1,hh_shortcut(hh_shortcut(0)+1)-1
       if(pp_exists(s, i) == 0) then
         pp_exists(s-1, i) = 0

plugins/mrcepa0/dressing.irp.f

@@ -232,12 +232,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
         !hIk = hij_mrcc(idx_alpha(k_sd),i_I)
         !         call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(k_sd)),Nint,hIk)
         
-        do i_state=1,N_states
-          dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state)
-          !dIk(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(k_sd)), N_int) !!hIk * lambda_mrcc(i_state,idx_alpha(k_sd))
-          !dIk(i_state) = psi_non_ref_coef(idx_alpha(k_sd), i_state) / psi_ref_coef(i_I, i_state)
-        enddo
-        
         
         ! |l> = Exc(k -> alpha) |I>
         call get_excitation(psi_non_ref(1,1,idx_alpha(k_sd)),tq(1,1,i_alpha),exc,degree,phase,Nint)
@@ -250,6 +244,10 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
         call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
         if(.not. ok) cycle
         
+        do i_state=1,N_states
+          dIK(i_state) = dij(i_I, idx_alpha(k_sd), i_state)
+        enddo
+        
         ! <I| \l/ |alpha>
         do i_state=1,N_states
           dka(i_state) = 0.d0
@@ -268,12 +266,8 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
             loop = .false.
             if (.not.loop) then
               call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
-              hIl = hij_mrcc(idx_alpha(l_sd),i_I)
-!                             call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
               do i_state=1,N_states
                 dka(i_state) = dij(i_I, idx_alpha(l_sd), i_state) * phase * phase2
-                !dka(i_state) = get_dij(psi_ref(1,1,i_I), psi_non_ref(1,1,idx_alpha(l_sd)), N_int) * phase * phase2 !hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
-                !dka(i_state) = psi_non_ref_coef(idx_alpha(l_sd), i_state) / psi_ref_coef(i_I, i_state) * phase * phase2 
               enddo
             endif
 
@@ -292,7 +286,6 @@ subroutine mrcc_part_dress(delta_ij_, delta_ii_,delta_ij_s2_, delta_ii_s2_,i_gen
         k_sd = idx_alpha(l_sd)
         hla = hij_cache(k_sd)
         sla = sij_cache(k_sd)
-!        call i_h_j(tq(1,1,i_alpha),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hla)
         do i_state=1,N_states
           dIa_hla(i_state,k_sd) = dIa(i_state) * hla
           dIa_sla(i_state,k_sd) = dIa(i_state) * sla
@@ -336,6 +329,7 @@ end
   integer                        :: i, j, i_state
   
   !mrmode : 1=mrcepa0, 2=mrsc2 add, 3=mrcc
+  PROVIDE psi_ref_lock
   
   if(mrmode == 3) then
     do i = 1, N_det_ref