Frozen core is working.

src/csf/sigma_vector.irp.f

@@ -837,6 +837,7 @@ subroutine calculate_preconditioner_cfg(diag_energies)
   real*8            :: hpp
   real*8            :: meCC
   real*8            :: ecore
+  real*8            :: core_act_contrib
 
   !PROVIDE h_core_ri
   PROVIDE core_fock_operator
@@ -863,6 +864,8 @@ subroutine calculate_preconditioner_cfg(diag_energies)
      starti = psi_config_data(i,1)
      endi   = psi_config_data(i,2)
 
+     core_act_contrib = 0.0d0
+
      ! find out all pq holes possible
      nholes = 0
      ! holes in SOMO
@@ -915,6 +918,13 @@ subroutine calculate_preconditioner_cfg(diag_energies)
         p = listholes(k)
         noccp = holetype(k)
 
+
+        ! core-active
+        do l = 1, n_core_orb
+         jj = list_core(l)
+         core_act_contrib += noccp * (2.d0 * mo_two_e_integrals_jj(jj,p) - mo_two_e_integrals_jj_exchange(jj,p))
+        enddo
+
         ! Calculate one-electron
         ! and two-electron coulomb terms
         do l=1,nholes
@@ -944,6 +954,10 @@ subroutine calculate_preconditioner_cfg(diag_energies)
         enddo
 
      enddo
+     !print *,"I=",i," core_act=",core_act_contrib
+     do j=starti,endi
+       diag_energies(j) += core_act_contrib
+     end do
   enddo
 
 end subroutine calculate_preconditioner_cfg
@@ -1345,6 +1359,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   real*8, external               :: get_two_e_integral
   real*8,dimension(:),allocatable:: diag_energies
   real*8                         :: tmpvar, tmptot
+  real*8                         :: core_act_contrib
 
   integer(omp_lock_kind), allocatable :: lock(:)
   call omp_set_max_active_levels(1)
@@ -1404,9 +1419,10 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
       !$OMP   diagfac, tmpvar, diagfactors_0)                                            &
       !$OMP shared(istart_cfg, iend_cfg, psi_configuration, mo_num, psi_config_data,&
       !$OMP    N_int, N_st, psi_out, psi_in, h_core_ri, core_energy, h_act_ri, AIJpqContainer,&
-      !$OMP     sze, NalphaIcfg_list,alphasIcfg_list, bit_tmp,       &
+      !$OMP     pp, sze, NalphaIcfg_list,alphasIcfg_list, bit_tmp,       &
       !$OMP     AIJpqMatrixDimsList, diag_energies, n_CSF, lock, NBFmax,nconnectedtotalmax, nconnectedmaxJ,maxnalphas,&
-      !$OMP     n_core_orb, n_act_orb, list_act, num_threads_max)
+      !$OMP     n_core_orb, n_act_orb, list_act, n, list_core,  list_core_is_built,core_act_contrib, num_threads_max,&
+      !$OMP     n_core_orb_is_built, mo_integrals_map, mo_integrals_map_is_built)
 
   allocate(singlesI(N_INT,2,max(sze,10000)))
   allocate(idxs_singlesI(max(sze,10000)))
@@ -1520,8 +1536,18 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
         do jj = startj, endj
           cntj = jj-startj+1
           !meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* h_core_ri(p,q)
-          meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* h_act_ri(p,q)
+          core_act_contrib = 0.0d0
+          if(p.ne.q)then
+          do pp=1,n_core_orb
+            n=list_core(pp)
+            core_act_contrib += 2.d0 * get_two_e_integral(p,n,q,n,mo_integrals_map) - get_two_e_integral(p,n,n,q,mo_integrals_map)
+          end do
+          endif
+          meCC1 = AIJpqContainer(cnti,cntj,pmodel,qmodel,extype,NSOMOI)* (h_act_ri(p,q) + core_act_contrib)
           !print *,"jj = ",jj
+          !if(ii.eq.1 .and. jj.eq.177 )then
+          !  print *,"p=",p," q=",q," hact=",h_act_ri(p,q), " core_act=",core_act_contrib
+          !endif
           call omp_set_lock(lock(jj))
           do kk = 1,n_st
             psi_out(kk,jj) = psi_out(kk,jj) + meCC1 * psi_in(kk,ii)
@@ -1545,7 +1571,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   deallocate(excitationIds_single)
   deallocate(excitationTypes_single)
 
-  !print *," psi(60,1)=",psi_out(1,60)
+  !print *," singles part psi(1,177)=",psi_out(1,177)
   
   allocate(listconnectedJ(N_INT,2,max(sze,10000)))
   allocate(alphas_Icfg(N_INT,2,max(sze,10000)))
@@ -1691,6 +1717,9 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
 
            do m = 1,colsikpq
               call omp_set_lock(lock(idxs_connectedI_alpha(j)+m-1))
+              !if((idxs_connectedI_alpha(j)+m-1).eq.177)then
+              !  print *,"CC=",CCmattmp(1,m)
+              !endif
               do kk = 1,n_st
                  psi_out(kk,idxs_connectedI_alpha(j)+m-1) += CCmattmp(kk,m)
               enddo
@@ -1714,6 +1743,7 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
   deallocate(excitationTypes)
   deallocate(diagfactors)
 
+  !print *," psi(1,177)=",psi_out(1,177)
 
   ! Add the diagonal contribution
   !$OMP DO
