From 3b2735886eb8f6e668d21ffcaa5a37d99b933a19 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 1 Apr 2016 00:05:44 +0200 Subject: [PATCH 1/5] Test in lambda --- plugins/MRCC_Utils/mrcc_utils.irp.f | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 80fdd4c7..25e6af77 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -21,7 +21,8 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] endif lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) ) lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - if (lambda_pert / lambda_mrcc(k,i) < 0.5d0 ) then + if ((lambda_pert < 0.5d0 * lambda_mrcc(k,i)).or. & + (lambda_pert > 2.0d0 * lambda_mrcc(k,i)) ) then i_pert_count += 1 lambda_mrcc(k,i) = 0.d0 endif From dd8647d5aa2e2596987c76596bb124ee7dcbfd77 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 1 Apr 2016 01:31:11 +0200 Subject: [PATCH 2/5] Introduced provider for psi_ref_coef_normalized --- plugins/MRCC_CASSD/mrcc_noiter.irp.f | 3 ++- plugins/MRCC_Utils/mrcc_utils.irp.f | 5 ++--- plugins/Psiref_CAS/psi_ref.irp.f | 16 +++++++++++++++- 3 files changed, 19 insertions(+), 5 deletions(-) diff --git a/plugins/MRCC_CASSD/mrcc_noiter.irp.f b/plugins/MRCC_CASSD/mrcc_noiter.irp.f index 949a0dab..d36bedef 100644 --- a/plugins/MRCC_CASSD/mrcc_noiter.irp.f +++ b/plugins/MRCC_CASSD/mrcc_noiter.irp.f @@ -1,7 +1,8 @@ program mrcc_noiter implicit none read_wf = .True. - SOFT_TOUCH read_wf + threshold_generators = .9999d0 + SOFT_TOUCH read_wf threshold_generators call print_cas_coefs call set_generators_bitmasks_as_holes_and_particles call run diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index 25e6af77..d83dc239 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] lambda_mrcc = 0.d0 do i=1,N_det_non_ref - call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,& + call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_normalized, N_int, N_det_ref,& size(psi_ref_coef,1), N_states,ihpsi_current) call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii) do k=1,N_states @@ -21,8 +21,7 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] endif lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) ) lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii) - if ((lambda_pert < 0.5d0 * lambda_mrcc(k,i)).or. & - (lambda_pert > 2.0d0 * lambda_mrcc(k,i)) ) then + if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then i_pert_count += 1 lambda_mrcc(k,i) = 0.d0 endif diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index d3dbda08..c7b9867f 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -23,7 +23,21 @@ use bitmasks psi_ref_coef(i,k) = psi_cas_coef(i,k) enddo enddo - call normalize(psi_ref_coef, N_det_ref) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_ref_coef_normalized, (psi_det_size,n_states) ] + implicit none + BEGIN_DOC +! Normalized coefficients of the reference + END_DOC + integer :: i,j,k + do k=1,N_states + do j=1,N_det_ref + psi_ref_coef_normalized(j,k) = psi_ref_coef(j,k) + enddo + call normalize(psi_ref_coef_normalized(1,k), N_det_ref) + enddo END_PROVIDER From 27b7c341066514c5c8701b2a4cf57af2da28b0f3 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 1 Apr 2016 23:33:58 +0200 Subject: [PATCH 3/5] MRCC+PT2 --- plugins/MRCC_CASSD/mrcc_cassd.irp.f | 58 +++++++++++++++++++- plugins/MRCC_CASSD/mrcc_noiter.irp.f | 55 ++++++++++++++++++- plugins/MRCC_Utils/davidson.irp.f | 5 +- plugins/MRCC_Utils/mrcc_utils.irp.f | 15 ++++- plugins/Perturbation/perturbation.template.f | 14 +++-- plugins/Perturbation/pt2_equations.irp.f | 38 ++++++------- scripts/generate_h_apply.py | 9 +-- 7 files changed, 155 insertions(+), 39 deletions(-) diff --git a/plugins/MRCC_CASSD/mrcc_cassd.irp.f b/plugins/MRCC_CASSD/mrcc_cassd.irp.f index c5ed175f..38cd3c55 100644 --- a/plugins/MRCC_CASSD/mrcc_cassd.irp.f +++ b/plugins/MRCC_CASSD/mrcc_cassd.irp.f @@ -1,15 +1,25 @@ program mrcc implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) + read_wf = .True. SOFT_TOUCH read_wf call print_cas_coefs call set_generators_bitmasks_as_holes_and_particles - call run + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) end -subroutine run +subroutine run(N_st,energy) implicit none + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) + integer :: i double precision :: E_new, E_old, delta_e @@ -37,10 +47,54 @@ subroutine run endif enddo call write_double(6,ci_energy_dressed(1),"Final MRCC energy") + energy(:) = ci_energy_dressed(:) end +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + pt2 = 0.d0 + + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + + N_det_generators = lambda_mrcc_pt2(0) + do i=1,N_det_generators + j = lambda_mrcc_pt2(i) + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + + + call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', energy + print *, 'E+PT2 = ', energy+pt2 + print *, '-----' + + + call ezfio_set_full_ci_energy_pt2(energy+pt2) + deallocate(pt2,norm_pert) + +end + + subroutine print_cas_coefs implicit none diff --git a/plugins/MRCC_CASSD/mrcc_noiter.irp.f b/plugins/MRCC_CASSD/mrcc_noiter.irp.f index d36bedef..8d95cea9 100644 --- a/plugins/MRCC_CASSD/mrcc_noiter.irp.f +++ b/plugins/MRCC_CASSD/mrcc_noiter.irp.f @@ -1,17 +1,24 @@ program mrcc_noiter implicit none + double precision, allocatable :: energy(:) + allocate (energy(N_states)) read_wf = .True. threshold_generators = .9999d0 SOFT_TOUCH read_wf threshold_generators call print_cas_coefs call set_generators_bitmasks_as_holes_and_particles - call run + call run(N_states,energy) + if(do_pt2_end)then + call run_pt2(N_states,energy) + endif + deallocate(energy) end -subroutine run +subroutine run(N_st,energy) implicit none - double precision :: lambda + integer, intent(in) :: N_st + double precision, intent(out) :: energy(N_st) integer :: i,j do j=1,N_states_diag do i=1,N_det @@ -22,6 +29,48 @@ subroutine run call write_double(6,ci_energy_dressed(1),"Final MRCC energy") call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1)) call save_wavefunction + energy(:) = ci_energy_dressed(:) +end + + +subroutine run_pt2(N_st,energy) + implicit none + integer :: i,j,k + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer, intent(in) :: N_st + double precision, intent(in) :: energy(N_st) + allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) + pt2 = 0.d0 + + print*,'Last iteration only to compute the PT2' + threshold_selectors = 1.d0 + threshold_generators = 0.999d0 + + N_det_generators = lambda_mrcc_pt2(0) + do i=1,N_det_generators + j = lambda_mrcc_pt2(i) + do k=1,N_int + psi_det_generators(k,1,i) = psi_non_ref(k,1,j) + psi_det_generators(k,2,i) = psi_non_ref(k,2,j) + enddo + do k=1,N_st + psi_coef_generators(i,k) = psi_non_ref_coef(j,k) + enddo + enddo + SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed + + + call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st) + print *, 'Final step' + print *, 'N_det = ', N_det + print *, 'N_states = ', N_states + print *, 'PT2 = ', pt2 + print *, 'E = ', energy + print *, 'E+PT2 = ', energy+pt2 + print *, '-----' + + call ezfio_set_full_ci_energy_pt2(energy+pt2) + deallocate(pt2,norm_pert) end diff --git a/plugins/MRCC_Utils/davidson.irp.f b/plugins/MRCC_Utils/davidson.irp.f index d278ba13..354ea389 100644 --- a/plugins/MRCC_Utils/davidson.irp.f +++ b/plugins/MRCC_Utils/davidson.irp.f @@ -189,10 +189,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin ! Davidson iterations ! =================== - converged = .False. + integer :: iteration + converged = .False. do while (.not.converged) - !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st) do k=1,N_st @@ -206,6 +206,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin do iter=1,davidson_sze_max-1 + ! Compute W_k = H |u_k> ! ---------------------- diff --git a/plugins/MRCC_Utils/mrcc_utils.irp.f b/plugins/MRCC_Utils/mrcc_utils.irp.f index d83dc239..e6c84a06 100644 --- a/plugins/MRCC_Utils/mrcc_utils.irp.f +++ b/plugins/MRCC_Utils/mrcc_utils.irp.f @@ -1,15 +1,19 @@ -BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] + BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] +&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ] implicit none BEGIN_DOC ! cm/ or perturbative 1/Delta_E(m) END_DOC integer :: i,k double precision :: ihpsi_current(N_states) - integer :: i_pert_count - double precision :: hii, lambda_pert + integer :: i_pert_count + double precision :: hii, lambda_pert + integer :: N_lambda_mrcc_pt2 i_pert_count = 0 lambda_mrcc = 0.d0 + N_lambda_mrcc_pt2 = 0 + lambda_mrcc_pt2(0) = 0 do i=1,N_det_non_ref call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_normalized, N_int, N_det_ref,& @@ -24,9 +28,14 @@ BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ] if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then i_pert_count += 1 lambda_mrcc(k,i) = 0.d0 + if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then + N_lambda_mrcc_pt2 += 1 + lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i + endif endif enddo enddo + lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2 print*,'N_det_non_ref = ',N_det_non_ref print*,'Number of ignored determinants = ',i_pert_count diff --git a/plugins/Perturbation/perturbation.template.f b/plugins/Perturbation/perturbation.template.f index 94b6b8b0..7bb08c21 100644 --- a/plugins/Perturbation/perturbation.template.f +++ b/plugins/Perturbation/perturbation.template.f @@ -3,7 +3,7 @@ import perturbation END_SHELL -subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp) +subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -14,6 +14,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) integer(bit_kind),intent(in) :: key_mask(Nint,2) double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num) + double precision, intent(in) :: electronic_energy(N_st) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) @@ -151,7 +152,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c idx_microlist_zero(ptr_microlist(1)+l) = idx_microlist(ptr_microlist(smallerlist)+l) enddo end if - call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_microlist(smallerlist)+N_microlist(0), & n_st,microlist_zero,idx_microlist_zero,N_microlist(smallerlist)+N_microlist(0)) else @@ -160,11 +161,11 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c cycle end if - call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) end if -! call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & +! call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & ! c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st @@ -182,7 +183,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c end -subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp) +subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy) implicit none BEGIN_DOC ! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply @@ -193,6 +194,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size) integer(bit_kind),intent(in) :: key_mask(Nint,2) double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num) + double precision, intent(in) :: electronic_energy(N_st) double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st) double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st) double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st) @@ -241,7 +243,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_ cycle endif - call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & + call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, & c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist) do k = 1,N_st diff --git a/plugins/Perturbation/pt2_equations.irp.f b/plugins/Perturbation/pt2_equations.irp.f index 72d03808..e990a37c 100644 --- a/plugins/Perturbation/pt2_equations.irp.f +++ b/plugins/Perturbation/pt2_equations.irp.f @@ -29,11 +29,11 @@ subroutine pt2_epstein_nesbet ($arguments) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st - if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then + if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then c_pert(i) = -1.d0 e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h) + else if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h) H_pert_diag(i) = h*c_pert(i)*c_pert(i) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) else @@ -66,7 +66,6 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) double precision :: i_H_psi_array(N_st) ASSERT (Nint == N_int) ASSERT (Nint > 0) - PROVIDE CI_electronic_energy !call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array) call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array) @@ -74,7 +73,7 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st if (i_H_psi_array(i) /= 0.d0) then - delta_e = h - CI_electronic_energy(i) + delta_e = h - electronic_energy(i) if (delta_e > 0.d0) then e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i))) else @@ -165,7 +164,7 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments) ! ! that can be repeated by repeating all the double excitations ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! : you repeat all the correlation energy already taken into account in electronic_energy(1) ! ! that could be repeated to this determinant. ! @@ -195,16 +194,16 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments) enddo h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) h = h + accu_e_corr - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + delta_E = 1.d0/(electronic_energy(1) - h) - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h) e_2_pert(1) = i_H_psi_array(1) * c_pert(1) do i =2,N_st H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h)) e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) else c_pert(i) = i_H_psi_array(i) @@ -248,7 +247,7 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments) ! ! that can be repeated by repeating all the double excitations ! - ! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1) + ! : you repeat all the correlation energy already taken into account in electronic_energy(1) ! ! that could be repeated to this determinant. ! @@ -278,16 +277,16 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments) enddo h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) h = h + accu_e_corr - delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h) + delta_E = 1.d0/(electronic_energy(1) - h) - c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h) + c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h) e_2_pert(1) = i_H_psi_array(1) * c_pert(1) do i =2,N_st H_pert_diag(i) = h - if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h)) + if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h)) e_2_pert(i) = (c_pert(i) * i_H_psi_array(i)) else c_pert(i) = i_H_psi_array(i) @@ -328,11 +327,11 @@ subroutine pt2_epstein_nesbet_sc2 ($arguments) h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) do i =1,N_st - if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then + if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then c_pert(i) = -1.d0 e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0 - else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then - c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h) + else if (dabs(electronic_energy(i) - h) > 1.d-6) then + c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h) H_pert_diag(i) = h*c_pert(i)*c_pert(i) e_2_pert(i) = c_pert(i) * i_H_psi_array(i) else @@ -348,7 +347,7 @@ end SUBST [ arguments, declarations ] -det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ; +electronic_energy,det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ; integer, intent(in) :: Nint integer, intent(in) :: ndet @@ -357,6 +356,7 @@ det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minili integer(bit_kind), intent(in) :: det_ref (Nint,2) integer(bit_kind), intent(in) :: det_pert(Nint,2) double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1) + double precision , intent(in) :: electronic_energy(N_st) double precision , intent(out) :: c_pert(N_st) double precision , intent(out) :: e_2_pert(N_st) double precision, intent(out) :: H_pert_diag(N_st) diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index c6466569..436f092d 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -61,6 +61,7 @@ class H_apply(object): s["params_post"] = "" self.selection_pt2 = None + self.energy = "CI_electronic_energy" self.perturbation = None self.do_double_exc = do_double_exc #s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) & @@ -264,13 +265,13 @@ class H_apply(object): self.data["keys_work"] = """ ! if(check_double_excitation)then call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & - sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) - """%(pert) + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s) + """%(pert,self.energy) else: self.data["keys_work"] = """ call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, & - sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp) - """%(pert) + sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s) + """%(pert,self.energy) self.data["finalization"] = """ From d20c63f0fe3e14996ae719e455b30470e0fc6a18 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Fri, 1 Apr 2016 23:46:24 +0200 Subject: [PATCH 4/5] Forgot H_apply --- plugins/MRCC_Utils/H_apply.irp.f | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index decc5a75..df2b67a0 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -24,5 +24,12 @@ s.data["size_max"] = "3072" print s + +s = H_apply_zmq("mrcc_PT2") +s.energy = "ci_electronic_energy_dressed" +s.set_perturbation("epstein_nesbet_2x2") +s.unset_openmp() +print s + END_SHELL From f0fc7dd39f5b8b7df9e706318d79cb2496506f09 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Sat, 2 Apr 2016 00:10:27 +0200 Subject: [PATCH 5/5] Fixed tests --- tests/bats/qp.bats | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/bats/qp.bats b/tests/bats/qp.bats index 6f1c916b..1ced9e1d 100644 --- a/tests/bats/qp.bats +++ b/tests/bats/qp.bats @@ -155,7 +155,7 @@ function run_all_1h_1p() { ezfio set determinants read_wf True qp_run mrcc_cassd $INPUT energy="$(ezfio get mrcc_cassd energy)" - eq $energy -76.2286679037553 1.e-4 + eq $energy -76.2284994316618 1.e-4 }