diff --git a/src/CIS_dressed/CIS_DT_lapack.irp.f b/src/CIS_dressed/CIS_DT_lapack.irp.f new file mode 100644 index 00000000..9d8ba522 --- /dev/null +++ b/src/CIS_dressed/CIS_DT_lapack.irp.f @@ -0,0 +1,68 @@ +program CIS_DT + implicit none + integer :: i + print*,'MP2_dresssing=',mp2_dressing + print*,'standard_doubles=',standard_doubles + print*,'n_state_CIS=',n_state_CIS + print*,'n_core_cis=',n_core_cis + print*,'n_act_cis=',n_act_cis + print*,'' + print*,'nuc repulsion E=',nuclear_repulsion +if (mp2_dressing==.true.) then + print*,'correlation E=',MP2_corr_energy +else + print*,'correlation E=',EN2_corr_energy +endif + do i = 1,n_state_CIS + print*,'' + print*,'i = ',i + print*,'CIS = ',eigenvalues_CIS(i) + print*,'CIS(DdT)= ',eigenvalues_CIS_dress_D_dt(i) + print*,'s2(DdT) = ',s_2_CIS_dress_D_dt(i) + print*,' = ',CIS_states_properties(1,i) + print*,' = ',CIS_states_properties(2,i) + print*,' = ',CIS_states_properties(3,i) + print*,' = ',CIS_states_properties(4,i) + print*,' = ',CIS_states_properties(5,i) + print*,' = ',CIS_states_properties(6,i) + print*,'' + enddo + double precision :: delta_E_CIS,delta_E_CIS_DT,convert + + convert = 1.d0 + print*,'Excitation energies : CIS CIS_DT (Hartree)' + do i = 2, n_state_CIS + delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1) + delta_E_CIS_DT = eigenvalues_CIS_dress_D_dt(i) - eigenvalues_CIS_dress_D_dt(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + enddo + + convert = 27.2114d0 + print*,'Excitation energies : CIS CIS_DT (eV)' + do i = 2, n_state_CIS + delta_E_CIS = eigenvalues_CIS(i) - eigenvalues_CIS(1) + delta_E_CIS_DT = eigenvalues_CIS_dress_D_dt(i) - eigenvalues_CIS_dress_D_dt(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.6,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + enddo + + + convert = 219475d0 + print*,'Excitation energies : CIS CIS_DT (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_dt(i) - eigenvalues_CIS_dress_D_dt(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.1,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + enddo + + convert = 627.51d0 + print*,'Excitation energies : CIS CIS_DT (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_dt(i) - eigenvalues_CIS_dress_D_dt(1) + write(*,'(I3,xxxxxxxxxxxxxxxx,5(F16.5,x))')i,delta_E_CIS*convert,delta_E_CIS_DT*convert + enddo + +!if(save_all_dm_cis)then +! call save_all_density_matrix +!endif +end diff --git a/src/CIS_dressed/CIS_providers.irp.f b/src/CIS_dressed/CIS_providers.irp.f index 9048aec3..f4fc6266 100644 --- a/src/CIS_dressed/CIS_providers.irp.f +++ b/src/CIS_dressed/CIS_providers.irp.f @@ -207,13 +207,13 @@ implicit none double precision,allocatable :: delta_H_matrix_doub(:,:) double precision,allocatable :: eigvalues(:),eigvectors(:,:) - double precision :: overlap,max_overlap,s2 + double precision :: overlap,max_overlap,s2,e_corr integer :: i_overlap,i,j,k 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 - delta_H_matrix_doub = 0.d0 + 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 @@ -246,7 +246,6 @@ enddo 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) diff --git a/src/CIS_dressed/EN2.irp.f b/src/CIS_dressed/EN2.irp.f index d89fbb9f..2b3af615 100644 --- a/src/CIS_dressed/EN2.irp.f +++ b/src/CIS_dressed/EN2.irp.f @@ -1,29 +1,33 @@ - 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 - + 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 @@ -31,222 +35,209 @@ 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 + if(EN_2_2)then + do i=n_core_cis+1,elec_alpha_num + + e_i=Fock_matrix_diag_mo(i) + + do k=elec_alpha_num+1,n_act_cis + + e_k=Fock_matrix_diag_mo(k) + delta_e_ik=e_i-e_k - 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 + !same spin contribution for EN2_corr_energy + do j=i+1,elec_alpha_num + e_j=Fock_matrix_diag_mo(j) + delta_e_ikj=delta_e_ik+e_j - 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 - + !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=Fock_matrix_diag_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=Fock_matrix_diag_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=Fock_matrix_diag_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=Fock_matrix_diag_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=Fock_matrix_diag_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 - !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 + else + do i=n_core_cis+1,elec_alpha_num + + e_i=Fock_matrix_diag_mo(i) + + do k=elec_alpha_num+1,n_act_cis + + e_k=Fock_matrix_diag_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=Fock_matrix_diag_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=Fock_matrix_diag_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/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=Fock_matrix_diag_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=Fock_matrix_diag_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/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=Fock_matrix_diag_mo(j) + delta_e_ikj=delta_e_ik+e_j + + do l=elec_beta_num+1,n_act_cis + e_l=Fock_matrix_diag_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/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 - 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 - - + 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 72012fde..0d9a8bb0 100644 --- a/src/CIS_dressed/MP2.irp.f +++ b/src/CIS_dressed/MP2.irp.f @@ -29,22 +29,22 @@ do i=n_core_cis+1,elec_alpha_num h_imp(i)=0.d0 - e_i=diagonal_Fock_matrix_mo(i) + e_i=Fock_matrix_diag_mo(i) do k=elec_alpha_num+1,n_act_cis hp_imp(i,k)=0.d0 - e_k=diagonal_Fock_matrix_mo(k) + e_k=Fock_matrix_diag_mo(k) delta_e_ik=e_i-e_k !same spin contribution for MP2_corr_energy do j=i+1,elec_alpha_num - e_j=diagonal_Fock_matrix_mo(j) + e_j=Fock_matrix_diag_mo(j) delta_e_ikj=delta_e_ik+e_j !same spin contribution for MP2_corr_energy and a part of the contribution to p_imp and h_imp do l=k+1,n_act_cis - e_l=diagonal_Fock_matrix_mo(l) + e_l=Fock_matrix_diag_mo(l) delta_e=delta_e_ikj-e_l direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) @@ -60,7 +60,7 @@ !remaining same spin contribution for p_imp do l=elec_alpha_num+1,k-1 - e_l=diagonal_Fock_matrix_mo(l) + e_l=Fock_matrix_diag_mo(l) delta_e=delta_e_ikj-e_l direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) @@ -75,11 +75,11 @@ !remaining same spin contribution for h_imp do j=n_core_cis+1,i-1 - e_j=diagonal_Fock_matrix_mo(j) + e_j=Fock_matrix_diag_mo(j) delta_e_ikj=delta_e_ik+e_j do l=k+1,n_act_cis - e_l=diagonal_Fock_matrix_mo(l) + e_l=Fock_matrix_diag_mo(l) delta_e=delta_e_ikj-e_l direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) @@ -94,11 +94,11 @@ !same spin contribution for hp_imp do j=n_core_cis+1,elec_alpha_num - e_j=diagonal_Fock_matrix_mo(j) + e_j=Fock_matrix_diag_mo(j) delta_e_ikj=delta_e_ik+e_j do l=elec_alpha_num+1,n_act_cis - e_l=diagonal_Fock_matrix_mo(l) + e_l=Fock_matrix_diag_mo(l) delta_e=delta_e_ikj-e_l direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) @@ -113,11 +113,11 @@ !different spin contribution do j=n_core_cis+1,elec_beta_num - e_j=diagonal_Fock_matrix_mo(j) + e_j=Fock_matrix_diag_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) + e_l=Fock_matrix_diag_mo(l) delta_e=delta_e_ikj-e_l direct=get_mo_bielec_integral(i,j,k,l,mo_integrals_map) @@ -385,25 +385,20 @@ 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 -! 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) + 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) -! 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 @@ -809,24 +804,24 @@ delta_H_trip=0.d0 !do i=5,6 do i=n_core_cis+1,elec_alpha_num - e_i=diagonal_Fock_matrix_mo(i) + e_i=Fock_matrix_diag_mo(i) !r=7 do r=elec_alpha_num+1,n_act_cis - e_r=diagonal_Fock_matrix_mo(r) + e_r=Fock_matrix_diag_mo(r) delta_e_ir=e_i-e_r key_ir=psi_CIS_adress(i,r) do j=i+1,elec_alpha_num - e_j=diagonal_Fock_matrix_mo(j) + e_j=Fock_matrix_diag_mo(j) delta_e_irj=delta_e_ir+e_j do s=r+1,n_act_cis - e_s=diagonal_Fock_matrix_mo(s) + e_s=Fock_matrix_diag_mo(s) delta_e_irjs=delta_e_irj-e_s !alpha-alpha-alpha do k=j+1,elec_alpha_num - e_k=diagonal_Fock_matrix_mo(k) + e_k=Fock_matrix_diag_mo(k) delta_e_irjsk=delta_e_irjs+e_k do t=s+1,n_act_cis - e_t=diagonal_Fock_matrix_mo(t) + e_t=Fock_matrix_diag_mo(t) delta_e=delta_e_irjsk-e_t energy=1.d0/(ref_energy-delta_e) occ=i @@ -898,10 +893,10 @@ enddo !alpha-alpha-beta do k=n_core_cis+1,elec_alpha_num - e_k=diagonal_Fock_matrix_mo(k) + e_k=Fock_matrix_diag_mo(k) delta_e_irjsk=delta_e_irjs+e_k do t=elec_alpha_num+1,n_act_cis - e_t=diagonal_Fock_matrix_mo(t) + e_t=Fock_matrix_diag_mo(t) delta_e=delta_e_irjsk-e_t energy=1.d0/(ref_energy-delta_e) occ=i @@ -992,23 +987,23 @@ !!generating the Singles included in the Triples !! do m=1,n_state_CIS !do i=n_core_cis+1,elec_alpha_num -! e_i=diagonal_Fock_matrix_mo(i) +! e_i=Fock_matrix_diag_mo(i) !!do r=elec_alpha_num+1,n_state_CIS ! do r=elec_alpha_num+1,n_act_cis -! e_r=diagonal_Fock_matrix_mo(r) +! e_r=Fock_matrix_diag_mo(r) ! delta_e_ir=e_i-e_r !key_ir=psi_CIS_adress(i,r) ! !alpha-alpha-x (=beta-beta-x) ! do j=i+1,elec_alpha_num -! e_j=diagonal_Fock_matrix_mo(j) +! e_j=Fock_matrix_diag_mo(j) ! delta_e_irj=delta_e_ir+e_j !key_jr=psi_CIS_adress(j,r) ! do s=r+1,n_act_cis -! e_s=diagonal_Fock_matrix_mo(s) +! e_s=Fock_matrix_diag_mo(s) ! delta_e_irjs=delta_e_irj-e_s ! key_is=psi_CIS_adress(i,s) @@ -1017,13 +1012,13 @@ ! !alpha-alpha-alpha (=beta-beta-beta) ! do k=j+1,elec_alpha_num -! e_k=diagonal_Fock_matrix_mo(k) +! e_k=Fock_matrix_diag_mo(k) ! delta_e_irjsk=delta_e_irjs+e_k ! key_kr=psi_CIS_adress(k,r) ! key_ks=psi_CIS_adress(k,s) ! do t=s+1,n_act_cis -! e_t=diagonal_Fock_matrix_mo(t) +! e_t=Fock_matrix_diag_mo(t) ! delta_e=delta_e_irjsk-e_t ! key_it=psi_CIS_adress(i,t) @@ -1157,7 +1152,7 @@ ! !alpha-alpha-beta (=beta-beta-alpha) ! do k=n_core_cis+1,elec_beta_num -! e_k=diagonal_Fock_matrix_mo(k) +! e_k=Fock_matrix_diag_mo(k) ! delta_e_irjsk=delta_e_irjs+e_k ! key_kr=psi_CIS_adress(k,r) @@ -1166,7 +1161,7 @@ ! do t=elec_beta_num+1,n_act_cis -! e_t=diagonal_Fock_matrix_mo(t) +! e_t=Fock_matrix_diag_mo(t) ! delta_e=delta_e_irjsk-e_t ! key_it=psi_CIS_adress(i,t) @@ -1298,53 +1293,3 @@ !enddo END - -subroutine diexcitation(i,j,k,l,ispin1,ispin2,key_in,key_out,i_ok,Nint) - implicit none - use bitmasks - ! realize the double excitation i-->k (ispin1) + j-->l (ispin2) on key_in - ! returns key_out and i_ok (i_ok = 0 means not possible, i_ok = 1 means the excitation was possible) - integer, intent(in) :: ispin1,ispin2,i,j,k,l,Nint - integer(bit_kind), intent(in) :: key_in(Nint,2) - integer, intent(out):: i_ok - integer(bit_kind), intent(out):: key_out(Nint,2) - integer :: k_hole,j_hole,k_particl,j_particl,i_nint,Nelec_alpha,Nelec_beta - integer :: i_test_hole,i_test_particl - key_out = key_in - - k_hole = ishft(i-1,-bit_kind_shift)+1 - j_hole = i-ishft(k_hole-1,bit_kind_shift)-1 - i_test_hole = ibset(0,j_hole) - if(iand(key_in(k_hole,ispin1),i_test_hole).ne.i_test_hole)then - i_ok = 0 - return - endif - key_out(k_hole,ispin1) = ibclr(key_out(k_hole,ispin1),j_hole) - k_particl = ishft(k-1,-bit_kind_shift)+1 - j_particl = k-ishft(k_particl-1,bit_kind_shift)-1 - i_test_particl= ibset(0,j_particl) - if(iand(key_in(k_particl,ispin1),i_test_particl).ne.0)then - i_ok = 0 - return - endif - key_out(k_particl,ispin1) = ibset(key_out(k_particl,ispin1),j_particl) - - k_hole = ishft(j-1,-bit_kind_shift)+1 - j_hole = j-ishft(k_hole-1,bit_kind_shift)-1 - i_test_hole = ibset(0,j_hole) - if(iand(key_in(k_hole,ispin2),i_test_hole).ne.i_test_hole)then - i_ok = 0 - return - endif - key_out(k_hole,ispin2) = ibclr(key_out(k_hole,ispin2),j_hole) - k_particl = ishft(l-1,-bit_kind_shift)+1 - j_particl = l-ishft(k_particl-1,bit_kind_shift)-1 - i_test_particl = ibset(0,j_particl) - if(iand(key_in(k_particl,ispin2),i_test_particl).ne.0)then - i_ok = 0 - return - endif - key_out(k_particl,ispin2) = ibset(key_out(k_particl,ispin2),j_particl) - i_ok = 1 - end - diff --git a/src/CIS_dressed/natural_particl_orbitals.irp.f b/src/CIS_dressed/natural_particl_orbitals.irp.f index 532f640f..ba6d4554 100644 --- a/src/CIS_dressed/natural_particl_orbitals.irp.f +++ b/src/CIS_dressed/natural_particl_orbitals.irp.f @@ -56,7 +56,7 @@ enddo do k = 1, mo_tot_num do l = 1, mo_tot_num - c_k = eigvectors(k,j) * eigvectors(l,j) + c_k = eigvectors(k,i) * eigvectors(l,i) particle_natural_orb_CIS_properties(1,i) += c_k * mo_dipole_x(k,l) particle_natural_orb_CIS_properties(2,i) += c_k * mo_dipole_y(k,l) particle_natural_orb_CIS_properties(3,i) += c_k * mo_dipole_z(k,l) diff --git a/src/CIS_dressed/options.irp.f b/src/CIS_dressed/options.irp.f index 15643737..2dc27eb9 100644 --- a/src/CIS_dressed/options.irp.f +++ b/src/CIS_dressed/options.irp.f @@ -10,7 +10,8 @@ BEGIN_PROVIDER [ integer , n_state_cis ] if (has) then call ezfio_get_cis_dressed_n_state_cis(n_state_cis) else - n_state_cis = 5 + n_state_cis = 10 + call ezfio_set_cis_dressed_n_state_cis(n_state_cis) endif END_PROVIDER diff --git a/src/CIS_dressed/repeat_all_doubles.irp.f b/src/CIS_dressed/repeat_all_doubles.irp.f new file mode 100644 index 00000000..21b2dabb --- /dev/null +++ b/src/CIS_dressed/repeat_all_doubles.irp.f @@ -0,0 +1,178 @@ +subroutine repeat_all_doubles(key_in, e_corr) + implicit none + integer(bit_kind), intent(in) :: key_in(N_int,2) + double precision, intent(out) :: e_corr + integer :: i,j,k,l + integer :: s1,s2 + integer :: i_ok + double precision :: hij,get_mo_bielec_integral,diag_H_mat_elem,delta_e + integer(bit_kind) :: key_out(N_int,2) + ! same spin (alpha) + if(mp2_dressing)then + s1 = 1 + e_corr = 0.d0 + do i = n_core_cis + 1, elec_alpha_num + do j = i + 1, elec_alpha_num + do k = elec_alpha_num + 1, n_act_cis + do l = k+1, n_act_cis + ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in> + call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int) + if(i_ok.ne.0)then + hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & + -get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + e_corr += hij*hij/(Fock_matrix_diag_mo(i)+Fock_matrix_diag_mo(j) & + -Fock_matrix_diag_mo(l)-Fock_matrix_diag_mo(k)) + endif + enddo + enddo + enddo + enddo + + s1 = 2 + do i = n_core_cis + 1, elec_alpha_num + do j = i + 1, elec_alpha_num + do k = elec_alpha_num + 1, n_act_cis + do l = k+1, n_act_cis + ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in> + call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int) + if(i_ok.ne.0)then + hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & + -get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + e_corr += hij*hij/(Fock_matrix_diag_mo(i)+Fock_matrix_diag_mo(j) & + -Fock_matrix_diag_mo(l)-Fock_matrix_diag_mo(k)) + endif + enddo + enddo + enddo + enddo + + s1 = 1 + s2 = 2 + do i = n_core_cis + 1, elec_alpha_num + do j = n_core_cis + 1, elec_alpha_num + do k = elec_alpha_num + 1, n_act_cis + do l = elec_alpha_num + 1, n_act_cis + ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in> + call diexcitation(i,j,k,l,s1,s2,key_in,key_out,i_ok,N_int) + if(i_ok.ne.0)then + hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + e_corr += hij*hij/(Fock_matrix_diag_mo(i)+Fock_matrix_diag_mo(j) & + -Fock_matrix_diag_mo(l)-Fock_matrix_diag_mo(k)) + endif + enddo + enddo + enddo + enddo + + else + ! same spin (alpha) + s1 = 1 + e_corr = 0.d0 + do i = n_core_cis + 1, elec_alpha_num + do j = i + 1, elec_alpha_num + do k = elec_alpha_num + 1, n_act_cis + do l = k+1, n_act_cis + ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in> + call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int) + if(i_ok.ne.0)then + hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & + -get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + call diexcitation(i,j,k,l,s1,s1,ref_bitmask,key_out,i_ok,N_int) + delta_e = HF_energy - diag_H_mat_elem(key_out,N_int) + e_corr += hij*hij/delta_e + endif + enddo + enddo + enddo + enddo + + s1 = 2 + do i = n_core_cis + 1, elec_alpha_num + do j = i + 1, elec_alpha_num + do k = elec_alpha_num + 1, n_act_cis + do l = k+1, n_act_cis + ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in> + call diexcitation(i,j,k,l,s1,s1,key_in,key_out,i_ok,N_int) + if(i_ok.ne.0)then + hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) & + -get_mo_bielec_integral(i,j,l,k,mo_integrals_map) + call diexcitation(i,j,k,l,s1,s1,ref_bitmask,key_out,i_ok,N_int) + delta_e = HF_energy - diag_H_mat_elem(key_out,N_int) + e_corr += hij*hij/delta_e + endif + enddo + enddo + enddo + enddo + + s1 = 1 + s2 = 2 + do i = n_core_cis + 1, elec_alpha_num + do j = n_core_cis + 1, elec_alpha_num + do k = elec_alpha_num + 1, n_act_cis + do l = elec_alpha_num + 1, n_act_cis + ! a^+ k(s1) a^+ l(s1) a_j(s1) a_i(s1) |key_in> + call diexcitation(i,j,k,l,s1,s2,key_in,key_out,i_ok,N_int) + if(i_ok.ne.0)then + hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) + call diexcitation(i,j,k,l,s1,s2,ref_bitmask,key_out,i_ok,N_int) + delta_e = HF_energy - diag_H_mat_elem(key_out,N_int) + e_corr += hij*hij/delta_e + endif + enddo + enddo + enddo + enddo + endif +end + +subroutine diexcitation(i,j,k,l,ispin1,ispin2,key_in,key_out,i_ok,Nint) + implicit none + use bitmasks + ! realize the double excitation i-->k (ispin1) + j-->l (ispin2) on key_in + ! returns key_out and i_ok (i_ok = 0 means not possible, i_ok = 1 means the excitation was possible) + integer, intent(in) :: ispin1,ispin2,i,j,k,l,Nint + integer(bit_kind), intent(in) :: key_in(Nint,2) + integer, intent(out):: i_ok + integer(bit_kind), intent(out):: key_out(Nint,2) + integer :: k_hole,j_hole,k_particl,j_particl,i_nint,Nelec_alpha,Nelec_beta + integer(bit_kind) :: i_test_hole,i_test_particl + character*(512) :: output(2) + key_out = key_in + + k_hole = ishft(i-1,-bit_kind_shift)+1 + j_hole = i-ishft(k_hole-1,bit_kind_shift)-1 + i_test_hole = ibset(0_bit_kind,j_hole) + if(iand(key_in(k_hole,ispin1),i_test_hole).ne.i_test_hole)then + i_ok = 0 + return + endif + key_out(k_hole,ispin1) = ibclr(key_out(k_hole,ispin1),j_hole) + k_particl = ishft(k-1,-bit_kind_shift)+1 + j_particl = k-ishft(k_particl-1,bit_kind_shift)-1 + i_test_particl= ibset(0_bit_kind,j_particl) + if(iand(key_in(k_particl,ispin1),i_test_particl).ne.0_bit_kind)then + i_ok = 0 + return + endif + key_out(k_particl,ispin1) = ibset(key_out(k_particl,ispin1),j_particl) + + k_hole = ishft(j-1,-bit_kind_shift)+1 + j_hole = j-ishft(k_hole-1,bit_kind_shift)-1 + i_test_hole = ibset(0_bit_kind,j_hole) + if(iand(key_in(k_hole,ispin2),i_test_hole).ne.i_test_hole)then + i_ok = 0 + return + endif + key_out(k_hole,ispin2) = ibclr(key_out(k_hole,ispin2),j_hole) + k_particl = ishft(l-1,-bit_kind_shift)+1 + j_particl = l-ishft(k_particl-1,bit_kind_shift)-1 + i_test_particl = ibset(0_bit_kind,j_particl) + if(iand(key_in(k_particl,ispin2),i_test_particl).ne.0_bit_kind)then + i_ok = 0 + return + endif + key_out(k_particl,ispin2) = ibset(key_out(k_particl,ispin2),j_particl) + i_ok = 1 + end + diff --git a/src/Full_CI/parameters.irp.f b/src/Full_CI/parameters.irp.f new file mode 100644 index 00000000..4d70df93 --- /dev/null +++ b/src/Full_CI/parameters.irp.f @@ -0,0 +1,51 @@ +BEGIN_PROVIDER [ integer, N_det_max_fci ] + implicit none + BEGIN_DOC +! Max number od determinants in the wave function + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_full_ci_n_det_max_fci(exists) + if (exists) then + call ezfio_get_full_ci_n_det_max_fci(n_det_max_fci) + else + n_det_max_fci = 10000 + endif + call write_int(output_full_ci,n_det_max_fci,'Max number of determinants ') + ASSERT (n_det_max_fci > 0) +END_PROVIDER + + +BEGIN_PROVIDER [ logical , do_pt2_end ] + implicit none + BEGIN_DOC +! if True then compute the PT2 when the selection process is finished + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_full_ci_do_pt2_end(exists) + if (exists) then + call ezfio_get_full_ci_do_pt2_end(do_pt2_end) + else + do_pt2_end = .True. + endif +!call write_i(output_full_ci,do_pt2_end,' computes the PT2 at the end of the selection ') + ASSERT (do_pt2_end > 0) +END_PROVIDER + + +BEGIN_PROVIDER [ double precision , pt2_max ] + implicit none + BEGIN_DOC +! The selection process stops when the largest PT2 (for all the states) is lower than pt2_max +! in absolute value + END_DOC + logical :: exists + PROVIDE ezfio_filename + call ezfio_has_full_ci_pt2_max(exists) + if (exists) then + call ezfio_get_full_ci_pt2_max(pt2_max) + else + pt2_max = 0.1d0 + endif +END_PROVIDER diff --git a/src/Hartree_Fock/Fock_matrix.irp.f b/src/Hartree_Fock/Fock_matrix.irp.f index 4ff99248..70eb1afa 100644 --- a/src/Hartree_Fock/Fock_matrix.irp.f +++ b/src/Hartree_Fock/Fock_matrix.irp.f @@ -115,6 +115,7 @@ END_PROVIDER if (do_direct_integrals) then PROVIDE all_utils ao_overlap_abs ao_integrals_threshold PROVIDE HF_density_matrix_ao_alpha HF_density_matrix_ao_beta + PROVIDE ao_bi_elec_integral_alpha !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral) & !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,& @@ -147,6 +148,8 @@ END_PROVIDER !$OMP END PARALLEL else + PROVIDE ao_bielec_integrals_in_map + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ao_ints_val,ao_ints_idx,kmax) & !$OMP SHARED(ao_num,ao_bi_elec_integral_alpha,ao_mono_elec_integral,&