From 9121d1a60476bef94a10497ee399e0342a8d597c Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 8 Sep 2016 15:48:52 +0200 Subject: [PATCH 1/3] corrected bugs in fock for MRPT --- plugins/MRPT_Utils/fock_like_operators.irp.f | 10 +++---- plugins/MRPT_Utils/mrpt_utils.irp.f | 29 ++++++++++++++------ 2 files changed, 24 insertions(+), 15 deletions(-) diff --git a/plugins/MRPT_Utils/fock_like_operators.irp.f b/plugins/MRPT_Utils/fock_like_operators.irp.f index 5900516e..2074daf6 100644 --- a/plugins/MRPT_Utils/fock_like_operators.irp.f +++ b/plugins/MRPT_Utils/fock_like_operators.irp.f @@ -11,8 +11,6 @@ accu = 0.d0 i_inact_core_orb = list_core_inact(i) do j = 1, n_core_inact_orb -! do j = 1, elec_alpha_num -! j_inact_core_orb = j j_inact_core_orb = list_core_inact(j) accu += 2.d0 * mo_bielec_integral_jj(i_inact_core_orb,j_inact_core_orb) & - mo_bielec_integral_jj_exchange(i_inact_core_orb,j_inact_core_orb) @@ -84,8 +82,8 @@ accu_exchange(2) += 2.d0 * nb * exchange enddo enddo - fock_core_inactive_from_act(i_inact_core_orb,1,i_state) = accu_coulomb + accu_exchange(1) - fock_core_inactive_from_act(i_inact_core_orb,2,i_state) = accu_coulomb + accu_exchange(2) + fock_core_inactive_from_act(i_inact_core_orb,1,i_state) = accu_coulomb - accu_exchange(1) + fock_core_inactive_from_act(i_inact_core_orb,2,i_state) = accu_coulomb - accu_exchange(2) enddo enddo END_PROVIDER @@ -131,8 +129,8 @@ accu_exchange(2) += 2.d0 * nb * exchange enddo enddo - fock_virt_from_act(i_virt_orb,1,i_state) = accu_coulomb + accu_exchange(1) - fock_virt_from_act(i_virt_orb,2,i_state) = accu_coulomb + accu_exchange(2) + fock_virt_from_act(i_virt_orb,1,i_state) = accu_coulomb - accu_exchange(1) + fock_virt_from_act(i_virt_orb,2,i_state) = accu_coulomb - accu_exchange(2) enddo enddo END_PROVIDER diff --git a/plugins/MRPT_Utils/mrpt_utils.irp.f b/plugins/MRPT_Utils/mrpt_utils.irp.f index 3b832b58..e298ae67 100644 --- a/plugins/MRPT_Utils/mrpt_utils.irp.f +++ b/plugins/MRPT_Utils/mrpt_utils.irp.f @@ -130,19 +130,30 @@ print*, '2h1p = ',accu ! 2h2p - delta_ij_tmp = 0.d0 - call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) - accu = 0.d0 +!delta_ij_tmp = 0.d0 +!call H_apply_mrpt_2h2p(delta_ij_tmp,N_det) +!accu = 0.d0 +!do i_state = 1, N_states +!do i = 1, N_det +! do j = 1, N_det +! accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) +! delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) +! enddo +!enddo +!second_order_pt_new_2h2p(i_state) = accu(i_state) +!enddo +!print*, '2h2p = ',accu + + double precision :: contrib_2h2p(N_states) + call give_2h2p(contrib_2h2p) do i_state = 1, N_states do i = 1, N_det - do j = 1, N_det - accu(i_state) += delta_ij_tmp(j,i,i_state) * psi_coef(i,i_state) * psi_coef(j,i_state) - delta_ij(j,i,i_state) += delta_ij_tmp(j,i,i_state) - enddo + delta_ij(i,i,i_state) += contrib_2h2p(i_state) enddo - second_order_pt_new_2h2p(i_state) = accu(i_state) + second_order_pt_new_2h2p(i_state) = contrib_2h2p(i_state) enddo - print*, '2h2p = ',accu + print*, '2h2p = ',contrib_2h2p(1) + ! total accu = 0.d0 From 9a152ca037c246284c494f5c619e154124a7ce82 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 10 Sep 2016 12:32:33 +0200 Subject: [PATCH 2/3] Beginning new way for computing pt2 --- config/ifort.cfg | 2 +- plugins/MRPT_Utils/excitations_cas.irp.f | 2 -- plugins/MRPT_Utils/mrpt_dress.irp.f | 15 +++++++++++++++ plugins/MRPT_Utils/psi_active_prov.irp.f | 6 ++++-- src/Determinants/slater_rules.irp.f | 12 ++++++++++++ 5 files changed, 32 insertions(+), 5 deletions(-) diff --git a/config/ifort.cfg b/config/ifort.cfg index da414912..a738a83c 100644 --- a/config/ifort.cfg +++ b/config/ifort.cfg @@ -18,7 +18,7 @@ IRPF90_FLAGS : --ninja --align=32 # 0 : Deactivate # [OPTION] -MODE : OPT ; [ OPT | PROFILE | DEBUG ] : Chooses the section below +MODE : DEBUG ; [ OPT | PROFILE | DEBUG ] : Chooses the section below CACHE : 1 ; Enable cache_compile.py OPENMP : 1 ; Append OpenMP flags diff --git a/plugins/MRPT_Utils/excitations_cas.irp.f b/plugins/MRPT_Utils/excitations_cas.irp.f index 5024967d..03571d5d 100644 --- a/plugins/MRPT_Utils/excitations_cas.irp.f +++ b/plugins/MRPT_Utils/excitations_cas.irp.f @@ -398,8 +398,6 @@ subroutine u0_H_dyall_u0(energies,psi_in,psi_in_coef,ndet,dim_psi_in,dim_psi_coe do j = 1, ndet if(psi_coef_tmp(j)==0.d0)cycle call i_H_j_dyall(psi_in(1,1,i),psi_in(1,1,j),N_int,hij) -! call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij_bis) -! print*, hij_bis,hij accu += psi_coef_tmp(i) * psi_coef_tmp(j) * hij enddo enddo diff --git a/plugins/MRPT_Utils/mrpt_dress.irp.f b/plugins/MRPT_Utils/mrpt_dress.irp.f index 512158cf..0c8ec98c 100644 --- a/plugins/MRPT_Utils/mrpt_dress.irp.f +++ b/plugins/MRPT_Utils/mrpt_dress.irp.f @@ -44,6 +44,8 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip integer :: N_miniList, leng double precision :: delta_e(N_states),hij_tmp integer :: index_i,index_j + double precision :: phase_array(N_det),phase + integer :: exc(0:2,2,2),degree leng = max(N_det_generators, N_det) @@ -74,11 +76,14 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip ! double precision :: ihpsi0,coef_pert ! ihpsi0 = 0.d0 ! coef_pert = 0.d0 + phase_array =0.d0 do i = 1,idx_alpha(0) index_i = idx_alpha(i) call get_delta_e_dyall(psi_det(1,1,index_i),tq(1,1,i_alpha),delta_e) call i_h_j(tq(1,1,i_alpha),psi_det(1,1,index_i),Nint,hialpha) hij_array(index_i) = hialpha + call get_excitation(psi_det(1,1,index_i),tq(1,1,i_alpha),exc,degree,phase,N_int) +! phase_array(index_i) = phase do i_state = 1,N_states delta_e_inv_array(index_i,i_state) = 1.d0/delta_e(i_state) enddo @@ -90,6 +95,16 @@ subroutine mrpt_dress(delta_ij_, Ndet,i_generator,n_selected,det_buffer,Nint,ip call omp_set_lock( psi_ref_bis_lock(index_i) ) do j = 1, idx_alpha(0) index_j = idx_alpha(j) +! call get_excitation(psi_det(1,1,index_i),psi_det(1,1,index_i),exc,degree,phase,N_int) +! if(index_j.ne.index_i)then +! if(phase_array(index_j) * phase_array(index_i) .ne. phase)then +! print*, phase_array(index_j) , phase_array(index_i) ,phase +! call debug_det(psi_det(1,1,index_i),N_int) +! call debug_det(psi_det(1,1,index_j),N_int) +! call debug_det(tq(1,1,i_alpha),N_int) +! stop +! endif +! endif do i_state=1,N_states ! standard dressing first order delta_ij_(index_i,index_j,i_state) += hij_array(index_j) * hij_tmp * delta_e_inv_array(index_j,i_state) diff --git a/plugins/MRPT_Utils/psi_active_prov.irp.f b/plugins/MRPT_Utils/psi_active_prov.irp.f index 8d705deb..6fb8219e 100644 --- a/plugins/MRPT_Utils/psi_active_prov.irp.f +++ b/plugins/MRPT_Utils/psi_active_prov.irp.f @@ -180,6 +180,8 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) double precision :: delta_e_inactive(N_states) integer :: i_hole_inact + + call give_holes_in_inactive_space(det_2,n_holes_spin,n_holes,holes_list) delta_e_inactive = 0.d0 do i = 1, n_holes_spin(1) @@ -429,7 +431,7 @@ subroutine get_delta_e_dyall(det_1,det_2,delta_e_final) do i_state = 1, n_states - delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state) - enddo + delta_e_final(i_state) = delta_e_act(i_state) + delta_e_inactive(i_state) - delta_e_virt(i_state) + enddo end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index 735ca8e7..f70fa594 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -749,6 +749,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) exc(1,2,2) ,mo_integrals_map) else if (exc(0,1,1) == 2) then ! Double alpha + print*,'phase hij = ',phase hij = phase*(get_mo_bielec_integral_schwartz( & exc(1,1,1), & exc(2,1,1), & @@ -759,6 +760,17 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) exc(2,1,1), & exc(2,2,1), & exc(1,2,1) ,mo_integrals_map) ) + print*,get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) + print*,get_mo_bielec_integral_schwartz( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) + else if (exc(0,1,2) == 2) then ! Double beta hij = phase*(get_mo_bielec_integral_schwartz( & From 156cbbdb27cb910f5fc0f9c5cb7ea1bfa79ac004 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Sat, 10 Sep 2016 16:51:09 +0200 Subject: [PATCH 3/3] New way pt2 is ok for 2h1p --- plugins/MRPT_Utils/new_way.irp.f | 192 ++++++++++++++++++++++++++++ src/Determinants/slater_rules.irp.f | 91 +++++++++++++ 2 files changed, 283 insertions(+) create mode 100644 plugins/MRPT_Utils/new_way.irp.f diff --git a/plugins/MRPT_Utils/new_way.irp.f b/plugins/MRPT_Utils/new_way.irp.f new file mode 100644 index 00000000..83831e41 --- /dev/null +++ b/plugins/MRPT_Utils/new_way.irp.f @@ -0,0 +1,192 @@ +subroutine give_2h1p_contrib(matrix_2h1p) + use bitmasks + implicit none + double precision , intent(inout) :: matrix_2h1p(N_det,N_det,*) + integer :: i,j,r,a,b + integer :: iorb, jorb, rorb, aorb, borb + integer :: ispin,jspin + integer :: idet,jdet + integer(bit_kind) :: perturb_dets(N_int,2,n_act_orb,2,2) + double precision :: perturb_dets_phase(n_act_orb,2,2) + double precision :: perturb_dets_hij(n_act_orb,2,2) + double precision :: coef_perturb_from_idet(n_act_orb,2,2,N_states) + integer :: inint + integer :: elec_num_tab_local(2),acu_elec + integer(bit_kind) :: det_tmp(N_int,2) + integer :: exc(0:2,2,2) + integer :: accu_elec + double precision :: get_mo_bielec_integral_schwartz + double precision :: active_int(n_act_orb,2) + double precision :: hij,phase +!matrix_2h1p = 0.d0 + + elec_num_tab_local = 0 + do inint = 1, N_int + elec_num_tab_local(1) += popcnt(psi_det(inint,1,1)) + elec_num_tab_local(2) += popcnt(psi_det(inint,2,1)) + enddo + do i = 1, n_inact_orb ! First inactive + iorb = list_inact(i) + do j = 1, n_inact_orb ! Second inactive + jorb = list_inact(j) + do r = 1, n_virt_orb ! First virtual + rorb = list_virt(r) + ! take all the integral you will need for i,j,r fixed + do a = 1, n_act_orb + aorb = list_act(a) + active_int(a,1) = get_mo_bielec_integral_schwartz(iorb,jorb,rorb,aorb,mo_integrals_map) ! direct + active_int(a,2) = get_mo_bielec_integral_schwartz(iorb,jorb,aorb,rorb,mo_integrals_map) ! exchange + enddo + + integer :: degree(N_det) + integer :: idx(0:N_det) + double precision :: delta_e(n_act_orb,2,N_states) + integer :: istate + integer :: index_orb_act_mono(N_det,3) + + do idet = 1, N_det + call get_excitation_degree_vector_mono(psi_det,psi_det(1,1,idet),degree,N_int,N_det,idx) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Precomputation of matrix elements + do ispin = 1, 2 ! spin of the couple a-a^dagger (i,r) + do jspin = 1, 2 ! spin of the couple z-a^dagger (j,a) + if(ispin == jspin .and. iorb.le.jorb)cycle ! condition not to double count + do a = 1, n_act_orb ! First active + aorb = list_act(a) + do inint = 1, N_int + det_tmp(inint,1) = psi_det(inint,1,idet) + det_tmp(inint,2) = psi_det(inint,2,idet) + enddo + ! Do the excitation inactive -- > virtual + call clear_bit_to_integer(iorb,det_tmp(1,ispin),N_int) ! hole in "iorb" of spin Ispin + call set_bit_to_integer(rorb,det_tmp(1,ispin),N_int) ! particle in "rorb" of spin Ispin + + ! Do the excitation inactive -- > active + call clear_bit_to_integer(jorb,det_tmp(1,jspin),N_int) ! hole in "jorb" of spin Jspin + call set_bit_to_integer(aorb,det_tmp(1,jspin),N_int) ! particle in "aorb" of spin Jspin + + ! Check if the excitation is possible or not on psi_det(idet) + accu_elec= 0 + do inint = 1, N_int + accu_elec+= popcnt(det_tmp(inint,jspin)) + enddo + if(accu_elec .ne. elec_num_tab_local(jspin))then + perturb_dets_phase(a,jspin,ispin) = 0.0 + perturb_dets_hij(a,jspin,ispin) = 0.d0 + do istate = 1, N_states + coef_perturb_from_idet(a,jspin,ispin,istate) = 0.d0 + enddo + cycle + endif + do inint = 1, N_int + perturb_dets(inint,1,a,jspin,ispin) = det_tmp(inint,1) + perturb_dets(inint,2,a,jspin,ispin) = det_tmp(inint,2) + enddo + call get_double_excitation(psi_det(1,1,idet),det_tmp,exc,phase,N_int) + perturb_dets_phase(a,jspin,ispin) = phase + do istate = 1, N_states + delta_e(a,jspin,istate) = one_creat(a,jspin,istate) & + - fock_virt_total_spin_trace(rorb,istate) & + + fock_core_inactive_total_spin_trace(iorb,istate) & + + fock_core_inactive_total_spin_trace(jorb,istate) + enddo + if(ispin == jspin)then + perturb_dets_hij(a,jspin,ispin) = phase * (active_int(a,2) - active_int(a,1) ) + else + perturb_dets_hij(a,jspin,ispin) = phase * active_int(a,1) + endif +!!!!!!!!!!!!!!!!!!!!!1 Computation of the coefficient at first order coming from idet +!!!!!!!!!!!!!!!!!!!!! for the excitation (i,j)(ispin,jspin) ---> (r,a)(ispin,jspin) + do istate = 1, N_states + coef_perturb_from_idet(a,jspin,ispin,istate) = perturb_dets_hij(a,jspin,ispin) / delta_e(a,jspin,istate) + enddo + + enddo + enddo + enddo + +!!!!!!!!!!!!!!!!!!!!!!!!!!! determination of the connections between I and the other J determinants mono excited in the CAS +!!!!!!!!!!!!!!!!!!!!!!!!!!!! the determinants I and J must be connected by the following operator +!!!!!!!!!!!!!!!!!!!!!!!!!!!! + do jdet = 1, idx(0) + if(idx(jdet).ne.idet)then + call get_mono_excitation(psi_det(1,1,idet),psi_det(1,1,idx(jdet)),exc,phase,N_int) + if (exc(0,1,1) == 1) then + ! Mono alpha + index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,1)) !!! a^{\dagger}_a + index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,1)) !!! a_{b} + index_orb_act_mono(idx(jdet),3) = 1 + else + ! Mono beta + index_orb_act_mono(idx(jdet),1) = list_act_reverse(exc(1,2,2)) !!! a^{\dagger}_a + index_orb_act_mono(idx(jdet),2) = list_act_reverse(exc(1,1,2)) !!! a_{b} + index_orb_act_mono(idx(jdet),3) = 2 + endif + else + index_orb_act_mono(idx(jdet),1) = -1 + endif + enddo + + integer :: kspin + do jdet = 1, idx(0) + if(idx(jdet).ne.idet)then + ! two determinants | Idet > and | Jdet > which are connected throw a mono excitation operator + ! are connected by the presence of the perturbers determinants |det_tmp> + aorb = index_orb_act_mono(idx(jdet),1) ! a^{\dagger}_{aorb} + borb = index_orb_act_mono(idx(jdet),2) ! a_{borb} + kspin = index_orb_act_mono(idx(jdet),3) ! spin of the excitation + ! the determinants Idet and Jdet interact throw the following operator + ! | Jdet > = a_{borb,kspin} a^{\dagger}_{aorb, kspin} | Idet > + + do ispin = 1, 2 ! you loop on all possible spin for the excitation + ! a^{\dagger}_r a_{i} (ispin) + if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count + + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet > + do inint = 1, N_int + det_tmp(inint,1) = perturb_dets(inint,1,aorb,kspin,ispin) + det_tmp(inint,2) = perturb_dets(inint,2,aorb,kspin,ispin) + enddo + double precision :: hja + ! you determine the interaction between the excited determinant and the other parent | Jdet > + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{borb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Jdet > + ! hja = < det_tmp | H | Jdet > + call get_double_excitation(psi_det(1,1,idx(jdet)),det_tmp,exc,phase,N_int) + if(kspin == ispin)then + hja = phase * (active_int(borb,2) - active_int(borb,1) ) + else + hja = phase * active_int(borb,1) + endif + + do istate = 1, N_states + matrix_2h1p(idx(jdet),idet,istate) += hja * coef_perturb_from_idet(aorb,kspin,ispin,istate) + enddo + enddo ! ispin + + else + ! diagonal part of the dressing : interaction of | Idet > with all the perturbers generated by the excitations + ! + ! | det_tmp > = a^{\dagger}_{rorb,ispin} a^{\dagger}_{aorb,kspin} a_{jorb,kspin} a_{iorb,ispin} | Idet > + do ispin = 1, 2 + do kspin = 1, 2 + if(ispin == kspin .and. iorb.le.jorb)cycle ! condition not to double count + do a = 1, n_act_orb ! First active + do istate = 1, N_states + matrix_2h1p(idet,idet,istate) += coef_perturb_from_idet(a,kspin,ispin,istate) * perturb_dets_hij(a,kspin,ispin) + enddo + enddo + enddo + enddo + + endif + + enddo + enddo + enddo + enddo + enddo + + + + + +end diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index f70fa594..0357fb88 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -1237,6 +1237,97 @@ subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_a print*,'------' end +subroutine get_excitation_degree_vector_mono(key1,key2,degree,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Applies get_excitation_degree to an array of determinants and return only the mono excitations + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: degree(sze) + integer, intent(out) :: idx(0:sze) + + integer :: i,l,d,m + + ASSERT (Nint > 0) + ASSERT (sze > 0) + + l=1 + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (d > 2) then + cycle + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (d > 2) then + cycle + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (d > 2) then + cycle + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + d = 0 + !DIR$ LOOP COUNT MIN(4) + do m=1,Nint + d = d + popcnt(xor( key1(m,1,i), key2(m,1))) & + + popcnt(xor( key1(m,2,i), key2(m,2))) + enddo + if (d > 2) then + cycle + else + degree(l) = ishft(d,-1) + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)