From 92be856e50c5fdeee0d058434b7338d97d9bc4b5 Mon Sep 17 00:00:00 2001 From: v1j4y Date: Fri, 17 Jun 2022 22:05:54 +0200 Subject: [PATCH] Frozen core is working. --- src/csf/sigma_vector.irp.f | 39 +++++++++++++++++++++---- src/davidson/diagonalization_hcfg.irp.f | 21 +++++++++++++ src/davidson/diagonalize_ci.irp.f | 2 ++ src/mo_two_e_ints/core_quantities.irp.f | 32 +++++++------------- 4 files changed, 68 insertions(+), 26 deletions(-) diff --git a/src/csf/sigma_vector.irp.f b/src/csf/sigma_vector.irp.f index 4420faec..8f9f99d3 100644 --- a/src/csf/sigma_vector.irp.f +++ b/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) diff --git a/src/davidson/diagonalization_hcfg.irp.f b/src/davidson/diagonalization_hcfg.irp.f index bde2bdd9..193d103b 100644 --- a/src/davidson/diagonalization_hcfg.irp.f +++ b/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 diff --git a/src/davidson/diagonalize_ci.irp.f b/src/davidson/diagonalize_ci.irp.f index 78478c57..f5045913 100644 --- a/src/davidson/diagonalize_ci.irp.f +++ b/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) diff --git a/src/mo_two_e_ints/core_quantities.irp.f b/src/mo_two_e_ints/core_quantities.irp.f index e048f6a6..ef99bd5a 100644 --- a/src/mo_two_e_ints/core_quantities.irp.f +++ b/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 - j=list_core(p) - do k=1,mo_num - do q=1,n_act_orb - i=list_act(q) - h_act_ri(i,j) = h_act_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j) - enddo - enddo - enddo + !do p=1,n_act_orb + ! j=list_core(p) + ! do k=1,n_core_orb + ! do q=1,n_act_orb + ! i=list_act(q) + ! h_act_ri(i,j) = h_act_ri(i,j) - 0.5 * big_array_exchange_integrals(k,i,j) + ! enddo + ! enddo + !enddo + !print *,' Aft----hact(1,14)=',h_act_ri(4,14), mo_one_e_integrals(4,14) END_PROVIDER