From 6a6c197b9a47a6696ae4831b1d0a2458d1b5d08c Mon Sep 17 00:00:00 2001 From: Manu Date: Wed, 18 Mar 2015 11:35:16 +0100 Subject: [PATCH] updated NEEDED_MODULES in src repository --- scripts/generate_h_apply.py | 8 ++ src/BiInts/map_integrals.irp.f | 4 + src/Dets/H_apply_template.f | 5 + src/Dets/program_beginer_determinants.irp.f | 138 ++++++++++++++++++++ src/MonoInts/pot_mo_ints.irp.f | 6 +- src/NEEDED_MODULES | 2 +- 6 files changed, 159 insertions(+), 4 deletions(-) create mode 100644 src/Dets/program_beginer_determinants.irp.f diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 61dbf1b9..5fbf1c71 100755 --- a/scripts/generate_h_apply.py +++ b/scripts/generate_h_apply.py @@ -23,6 +23,7 @@ deinit_thread skip init_main filter_integrals +filter2h2p """.split() class H_apply(object): @@ -115,6 +116,13 @@ class H_apply(object): buffer = buffer.replace('$'+key, value) return buffer + def set_filter_2h_2p(self): + self["filter2h2p"] = """ +! ! DIR$ FORCEINLINE + if(is_a_two_holes_two_particles(key))cycle + """ + + def set_perturbation(self,pert): if self.perturbation is not None: raise diff --git a/src/BiInts/map_integrals.irp.f b/src/BiInts/map_integrals.irp.f index ff0a83f6..82d3b49b 100644 --- a/src/BiInts/map_integrals.irp.f +++ b/src/BiInts/map_integrals.irp.f @@ -276,6 +276,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map) implicit none BEGIN_DOC ! Returns one integral in the MO basis + ! i(1)j(1) 1/r12 k(2)l(2) END_DOC integer, intent(in) :: i,j,k,l integer*8 :: idx @@ -292,6 +293,7 @@ double precision function mo_bielec_integral(i,j,k,l) implicit none BEGIN_DOC ! Returns one integral in the MO basis + ! i(1)j(1) 1/r12 k(2)l(2) END_DOC integer, intent(in) :: i,j,k,l double precision :: get_mo_bielec_integral @@ -306,6 +308,7 @@ subroutine get_mo_bielec_integrals(j,k,l,sze,out_val,map) BEGIN_DOC ! Returns multiple integrals in the MO basis, all ! i for j,k,l fixed. + ! i(1)j(1) 1/r12 k(2)l(2) END_DOC integer, intent(in) :: j,k,l, sze real(integral_kind), intent(out) :: out_val(sze) @@ -327,6 +330,7 @@ subroutine get_mo_bielec_integrals_existing_ik(j,l,sze,out_array,map) implicit none BEGIN_DOC ! Returns multiple integrals in the MO basis, all + ! i(1)j(1) 1/r12 k(2)l(2) ! i for j,k,l fixed. END_DOC integer, intent(in) :: j,l, sze diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index 90959070..557a6325 100644 --- a/src/Dets/H_apply_template.f +++ b/src/Dets/H_apply_template.f @@ -26,6 +26,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) double precision :: mo_bielec_integral + logical :: is_a_two_holes_two_particles integer, allocatable :: ia_ja_pairs(:,:,:) integer, allocatable :: ib_jb_pairs(:,:) double precision :: diag_H_mat_elem @@ -162,6 +163,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,other_spin) = ibset(key(k,other_spin),l) + $filter2h2p key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -210,6 +212,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene k = ishft(j_b-1,-bit_kind_shift)+1 l = j_b-ishft(k-1,bit_kind_shift)-1 key(k,ispin) = ibset(key(k,ispin),l) + $filter2h2p key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = key(k,1) @@ -267,6 +270,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param integer :: kk,pp,other_spin,key_idx integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) + logical :: is_a_two_holes_two_particles integer, allocatable :: ia_ja_pairs(:,:,:) logical, allocatable :: array_pairs(:,:) @@ -333,6 +337,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc $param k_a = ishft(j_a-1,-bit_kind_shift)+1 l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) + $filter2h2p key_idx += 1 do k=1,N_int keys_out(k,1,key_idx) = hole(k,1) diff --git a/src/Dets/program_beginer_determinants.irp.f b/src/Dets/program_beginer_determinants.irp.f new file mode 100644 index 00000000..6375af22 --- /dev/null +++ b/src/Dets/program_beginer_determinants.irp.f @@ -0,0 +1,138 @@ +program pouet + implicit none + print*,'HF energy = ',ref_bitmask_energy + nuclear_repulsion + call routine + +end +subroutine routine + use bitmasks + implicit none + integer :: i,j,k,l + double precision :: hij,get_mo_bielec_integral + double precision :: hmono,h_bi_ispin,h_bi_other_spin + integer(bit_kind),allocatable :: key_tmp(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + ! First checks + print*,'N_int = ',N_int + print*,'mo_tot_num = ',mo_tot_num + print*,'mo_tot_num / 64+1= ',mo_tot_num/64+1 + ! We print the HF determinant + do i = 1, N_int + print*,'ref_bitmask(i,1) = ',ref_bitmask(i,1) + print*,'ref_bitmask(i,2) = ',ref_bitmask(i,2) + enddo + print*,'' + print*,'Hartree Fock determinant ...' + call debug_det(ref_bitmask,N_int) + allocate(key_tmp(N_int,2)) + ! We initialize key_tmp to the Hartree Fock one + key_tmp = ref_bitmask + integer :: i_hole,i_particle,ispin,i_ok,other_spin + ! We do a mono excitation on the top of the HF determinant + write(*,*)'Enter the (hole, particle) couple for the mono excitation ...' + read(5,*)i_hole,i_particle +!!i_hole = 4 +!!i_particle = 20 + write(*,*)'Enter the ispin variable ...' + write(*,*)'ispin = 1 ==> alpha ' + write(*,*)'ispin = 2 ==> beta ' + read(5,*)ispin + if(ispin == 1)then + other_spin = 2 + else if(ispin == 2)then + other_spin = 1 + else + print*,'PB !! ' + print*,'ispin must be 1 or 2 !' + stop + endif +!!ispin = 1 + call do_mono_excitation(key_tmp,i_hole,i_particle,ispin,i_ok) + ! We check if it the excitation was possible with "i_ok" + if(i_ok == -1)then + print*,'i_ok = ',i_ok + print*,'You can not do this excitation because of Pauli principle ...' + print*,'check your hole particle couple, there must be something wrong ...' + stop + + endif + print*,'New det = ' + call debug_det(key_tmp,N_int) + call i_H_j(key_tmp,ref_bitmask,N_int,hij) + ! We calculate the H matrix element between the new determinant and HF + print*,' = ',hij + print*,'' + print*,'' + print*,'Recalculating it old school style ....' + print*,'' + print*,'' + ! We recalculate this old school style !!! + ! Mono electronic part + hmono = mo_mono_elec_integral(i_hole,i_particle) + print*,'' + print*,'Mono electronic part ' + print*,'' + print*,' = ',hmono + h_bi_ispin = 0.d0 + h_bi_other_spin = 0.d0 + print*,'' + print*,'Getting all the info for the calculation of the bi electronic part ...' + print*,'' + allocate (occ(N_int*bit_kind_size,2)) + ! We get the occupation of the alpha electrons in occ(:,1) + call bitstring_to_list(key_tmp(1,1), occ(1,1), n_occ_alpha, N_int) + print*,'n_occ_alpha = ',n_occ_alpha + print*,'elec_alpha_num = ',elec_alpha_num + ! We get the occupation of the beta electrons in occ(:,2) + call bitstring_to_list(key_tmp(1,2), occ(1,2), n_occ_beta, N_int) + print*,'n_occ_beta = ',n_occ_beta + print*,'elec_beta_num = ',elec_beta_num + ! We print the occupation of the alpha electrons + print*,'Alpha electrons !' + do i = 1, n_occ_alpha + print*,'i = ',i + print*,'occ(i,1) = ',occ(i,1) + enddo + ! We print the occupation of the beta electrons + print*,'Alpha electrons !' + do i = 1, n_occ_beta + print*,'i = ',i + print*,'occ(i,2) = ',occ(i,2) + enddo + integer :: exc(0:2,2,2),degree,h1,p1,h2,p2,s1,s2 + double precision :: phase + + call get_excitation_degree(key_tmp,ref_bitmask,degree,N_int) + print*,'degree = ',degree + call get_mono_excitation(ref_bitmask,key_tmp,exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + print*,'h1 = ',h1 + print*,'p1 = ',p1 + print*,'s1 = ',s1 + print*,'phase = ',phase + do i = 1, elec_num_tab(ispin) + integer :: orb_occupied + orb_occupied = occ(i,ispin) + h_bi_ispin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) & + -get_mo_bielec_integral(i_hole,i_particle,orb_occupied,orb_occupied,mo_integrals_map) + enddo + print*,'h_bi_ispin = ',h_bi_ispin + + do i = 1, elec_num_tab(other_spin) + orb_occupied = occ(i,other_spin) + h_bi_other_spin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) + enddo + print*,'h_bi_other_spin = ',h_bi_other_spin + print*,'h_bi_ispin + h_bi_other_spin = ',h_bi_ispin + h_bi_other_spin + + print*,'Total matrix element = ',phase*(h_bi_ispin + h_bi_other_spin + hmono) +!i = 1 +!j = 1 +!k = 1 +!l = 1 +!hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) +!print*,' = ',hij + + +end diff --git a/src/MonoInts/pot_mo_ints.irp.f b/src/MonoInts/pot_mo_ints.irp.f index fab0583f..6393e57e 100644 --- a/src/MonoInts/pot_mo_ints.irp.f +++ b/src/MonoInts/pot_mo_ints.irp.f @@ -11,11 +11,11 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_to do i = 1, mo_tot_num do j = 1, mo_tot_num do i1 = 1,ao_num - c_i1 = mo_coef(i1,i) + c_i1 = mo_coef(i1,i) ! do j1 = 1,ao_num - c_j1 = c_i1*mo_coef(j1,j) + c_j1 = c_i1*mo_coef(j1,j) ! mo_nucl_elec_integral(j,i) = mo_nucl_elec_integral(j,i) + & - c_j1 * ao_nucl_elec_integral(j1,i1) + c_j1 * ao_nucl_elec_integral(j1,i1) enddo enddo enddo diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index a0699a91..e8752243 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -1 +1 @@ -AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS +AOs BiInts Bitmask Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS Denstity_stuff CAS_SD_selected DDCI_selected Two_body_density_matrix