diff --git a/src/BiInts/mo_bi_integrals.irp.f b/src/BiInts/mo_bi_integrals.irp.f index 23565b66..9504c563 100644 --- a/src/BiInts/mo_bi_integrals.irp.f +++ b/src/BiInts/mo_bi_integrals.irp.f @@ -307,7 +307,9 @@ end &BEGIN_PROVIDER [ double precision, mo_bielec_integral_jj_anti, (mo_tot_num_align,mo_tot_num) ] implicit none BEGIN_DOC - ! Transform Bi-electronic integrals and + ! mo_bielec_integral_jj(i,j) = J_ij + ! mo_bielec_integral_jj_exchange(i,j) = J_ij + ! mo_bielec_integral_jj_anti(i,j) = J_ij - K_ij END_DOC integer :: i,j,p,q,r,s diff --git a/src/CIS_dressed/CIS_D_lapack.irp.f b/src/CIS_dressed/CIS_D_lapack.irp.f index 7ac8a229..9da91270 100644 --- a/src/CIS_dressed/CIS_D_lapack.irp.f +++ b/src/CIS_dressed/CIS_D_lapack.irp.f @@ -1,4 +1,4 @@ -program CIS_DT +program CIS_D implicit none integer :: i print*,'MP2_dresssing=',mp2_dressing @@ -18,7 +18,7 @@ endif print*,'i = ',i print*,'CIS = ',eigenvalues_CIS(i) print*,'CIS(DdT)= ',eigenvalues_CIS_dress_D(i) - print*,'s2(DdT) = ',s_2_CIS_dress_D_dt(i) + print*,'s2(DdT) = ',s_2_CIS_dress_D(i) print*,' = ',CIS_states_properties(1,i) print*,' = ',CIS_states_properties(2,i) print*,' = ',CIS_states_properties(3,i) @@ -27,39 +27,39 @@ endif print*,' = ',CIS_states_properties(6,i) print*,'' enddo - double precision :: delta_E_CIS,delta_E_CIS_DT,convert + double precision :: delta_E_CIS,delta_E_CIS_D,convert convert = 1.d0 - print*,'Excitation energies : CIS CIS_DT (Hartree)' + print*,'Excitation energies : CIS CIS_D (Hartree)' do i = 2, n_state_CIS delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1) - delta_E_CIS_DT = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) - write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + delta_E_CIS_D = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_D*convert enddo convert = 27.2114d0 - print*,'Excitation energies : CIS CIS_DT (eV)' + print*,'Excitation energies : CIS CIS_D (eV)' do i = 2, n_state_CIS delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1) - delta_E_CIS_DT = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) - write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.6,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + delta_E_CIS_D = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.6,x))')i,delta_E_CIS*convert,delta_E_CIS_D*convert enddo convert = 219475d0 - print*,'Excitation energies : CIS CIS_DT (cm-1)' + print*,'Excitation energies : CIS CIS_D (cm-1)' do i = 2, n_state_CIS delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1) - delta_E_CIS_DT = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) - write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.1,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + delta_E_CIS_D = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.1,x))')i,delta_E_CIS*convert,delta_E_CIS_D*convert enddo convert = 627.51d0 - print*,'Excitation energies : CIS CIS_DT (Kcal/mol)' + print*,'Excitation energies : CIS CIS_D (Kcal/mol)' do i = 2, n_state_CIS delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1) - delta_E_CIS_DT = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) - write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + delta_E_CIS_D = eigenvalues_CIS_dress_D(i) - eigenvalues_CIS_dress_D(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_D*convert enddo !if(save_all_dm_cis)then diff --git a/src/CIS_dressed/CIS_providers.irp.f b/src/CIS_dressed/CIS_providers.irp.f index 643c916e..9048aec3 100644 --- a/src/CIS_dressed/CIS_providers.irp.f +++ b/src/CIS_dressed/CIS_providers.irp.f @@ -212,7 +212,8 @@ allocate (delta_H_matrix_doub(size_psi_CIS,size_psi_CIS)) allocate(eigvalues(size_psi_CIS),eigvectors(size_psi_CIS,size_psi_CIS)) do i = 1,n_state_CIS - call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles +! call dress_by_doubles(eigenvalues_CIS(i),coefs_CIS(1,i),delta_H_matrix_doub,size_psi_CIS) !dressing of the Doubles + delta_H_matrix_doub = 0.d0 do j = 1,size_psi_CIS do k = 1,size_psi_CIS @@ -220,6 +221,14 @@ enddo delta_H_matrix_doub(j,j) += dress_T_discon_array_CIS(j) enddo + + double precision :: accu + accu = 0.d0 + do j = 1, size_psi_CIS + do k = 1, size_psi_CIS + accu += delta_H_matrix_doub(j,k) * coefs_CIS(j,i) * coefs_CIS(k,i) + enddo + enddo call lapack_diag(eigvalues,eigvectors,delta_H_matrix_doub,size_psi_CIS,size_psi_CIS) ! state following @@ -235,7 +244,9 @@ endif ! enddo - ! print*,'overlap = ',max_overlap + print*,i,i_overlap + print*,'overlap = ',max_overlap + i_overlap = i overlap_Ddt=max_overlap do k = 1,size_psi_CIS eigenvectors_CIS_dress_D_dt(k,i) = eigvectors(k,i_overlap) @@ -248,6 +259,9 @@ call get_s2_u0(psi_CIS,eigenvectors_CIS_dress_D_dt(1,i),size_psi_CIS,size_psi_CIS,s2) s_2_CIS_dress_D_dt(i) = s2 eigenvalues_CIS_dress_D_dt(i) = eigvalues(i_overlap) + print*,'eigenvalues_CIS_dress_D_dt(i)= ',eigenvalues_CIS_dress_D_dt(i) + print*,'accu = ',accu + print*,'eigenvalues_CIS = ',eigenvalues_CIS(i) enddo END_PROVIDER diff --git a/src/CIS_dressed/EN2.irp.f b/src/CIS_dressed/EN2.irp.f new file mode 100644 index 00000000..d89fbb9f --- /dev/null +++ b/src/CIS_dressed/EN2.irp.f @@ -0,0 +1,252 @@ + + BEGIN_PROVIDER [double precision, EN2_corr_energy] +&BEGIN_PROVIDER [double precision, p_imp_EN,(elec_alpha_num+1:n_act_cis)] +&BEGIN_PROVIDER [double precision, h_imp_EN,(n_core_cis+1:elec_alpha_num)] +&BEGIN_PROVIDER [double precision, hp_imp_EN,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] + + BEGIN_DOC + !Calculation of the EN2 correlation energy (EN2_corr_energy) + !and calculation of the contribution of the disconnected Triples on the + !Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN) + END_DOC + + implicit none + integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) + double precision :: direct,exchg,hij !calculating direct, exchange and total contribution + double precision :: get_mo_bielec_integral + double precision :: e_i,e_k,e_j,e_l !epsilons of i,j,k,l + double precision :: delta_e_ik,delta_e_ikj + double precision :: delta_e !delta epsilons + double precision :: delta_e_tmp,H_jj_total + integer :: ispin1,ispin2,i_ok + + print*,'EN2_corr_energy' + + EN2_corr_energy=0.d0 + + do i=n_core_cis+1,elec_alpha_num + h_imp_EN(i)=0.d0 + do k =elec_beta_num + 1, n_act_cis + p_imp_EN(k)=0.d0 + hp_imp_EN(i,k)=0.d0 + enddo + enddo + print*,'HF_energy = ',HF_energy + print*,'1' + + if(EN_2_2)then + do i=n_core_cis+1,elec_alpha_num + h_imp_EN(i)=0.d0 + + e_i=diagonal_Fock_matrix_mo(i) + + do k=elec_alpha_num+1,n_act_cis + hp_imp_EN(i,k)=0.d0 + + e_k=diagonal_Fock_matrix_mo(k) + delta_e_ik=e_i-e_k + + !same spin contribution for EN2_corr_energy + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) + ! + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + + EN2_corr_energy=EN2_corr_energy+2*hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + + p_imp_EN(l)=p_imp_EN(l)+hij + h_imp_EN(j)=h_imp_EN(j)+hij + + enddo + enddo + + !same spin contribution for hp_imp_EN + + do j=n_core_cis+1,elec_alpha_num + if(j==i)cycle + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_alpha_num+1,n_act_cis + if(l==k)cycle + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij = 0.5d0 * (-delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij)) + + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj(i,j) - & + mo_bielec_integral_jj(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(j,l) + & + mo_bielec_integral_jj(i,l) + & + mo_bielec_integral_jj(j,k) + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + + hij=direct*direct + hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) + + EN2_corr_energy=EN2_corr_energy+hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + enddo + enddo + print*,'2' + + else + do i=n_core_cis+1,elec_alpha_num + h_imp_EN(i)=0.d0 + + e_i=diagonal_Fock_matrix_mo(i) + + do k=elec_alpha_num+1,n_act_cis + hp_imp_EN(i,k)=0.d0 + + e_k=diagonal_Fock_matrix_mo(k) + delta_e_ik=e_i-e_k + + !same spin contribution for EN2_corr_energy + do j=i+1,elec_alpha_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN + do l=k+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) + ! + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij = hij*hij/delta_e + + EN2_corr_energy=EN2_corr_energy+2*hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + + p_imp_EN(l)=p_imp_EN(l)+hij + h_imp_EN(j)=h_imp_EN(j)+hij + + enddo + enddo + + !same spin contribution for hp_imp_EN + + do j=n_core_cis+1,elec_alpha_num + if(j==i)cycle + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_alpha_num+1,n_act_cis + if(l==k)cycle + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & + mo_bielec_integral_jj_anti(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(i,l) + & + mo_bielec_integral_jj_anti(j,k) + & + mo_bielec_integral_jj_anti(j,l) + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + + hij=(direct-exchg)*(direct-exchg) + hij = hij*hij/delta_e + + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + + !different spin contribution + do j=n_core_cis+1,elec_beta_num + e_j=diagonal_Fock_matrix_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=diagonal_Fock_matrix_mo(l) + delta_e=delta_e_ikj-e_l + delta_e=delta_e-mo_bielec_integral_jj(i,j) - & + mo_bielec_integral_jj(k,l) + delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & + mo_bielec_integral_jj_anti(j,l) + & + mo_bielec_integral_jj(i,l) + & + mo_bielec_integral_jj(j,k) + + direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + + hij=direct*direct + hij = hij*hij/delta_e + + EN2_corr_energy=EN2_corr_energy+hij + p_imp_EN(k)=p_imp_EN(k)+hij + h_imp_EN(i)=h_imp_EN(i)+hij + hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij + + enddo + enddo + enddo + enddo + + endif + + + print*,'EN correlation energy=',EN2_corr_energy + print*,'EN correlation energy=',EN2_corr_energy + + END_PROVIDER + + diff --git a/src/CIS_dressed/MP2.irp.f b/src/CIS_dressed/MP2.irp.f index b8a6931c..72012fde 100644 --- a/src/CIS_dressed/MP2.irp.f +++ b/src/CIS_dressed/MP2.irp.f @@ -139,269 +139,6 @@ END_PROVIDER - BEGIN_PROVIDER [double precision, EN2_corr_energy] - &BEGIN_PROVIDER [double precision, p_imp_EN,(elec_alpha_num+1:n_act_cis)] - &BEGIN_PROVIDER [double precision, h_imp_EN,(n_core_cis+1:elec_alpha_num)] - &BEGIN_PROVIDER [double precision, hp_imp_EN,(n_core_cis+1:elec_alpha_num,elec_alpha_num+1:n_act_cis)] - - BEGIN_DOC - !Calculation of the EN2 correlation energy (EN2_corr_energy) - !and calculation of the contribution of the disconnected Triples on the - !Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN) - END_DOC - - implicit none - integer :: i,j,k,l !variables for going over the occupied (i,j) and virutal (k,l) - double precision :: direct,exchg,hij !calculating direct, exchange and total contribution - double precision :: get_mo_bielec_integral - double precision :: e_i,e_k,e_j,e_l !epsilons of i,j,k,l - double precision :: delta_e_ik,delta_e_ikj - double precision :: delta_e !delta epsilons - double precision :: delta_e_tmp,H_jj_total - integer, allocatable :: key_in(:,:) - integer, allocatable :: key_out(:,:) - integer :: ispin1,ispin2,i_ok - allocate(key_in(N_int,2)) - allocate(key_out(N_int,2)) - - print*,'EN2_corr_energy' - - EN2_corr_energy=0.d0 - -do i=n_core_cis+1,elec_alpha_num - h_imp_EN(i)=0.d0 - do k =elec_beta_num + 1, n_act_cis - p_imp_EN(k)=0.d0 - hp_imp_EN(i,k)=0.d0 - enddo -enddo - - do i=n_core_cis+1,elec_alpha_num -! h_imp_EN(i)=0.d0 - - e_i=diagonal_Fock_matrix_mo(i) - - do k=elec_alpha_num+1,n_act_cis -! hp_imp_EN(i,k)=0.d0 - - e_k=diagonal_Fock_matrix_mo(k) - delta_e_ik=e_i-e_k - - !same spin contribution for EN2_corr_energy - do j=i+1,elec_alpha_num - e_j=diagonal_Fock_matrix_mo(j) - delta_e_ikj=delta_e_ik+e_j - - !same spin contribution for EN2_corr_energy and a part of the contribution to p_imp_EN,h_imp_EN - do l=k+1,n_act_cis - e_l=diagonal_Fock_matrix_mo(l) - delta_e=delta_e_ikj-e_l - delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & - mo_bielec_integral_jj_anti(k,l) - delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & - mo_bielec_integral_jj_anti(i,l) + & - mo_bielec_integral_jj_anti(j,k) + & - mo_bielec_integral_jj_anti(j,l) - ! ispin1 = 1 - ! ispin2 = 1 - ! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) - ! delta_e_tmp = HF_energy - H_jj_total(key_out) - ! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then - ! print*,'same spin first' - ! print*,'delta_e_SSS = ',delta_e_tmp-delta_e - ! stop - ! endif - ! - direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) - - hij=(direct-exchg)*(direct-exchg) - hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) - ! ispin1 = 1 - ! ispin2 = 1 - ! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) - ! delta_e_tmp = HF_energy - H_jj_total(key_out) - - EN2_corr_energy=EN2_corr_energy+2*hij - p_imp_EN(k)=p_imp_EN(k)+hij - h_imp_EN(i)=h_imp_EN(i)+hij - - p_imp_EN(l)=p_imp_EN(l)+hij - h_imp_EN(j)=h_imp_EN(j)+hij - - enddo -enddo -! !remaining same spin contribution for p_imp_EN - -! do l=elec_alpha_num+1,k-1 -! e_l=diagonal_Fock_matrix_mo(l) -! delta_e=delta_e_ikj-e_l -! delta_e=e_i + e_j - e_k -e_l -! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & -! mo_bielec_integral_jj_anti(k,l) -! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & -! mo_bielec_integral_jj_anti(i,l) + & -! mo_bielec_integral_jj_anti(j,k) + & -! mo_bielec_integral_jj_anti(j,l) -! e_l=diagonal_Fock_matrix_mo(l) -! delta_e=delta_e_ikj-e_l -! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & -! mo_bielec_integral_jj_anti(k,l) -! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & -! mo_bielec_integral_jj_anti(i,l) + & -! mo_bielec_integral_jj_anti(j,k) + & -! mo_bielec_integral_jj_anti(j,l) - - - - - -!! ispin1 = 1 -!! ispin2 = 1 -!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) -!! delta_e_tmp = HF_energy - H_jj_total(key_out) -!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then -!! call print_key(key_out) -!! call print_key(ref_bitmask) -!! print*,i,k,j,l -!! print*,'same spin middle' -!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e -!! print*,'delta_e_SSS = ',delta_e_tmp,delta_e -!! stop -!! endif - -! direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) -! exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) - -! hij=(direct-exchg)*(direct-exchg) -! hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) - -! p_imp_EN(k)=p_imp_EN(k)+hij - -! enddo -! enddo - -! !remaining same spin contribution for h_imp_EN -! do j=n_core_cis+1,i-1 -! e_j=diagonal_Fock_matrix_mo(j) -! delta_e_ikj=delta_e_ik+e_j - -! do l=k+1,n_act_cis -! e_l=diagonal_Fock_matrix_mo(l) -! delta_e=delta_e_ikj-e_l -! delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & -! mo_bielec_integral_jj_anti(k,l) -! delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & -! mo_bielec_integral_jj_anti(i,l) + & -! mo_bielec_integral_jj_anti(j,k) + & -! mo_bielec_integral_jj_anti(j,l) -!! ispin2 = 1 -!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) -!! delta_e_tmp = HF_energy - H_jj_total(key_out) -!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then -!! print*,'same spin third ' -!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e -!! stop -!! endif - -! direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) -! exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) - -! hij=(direct-exchg)*(direct-exchg) -! hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) - -! h_imp_EN(i)=h_imp_EN(i)+hij - -! enddo -! enddo - - !same spin contribution for hp_imp_EN - - do j=n_core_cis+1,elec_alpha_num - if(j==i)cycle - e_j=diagonal_Fock_matrix_mo(j) - delta_e_ikj=delta_e_ik+e_j - - do l=elec_alpha_num+1,n_act_cis - if(l==k)cycle - e_l=diagonal_Fock_matrix_mo(l) - delta_e=delta_e_ikj-e_l - delta_e=delta_e-mo_bielec_integral_jj_anti(i,j) - & - mo_bielec_integral_jj_anti(k,l) - delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & - mo_bielec_integral_jj_anti(i,l) + & - mo_bielec_integral_jj_anti(j,k) + & - mo_bielec_integral_jj_anti(j,l) -!!! ispin1 = 1 -!!! ispin2 = 1 -!!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) -!!! delta_e_tmp = HF_energy - H_jj_total(key_out) -!!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then -!!! print*,'same spin fourth' -!!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e -!!! call print_key(key_out) -!!! call print_key(ref_bitmask) -!!! print*,i,k,j,l -!!! stop -!!! endif - - direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - exchg=get_mo_bielec_integral(i,j,l,k,mo_integrals_map) - - hij=(direct-exchg)*(direct-exchg) - hij = 0.5d0 * (-delta_e - dsqrt(delta_e * delta_e + 4.d0 * hij)) - - hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij - - enddo - enddo - - !different spin contribution - do j=n_core_cis+1,elec_beta_num - e_j=diagonal_Fock_matrix_mo(j) - delta_e_ikj=delta_e_ik+e_j - - do l=elec_beta_num+1,n_act_cis - e_l=diagonal_Fock_matrix_mo(l) - delta_e=delta_e_ikj-e_l - delta_e=delta_e-mo_bielec_integral_jj(i,j) - & - mo_bielec_integral_jj(k,l) - delta_e=delta_e+mo_bielec_integral_jj_anti(i,k) + & - mo_bielec_integral_jj_anti(j,l) + & - mo_bielec_integral_jj(i,l) + & - mo_bielec_integral_jj(j,k) -!! ispin1 = 2 -!! ispin2 = 1 -!! call diexcitation(i,j,k,l,ispin1,ispin2,ref_bitmask,key_out,i_ok,N_int) -!! delta_e_tmp = HF_energy - H_jj_total(key_out) -!! if(dabs(delta_e_tmp-delta_e).gt.1.d-10)then -!! print*,'different spin ' -!! print*,'delta_e_SSS = ',delta_e_tmp-delta_e -!! stop -!! endif - - - direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) - - hij=direct*direct - hij=0.5d0*(-delta_e-dsqrt(delta_e*delta_e+4.d0*hij)) - - EN2_corr_energy=EN2_corr_energy+hij - p_imp_EN(k)=p_imp_EN(k)+hij - h_imp_EN(i)=h_imp_EN(i)+hij - hp_imp_EN(i,k)=hp_imp_EN(i,k)+hij - - enddo - enddo - enddo - enddo - - print*,'EN correlation energy=',EN2_corr_energy - print*,'EN correlation energy=',EN2_corr_energy - - END_PROVIDER - - BEGIN_PROVIDER [integer, size_psi_CIS] @@ -630,6 +367,7 @@ enddo integer :: key !key for CIS-matrix print*,'dress_T_discon' + print*,'mp2_dressing = ',mp2_dressing if(mp2_dressing)then dress_T_discon_array_CIS(1)=MP2_corr_energy @@ -647,24 +385,25 @@ enddo enddo else !EN Dressing + print*,'coucou !' dress_T_discon_array_CIS(1)=EN2_corr_energy + print*,'coucou !' - do i=n_core_cis+1,elec_alpha_num - do k=elec_alpha_num+1,n_act_cis - key=psi_CIS_adress(i,k) - - dress_T_discon(i,k)=EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k) +! do i=n_core_cis+1,elec_alpha_num +! print*,'i',i,n_core_cis +! do k=elec_alpha_num+1,n_act_cis +! print*,'k',k,n_act_cis +! key=psi_CIS_adress(i,k) +! +! dress_T_discon(i,k)=EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k) -!print*,'i,k,key',i,k,key -!print*,'EN2_corr_energy-p_imp_EN(k)-h_imp_EN(i)+hp_imp_EN(i,k)',EN2_corr_energy,p_imp_EN(k),h_imp_EN(i),hp_imp_EN(i,k) -!print*,'dress_T_discon',dress_T_discon(i,k) +! dress_T_discon_array_CIS(key) = dress_T_discon(i,k) +! dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k) - dress_T_discon_array_CIS(key) = dress_T_discon(i,k) - dress_T_discon_array_CIS(key+1) = dress_T_discon(i,k) - - enddo - enddo +! enddo +! enddo end if + print*,'end' END_PROVIDER diff --git a/src/CIS_dressed/cis_dressed.ezfio_config b/src/CIS_dressed/cis_dressed.ezfio_config index fd9f291a..3dda6ffe 100644 --- a/src/CIS_dressed/cis_dressed.ezfio_config +++ b/src/CIS_dressed/cis_dressed.ezfio_config @@ -4,3 +4,4 @@ cis_dressed n_act_cis integer mp2_dressing logical standard_doubles logical + en_2_2 logical diff --git a/src/CIS_dressed/options.irp.f b/src/CIS_dressed/options.irp.f index 50872abc..15643737 100644 --- a/src/CIS_dressed/options.irp.f +++ b/src/CIS_dressed/options.irp.f @@ -83,3 +83,21 @@ BEGIN_PROVIDER [ logical , standard_doubles] END_PROVIDER + +BEGIN_PROVIDER [ logical , en_2_2] + implicit none + BEGIN_DOC + ! Number of states asked for the CIS vector + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_cis_dressed_en_2_2(has) + if (has) then + call ezfio_get_cis_dressed_en_2_2(en_2_2) + else + en_2_2 = .False. + endif + +END_PROVIDER + diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 1cb5d4e7..9474b4f9 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -385,7 +385,10 @@ subroutine $subroutine($params_main) nmax = ( N_det_generators/nproc ) *nproc call wall_time(wall_1) + + !$ call omp_init_lock(lck) +IRP_IF I_LIKE_BUGS !$OMP PARALLEL DEFAULT(SHARED) & !$OMP PRIVATE(i_generator,wall_2,ispin,k,mask) allocate( mask(N_int,2,6) ) @@ -447,6 +450,13 @@ subroutine $subroutine($params_main) allocate( mask(N_int,2,6) ) do i_generator=nmax+1,N_det_generators + +IRP_ELSE + allocate( mask(N_int,2,6) ) + do i_generator=1,N_det_generators + +IRP_ENDIF + if (abort_here) then exit endif diff --git a/src/Full_CI/H_apply.irp.f b/src/Full_CI/H_apply.irp.f index 76c046ae..0bf120c0 100644 --- a/src/Full_CI/H_apply.irp.f +++ b/src/Full_CI/H_apply.irp.f @@ -3,11 +3,11 @@ BEGIN_SHELL [ /usr/bin/env python ] from generate_h_apply import * s = H_apply("FCI") -s.set_selection_pt2("epstein_nesbet_2x2") +s.set_selection_pt2("epstein_nesbet") print s s = H_apply("FCI_PT2") -s.set_perturbation("epstein_nesbet_2x2") +s.set_perturbation("epstein_nesbet") print s END_SHELL diff --git a/src/Full_CI/full_ci.ezfio_config b/src/Full_CI/full_ci.ezfio_config index 9d984059..022f40ac 100644 --- a/src/Full_CI/full_ci.ezfio_config +++ b/src/Full_CI/full_ci.ezfio_config @@ -1,3 +1,4 @@ full_ci n_det_max_fci integer + pt2_max double precision do_pt2_end logical diff --git a/src/Full_CI/full_ci.irp.f b/src/Full_CI/full_ci.irp.f index 6721bc0a..d1dbc2f5 100644 --- a/src/Full_CI/full_ci.irp.f +++ b/src/Full_CI/full_ci.irp.f @@ -12,7 +12,7 @@ program cisd pt2 = 1.d0 diag_algorithm = "Lapack" ! do while (maxval(abs(pt2(1:N_st))) > 1.d-4) - do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > 1.d-4) + do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max) call H_apply_FCI(pt2, norm_pert, H_pert_diag, N_st) call diagonalize_CI call save_wavefunction @@ -29,7 +29,7 @@ program cisd N_det = min(n_det_max_fci,N_det) if(do_pt2_end)then threshold_selectors = 1.d0 - threshold_generators = 0.999d0 + threshold_generators = 0.99d0 touch N_det psi_det psi_coef call diagonalize_CI call H_apply_FCI_PT2(pt2, norm_pert, H_pert_diag, N_st) diff --git a/src/Hartree_Fock/diagonalize_fock.irp.f b/src/Hartree_Fock/diagonalize_fock.irp.f index c5faea9c..d2e6d538 100644 --- a/src/Hartree_Fock/diagonalize_fock.irp.f +++ b/src/Hartree_Fock/diagonalize_fock.irp.f @@ -54,3 +54,28 @@ deallocate(work, iwork, F, S) END_PROVIDER +BEGIN_PROVIDER [double precision, diagonal_Fock_matrix_mo_sum, (mo_tot_num)] + implicit none + BEGIN_DOC + ! diagonal element of the fock matrix calculated as the sum over all the interactions + ! with all the electrons in the RHF determinant + ! diagonal_Fock_matrix_mo_sum(i) = sum_{j=1, N_elec} 2 J_ij -K_ij + END_DOC + integer :: i,j + double precision :: accu + do i = 1,elec_alpha_num + accu = 0.d0 + do j = 1, elec_alpha_num + accu += 2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i) + enddo + do i = elec_alpha_num+1,mo_tot_num + accu = 0.d0 + do j = 1, elec_alpha_num + accu += 2.d0 * mo_bielec_integral_jj(i,j) - mo_bielec_integral_jj_exchange(i,j) + enddo + diagonal_Fock_matrix_mo_sum(i) = accu + mo_mono_elec_integral(i,i) + enddo + +END_PROVIDER diff --git a/src/Utils/util.irp.f b/src/Utils/util.irp.f index b2618a60..7d81d49b 100644 --- a/src/Utils/util.irp.f +++ b/src/Utils/util.irp.f @@ -296,6 +296,27 @@ subroutine normalize(u,sze) endif end +double precision function approx_dble(a,n) + implicit none + integer, intent(in) :: n + double precision, intent(in) :: a + double precision :: f + integer :: i + + if (a == 0.d0) then + approx_dble = 0.d0 + return + endif + f = 1.d0 + do i=1,-int(dlog10(dabs(a)))+n + f = f*.1d0 + enddo + do i=1,int(dlog10(dabs(a)))-n + f = f*10.d0 + enddo + approx_dble = dnint(a/f)*f + +end