diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 6e23ebac..e0d0e02e 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 6e23ebac001acae91d1c762ca934e09a9b7d614a +Subproject commit e0d0e02e9f5ece138d1520106954a881ab0b8db2 diff --git a/src/cipsi_tc_bi_ortho/pt2.irp.f b/src/cipsi_tc_bi_ortho/pt2.irp.f index 13b4dff4..833cc0ea 100644 --- a/src/cipsi_tc_bi_ortho/pt2.irp.f +++ b/src/cipsi_tc_bi_ortho/pt2.irp.f @@ -1,4 +1,4 @@ -subroutine pt2_tc_bi_ortho +subroutine tc_pt2 use selection_types implicit none BEGIN_DOC @@ -15,7 +15,7 @@ subroutine pt2_tc_bi_ortho double precision, external :: memory_of_double double precision :: correlation_energy_ratio,E_denom,E_tc,norm double precision, allocatable :: ept2(:), pt1(:),extrap_energy(:) - PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map + PROVIDE H_apply_buffer_allocated distributed_davidson print*,'Diagonal elements of the Fock matrix ' do i = 1, mo_num @@ -44,24 +44,14 @@ subroutine pt2_tc_bi_ortho pt2_data % overlap= 0.d0 pt2_data % variance = huge(1.e0) - if (s2_eig) then - call make_s2_eigenfunction - endif + !!!! WARNING !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION +! if (s2_eig) then +! call make_s2_eigenfunction +! endif print_pt2 = .False. call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) ! call routine_save_right - if (N_det > N_det_max) then - psi_det(1:N_int,1:2,1:N_det) = psi_det_generators(1:N_int,1:2,1:N_det) - psi_coef(1:N_det,1:N_states) = psi_coef_sorted_tc_gen(1:N_det,1:N_states) - N_det = N_det_max - soft_touch N_det psi_det psi_coef - if (s2_eig) then - call make_s2_eigenfunction - endif - print_pt2 = .False. - call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) - endif allocate(ept2(1000),pt1(1000),extrap_energy(100)) @@ -71,18 +61,11 @@ subroutine pt2_tc_bi_ortho ! soft_touch thresh_it_dav print_pt2 = .True. - to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor) - to_select = max(N_states_diag, to_select) - - E_denom = E_tc ! TC Energy of the current wave function call pt2_dealloc(pt2_data) call pt2_dealloc(pt2_data_err) call pt2_alloc(pt2_data, N_states) call pt2_alloc(pt2_data_err, N_states) - call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection - - N_iter += 1 - + call ZMQ_pt2(E_tc, pt2_data, pt2_data_err, relative_error,0) ! Stochastic PT2 and selection call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) end diff --git a/src/fci_tc_bi/pt2_tc.irp.f b/src/fci_tc_bi/pt2_tc.irp.f new file mode 100644 index 00000000..96a54825 --- /dev/null +++ b/src/fci_tc_bi/pt2_tc.irp.f @@ -0,0 +1,31 @@ +program tc_pt2_prog + implicit none + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + pruning = -1.d0 + touch pruning +! pt2_relative_error = 0.01d0 +! touch pt2_relative_error + call run_pt2_tc + +end + + +subroutine run_pt2_tc + + implicit none + + PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e + if(elec_alpha_num+elec_beta_num.ge.3)then + if(three_body_h_tc)then + call provide_all_three_ints_bi_ortho + endif + endif + ! --- + + call tc_pt2 + + +end diff --git a/src/kohn_sham/print_mos.irp.f b/src/kohn_sham/print_mos.irp.f index 5e728444..19bb98bc 100644 --- a/src/kohn_sham/print_mos.irp.f +++ b/src/kohn_sham/print_mos.irp.f @@ -3,7 +3,7 @@ program print_mos integer :: i,nx double precision :: r(3), xmax, dx, accu double precision, allocatable :: mos_array(:) - double precision:: alpha,envelop + double precision:: alpha,envelop,dm_a,dm_b allocate(mos_array(mo_num)) xmax = 5.d0 nx = 1000 @@ -11,20 +11,14 @@ program print_mos r = 0.d0 alpha = 0.5d0 do i = 1, nx + call dm_dft_alpha_beta_at_r(r,dm_a,dm_b) call give_all_mos_at_r(r,mos_array) accu = mos_array(3)**2+mos_array(4)**2+mos_array(5)**2 accu = dsqrt(accu) envelop = (1.d0 - dexp(-alpha * r(3)**2)) - write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, envelop + write(33,'(100(F16.10,X))')r(3), mos_array(1), mos_array(2), accu, dm_a+dm_b, envelop r(3) += dx enddo end -double precision function f_mu(x) - implicit none - double precision, intent(in) :: x - - - -end diff --git a/src/tc_bi_ortho/e_corr_bi_ortho.irp.f b/src/tc_bi_ortho/e_corr_bi_ortho.irp.f index ec66a8b5..3a715b44 100644 --- a/src/tc_bi_ortho/e_corr_bi_ortho.irp.f +++ b/src/tc_bi_ortho/e_corr_bi_ortho.irp.f @@ -45,6 +45,9 @@ &BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj ] &BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth ] &BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth ] +&BEGIN_PROVIDER [ double precision, e_corr_bi_orth_proj_abs ] +&BEGIN_PROVIDER [ double precision, e_corr_single_bi_orth_abs ] +&BEGIN_PROVIDER [ double precision, e_corr_double_bi_orth_abs ] implicit none integer :: i,degree double precision :: hmono,htwoe,hthree,htilde_ij @@ -57,13 +60,15 @@ call htilde_mu_mat_bi_ortho(HF_bitmask,psi_det(1,1,i),N_int,hmono,htwoe,hthree,htilde_ij) if(degree == 1)then e_corr_single_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) + e_corr_single_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)) else if(degree == 2)then e_corr_double_bi_orth += reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) -! print*,'coef_wf , e_cor',reigvec_tc_bi_orth(i,1)/reigvec_tc_bi_orth(1,1), reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1) + e_corr_double_bi_orth_abs += dabs(reigvec_tc_bi_orth(i,1) * htilde_ij/reigvec_tc_bi_orth(1,1)) endif enddo e_corr_bi_orth_proj = e_corr_single_bi_orth + e_corr_double_bi_orth e_corr_bi_orth = eigval_right_tc_bi_orth(1) - e_tilde_bi_orth_00 + e_corr_bi_orth_proj_abs = e_corr_single_bi_orth_abs + e_corr_double_bi_orth_abs END_PROVIDER BEGIN_PROVIDER [ double precision, e_tc_left_right ] diff --git a/src/tc_bi_ortho/pt2_tc_cisd.irp.f b/src/tc_bi_ortho/pt2_tc_cisd.irp.f new file mode 100644 index 00000000..ecf5bb42 --- /dev/null +++ b/src/tc_bi_ortho/pt2_tc_cisd.irp.f @@ -0,0 +1,43 @@ +program pt2_tc_cisd + + BEGIN_DOC + ! + ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together + ! with the energy. Saves the left-right wave functions at the end. + ! + END_DOC + + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + read_wf = .True. + touch read_wf + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + print*, ' nb of states = ', N_states + print*, ' nb of det = ', N_det + + call routine +end + +subroutine routine + implicit none + integer :: i,h1,p1,h2,p2,s1,s2 + double precision :: h0i,hi0,e00,ei,delta_e + double precision :: norm,e_corr,coef + norm = 0.d0 + e_corr = 0.d0 + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,1), psi_det(1,1,1), N_int, e00) + do i = 2, N_det + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,1), N_int, hi0) + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,1), psi_det(1,1,i), N_int, h0i) + call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,i), psi_det(1,1,i), N_int, ei) + delta_e = e00 - ei + coef = hi0/delta_e + norm += coef*coef + e_corr += dabs(coef* h0i) + enddo + print*,'e_corr = ',e_corr + print*,'norm = ',norm + +end diff --git a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f index 00cebf3a..5a3f9935 100644 --- a/src/tc_bi_ortho/slater_tc_opt_diag.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_diag.irp.f @@ -156,7 +156,7 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo - if(three_body_h_tc)then + if(three_body_h_tc.and.elec_num.gt.2)then !!!!! 3-e part !! same-spin/same-spin do j = 1, na @@ -243,7 +243,7 @@ subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb) htwoe= htwoe- mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb) enddo - if(three_body_h_tc)then + if(three_body_h_tc.and.elec_num.gt.2)then !!!!! 3-e part !! same-spin/same-spin do j = 1, na diff --git a/src/tc_bi_ortho/slater_tc_opt_double.irp.f b/src/tc_bi_ortho/slater_tc_opt_double.irp.f index baca498c..1b0e43bb 100644 --- a/src/tc_bi_ortho/slater_tc_opt_double.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_double.irp.f @@ -41,15 +41,15 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, if(s1.ne.s2)then ! opposite spin two-body htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) - if(three_body_h_tc)then + if(three_body_h_tc.and.elec_num.gt.2)then if(.not.double_normal_ord)then if(degree_i>degree_j)then call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree) else call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) endif - elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then - htwoe += normal_two_body_bi_orth(p2,h2,p1,h1)!!! WTF ??? + elseif(double_normal_ord.and.elec_num.gt.2)then + htwoe += normal_two_body_bi_orth(p2,h2,p1,h1) endif endif else @@ -58,16 +58,16 @@ subroutine double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, htwoe = mo_bi_ortho_tc_two_e(p2,p1,h2,h1) ! exchange terms htwoe -= mo_bi_ortho_tc_two_e(p1,p2,h2,h1) - if(three_body_h_tc)then + if(three_body_h_tc.and.elec_num.gt.2)then if(.not.double_normal_ord)then if(degree_i>degree_j)then call three_comp_two_e_elem(key_j,h1,h2,p1,p2,s1,s2,hthree) else call three_comp_two_e_elem(key_i,h1,h2,p1,p2,s1,s2,hthree) endif - elseif(double_normal_ord.and.elec_num+elec_num.gt.2)then - htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2)!!! WTF ??? - htwoe += normal_two_body_bi_orth(h1,p1,h2,p2)!!! WTF ??? + elseif(double_normal_ord.and.elec_num.gt.2)then + htwoe -= normal_two_body_bi_orth(h2,p1,h1,p2) + htwoe += normal_two_body_bi_orth(h1,p1,h2,p2) endif endif endif diff --git a/src/tc_bi_ortho/slater_tc_opt_single.irp.f b/src/tc_bi_ortho/slater_tc_opt_single.irp.f index 7cff3c73..2f9d83bf 100644 --- a/src/tc_bi_ortho/slater_tc_opt_single.irp.f +++ b/src/tc_bi_ortho/slater_tc_opt_single.irp.f @@ -106,7 +106,7 @@ subroutine get_single_excitation_from_fock_tc(key_i,key_j,h,p,spin,phase,hmono,h htwoe -= buffer_x(i) enddo hthree = 0.d0 - if (three_body_h_tc)then + if (three_body_h_tc.and.elec_num.gt.2)then call three_comp_fock_elem(key_i,h,p,spin,hthree) endif diff --git a/src/tc_bi_ortho/tc_bi_ortho.irp.f b/src/tc_bi_ortho/tc_bi_ortho.irp.f index b69a2fe6..f69684c0 100644 --- a/src/tc_bi_ortho/tc_bi_ortho.irp.f +++ b/src/tc_bi_ortho/tc_bi_ortho.irp.f @@ -48,17 +48,20 @@ subroutine routine_diag() if(N_states .eq. 1) then - print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) - print*,'e_tc_left_right = ',e_tc_left_right - print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00 - print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth - print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single - print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double - print*,'***' - print*,'e_corr_bi_orth = ',e_corr_bi_orth - print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj - print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth - print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth + print*,'eigval_right_tc_bi_orth = ',eigval_right_tc_bi_orth(1) + print*,'e_tc_left_right = ',e_tc_left_right + print*,'e_tilde_bi_orth_00 = ',e_tilde_bi_orth_00 + print*,'e_pt2_tc_bi_orth = ',e_pt2_tc_bi_orth + print*,'e_pt2_tc_bi_orth_single = ',e_pt2_tc_bi_orth_single + print*,'e_pt2_tc_bi_orth_double = ',e_pt2_tc_bi_orth_double + print*,'***' + print*,'e_corr_bi_orth = ',e_corr_bi_orth + print*,'e_corr_bi_orth_proj = ',e_corr_bi_orth_proj + print*,'e_corr_bi_orth_proj_abs = ',e_corr_bi_orth_proj_abs + print*,'e_corr_single_bi_orth = ',e_corr_single_bi_orth + print*,'e_corr_double_bi_orth = ',e_corr_double_bi_orth + print*,'e_corr_single_bi_orth_abs = ',e_corr_single_bi_orth_abs + print*,'e_corr_double_bi_orth_abs = ',e_corr_double_bi_orth_abs print*,'Left/right eigenvectors' do i = 1,N_det write(*,'(I5,X,(100(F12.7,X)))')i,leigvec_tc_bi_orth(i,1),reigvec_tc_bi_orth(i,1),leigvec_tc_bi_orth(i,1)*reigvec_tc_bi_orth(i,1)