@@ -1726,7 +1756,6 @@ subroutine calculate_sigma_vector_cfg_nst_naive_store(psi_out, psi_in, n_st, sze
 
   !$OMP END PARALLEL
   call omp_set_max_active_levels(4)
-  !print *," psi(60,1)=",psi_out(1,60)
   !print *," diag_enregy=",diag_energies(1), " psi_out(1,1)=",psi_out(1,1)
 
   deallocate(diag_energies)

src/davidson/diagonalization_hcfg.irp.f

@@ -112,6 +112,7 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
   double precision, allocatable  :: U(:,:), U_csf(:,:), overlap(:,:)
   double precision, allocatable  :: tmpU(:,:), tmpW(:,:)
   double precision, pointer      :: W(:,:), W_csf(:,:)
+  double precision, pointer      :: W2(:,:), U2(:,:)
   logical                        :: disk_based
   double precision               :: energy_shift(N_st_diag_in*davidson_sze_max)
 
@@ -234,11 +235,13 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
     call c_f_pointer(ptr_w, W_csf, (/sze_csf,N_st_diag*itermax/))
   else
     allocate(W(sze,N_st_diag),W_csf(sze_csf,N_st_diag*itermax))
+    allocate(W2(sze,N_st_diag))
   endif
 
   allocate(                                                          &
       ! Large
       U(sze,N_st_diag),                                              &
+      U2(sze,N_st_diag),                                              &
       U_csf(sze_csf,N_st_diag*itermax),                              &
 
       ! Small
@@ -324,6 +327,10 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
                 tmpU(kk,ii) = U_csf(ii,shift+kk)
               enddo
             enddo
+            !do j=1,1
+            !  print *,"====> J=",j
+            !tmpU=0.0d0
+            !tmpU(1,j)=1.0d0
             call calculate_sigma_vector_cfg_nst_naive_store(tmpW,tmpU,N_st_diag,sze_csf,1,sze_csf,0,1)
             do kk=1,N_st_diag
               do ii=1,sze_csf
@@ -331,6 +338,20 @@ subroutine davidson_diag_cfg_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,sze_csf,N
               enddo
             enddo
 
+            !U_csf=0.0d0
+            !U_csf(j,1)=1.0d0
+            !call convertWFfromCSFtoDET(N_st_diag,U_csf(1,1),U(1,1))
+            !call H_u_0_nstates_openmp(U2,U,N_st_diag,sze)
+            !call convertWFfromDETtoCSF(N_st_diag,U2(1,1),W2(1,1))
+            !print *," w2=",W2(j,1)
+            !do i=1,sze_csf
+            !  print *, " i=",i,"qp=",W2(i,1)," my=",W_csf(i,1),dabs(dabs(W2(i,1))-dabs(W_csf(i,1)))
+            !  if(dabs(dabs(W2(i,1))-dabs(W_csf(i,1))) .ge. 1.0e-10)then
+            !    print *," somo=",psi_configuration(1,1,i)," domo=",psi_configuration(1,2,i)
+            !  endif
+            !end do
+            !end do
+            !stop
             deallocate(tmpW)
             deallocate(tmpU)
         endif

src/davidson/diagonalize_ci.irp.f

@@ -72,7 +72,9 @@ END_PROVIDER
    if (diag_algorithm == "Davidson") then
 
      if (do_csf) then
+     !if (.true.) then
        if (sigma_vector_algorithm == 'det') then
+       !if (.false.) then
          call davidson_diag_H_csf(psi_det,CI_eigenvectors, &
            size(CI_eigenvectors,1),CI_electronic_energy,               &
            N_det,N_csf,min(N_det,N_states),min(N_det,N_states_diag),N_int,0,converged)

src/mo_two_e_ints/core_quantities.irp.f

@@ -13,19 +13,7 @@ BEGIN_PROVIDER [double precision, core_energy]
    core_energy += 2.d0 * (2.d0 * mo_two_e_integrals_jj(j,l) - mo_two_e_integrals_jj_exchange(j,l))
   enddo
  enddo
- !print *," core no nucl=",core_energy
- ! core-active
- do i = 1, n_core_orb
-  j = list_core(i)
-  !!! VJ
-  !!! TODO: Correct the loop over active electrons
-  do k = 1, (elec_num - 2*n_core_orb)/2
-   l = k + n_core_orb
-   core_energy += 2.d0 * (2.d0 * mo_two_e_integrals_jj(j,l) - mo_two_e_integrals_jj_exchange(j,l))
-  enddo
- enddo
  core_energy += nuclear_repulsion
- !print *," core no nucl=",core_energy
 
 END_PROVIDER
 
@@ -84,6 +72,7 @@ BEGIN_PROVIDER [ double precision, h_act_ri, (mo_num, mo_num) ]
  integer :: p,q, r
  ! core-core contribution
  h_act_ri = core_fock_operator
+ !print *,' Bef----hact(1,14)=',h_act_ri(4,14)
  ! act-act contribution
  do p=1,n_act_orb
    j=list_act(p)
@@ -100,14 +89,15 @@ BEGIN_PROVIDER [ double precision, h_act_ri, (mo_num, mo_num) ]
    enddo
  enddo
  ! core-act contribution
- do p=1,n_core_orb
+ !do p=1,n_act_orb
-   j=list_core(p)
+ !  j=list_core(p)
-   do k=1,mo_num
+ !  do k=1,n_core_orb
-     do q=1,n_act_orb
+ !    do q=1,n_act_orb
-       i=list_act(q)
+ !      i=list_act(q)
-       h_act_ri(i,j) = h_act_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j)
+ !      h_act_ri(i,j) = h_act_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j)
-     enddo
+ !    enddo
-   enddo
+ !  enddo
- enddo
+ !enddo
+ !print *,' Aft----hact(1,14)=',h_act_ri(4,14), mo_one_e_integrals(4,14)
 END_PROVIDER