From 92f321e5941ea7e6441450f6b57cfb80ecea8874 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 24 Mar 2020 09:54:48 -0500 Subject: [PATCH] ci kpts --- src/determinants/single_excitations.irp.f | 175 +++++++++++++-------- src/mo_two_e_ints/integrals_3_index.irp.f | 37 +++++ src/mo_two_e_ints/map_integrals_cplx.irp.f | 38 +++++ 3 files changed, 181 insertions(+), 69 deletions(-) diff --git a/src/determinants/single_excitations.irp.f b/src/determinants/single_excitations.irp.f index d9ca40f0..6e394572 100644 --- a/src/determinants/single_excitations.irp.f +++ b/src/determinants/single_excitations.irp.f @@ -1,7 +1,7 @@ use bitmasks BEGIN_PROVIDER [integer(bit_kind), ref_closed_shell_bitmask, (N_int,2)] implicit none - integer :: i,i0 + integer :: i,i0,k integer :: n_occ_ab(2) integer :: occ(N_int*bit_kind_size,2) call bitstring_to_list_ab(ref_bitmask, occ, n_occ_ab, N_int) @@ -11,18 +11,22 @@ BEGIN_PROVIDER [integer(bit_kind), ref_closed_shell_bitmask, (N_int,2)] ref_closed_shell_bitmask(i,2) = ref_bitmask(i,2) enddo if (is_complex) then - do + !todo: check this + do k=1,kpt_num + call bitstring_to_list_ab(ref_bitmask_kpts(1,1,k),occ,n_occ_ab,N_int) + do i0=elec_beta_num_kpts(k)+1,elec_alpha_num_kpts(k) + i=occ(i0,1) + call clear_bit_to_integer(i,ref_closed_shell_bitmask(1,1),N_int) + enddo + enddo else do i0 = elec_beta_num+1, elec_alpha_num i=occ(i0,1) call clear_bit_to_integer(i,ref_closed_shell_bitmask(1,1),N_int) enddo endif - - END_PROVIDER - BEGIN_PROVIDER [double precision, fock_op_cshell_ref_bitmask, (mo_num, mo_num) ] implicit none integer :: i0,j0,i,j,k0,k @@ -318,85 +322,117 @@ end ! ! !============================================! +BEGIN_PROVIDER [integer(bit_kind), ref_closed_shell_bitmask_kpts, (N_int,2,kpt_num)] + implicit none + integer :: i,k + do k = 1, kpt_num + do i = 1, N_int + ref_closed_shell_bitmask_kpts(i,1,k) = iand(ref_closed_shell_bitmask(i,1),kpts_bitmask(i,k)) + ref_closed_shell_bitmask_kpts(i,2,k) = iand(ref_closed_shell_bitmask(i,2),kpts_bitmask(i,k)) + enddo + enddo +END_PROVIDER + BEGIN_PROVIDER [complex*16, fock_op_cshell_ref_bitmask_kpts, (mo_num_per_kpt, mo_num_per_kpt,kpt_num) ] implicit none - integer :: i0,j0,i,j,k0,k - integer :: n_occ_ab(2) - integer :: occ(N_int*bit_kind_size,2) + integer :: i0,j0,i,j,k0,k,kblock,kvirt + integer :: n_occ_ab(2,kpt_num) + integer :: occ(N_int*bit_kind_size,2,kpt_num) integer :: n_occ_ab_virt(2) integer :: occ_virt(N_int*bit_kind_size,2) integer(bit_kind) :: key_test(N_int) integer(bit_kind) :: key_virt(N_int,2) complex*16 :: accu + complex*16, allocatable :: array_coulomb(:,:),array_exchange(:,:) - call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) - do i = 1, N_int - key_virt(i,1) = full_ijkl_bitmask(i) - key_virt(i,2) = full_ijkl_bitmask(i) - key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) - key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) + do kblock = 1,kpt_num + call bitstring_to_list_ab(ref_closed_shell_bitmask_kpts(1,1,kblock), & + occ(1,1,kblock), n_occ_ab(1,kblock), N_int) enddo - complex*16, allocatable :: array_coulomb(:),array_exchange(:) - allocate (array_coulomb(mo_num),array_exchange(mo_num)) - call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) - ! docc ---> virt single excitations - do i0 = 1, n_occ_ab(1) - i=occ(i0,1) - do j0 = 1, n_occ_ab_virt(1) - j = occ_virt(j0,1) - ! - call get_mo_two_e_integrals_coulomb_ii_complex(i,j,mo_num,array_coulomb,mo_integrals_map,mo_integrals_map_2) - ! - call get_mo_two_e_integrals_exch_ii_complex(i,j,mo_num,array_exchange,mo_integrals_map,mo_integrals_map_2) - accu = (0.d0,0.d0) - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * array_coulomb(k) - array_exchange(k) + allocate (array_coulomb(mo_num_per_kpt),array_exchange(mo_num_per_kpt)) + do kblock = 1,kpt_num + ! get virt orbs for this kpt + do i = 1, N_int + key_virt(i,1) = iand(full_ijkl_bitmask(i),kpts_bitmask(i,kblock)) + key_virt(i,2) = iand(full_ijkl_bitmask(i),kpts_bitmask(i,kblock)) + key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask_kpts(i,1,kblock)) + key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask_kpts(i,2,kblock)) enddo - fock_op_cshell_ref_bitmask_cplx(i,j) = accu + mo_one_e_integrals_complex(i,j) - !fock_op_cshell_ref_bitmask_cplx(j,i) = dconjg(accu) + mo_one_e_integrals_complex(j,i) - fock_op_cshell_ref_bitmask_cplx(j,i) = dconjg(fock_op_cshell_ref_bitmask_cplx(i,j)) - enddo - enddo - - ! virt ---> virt single excitations - do i0 = 1, n_occ_ab_virt(1) - i=occ_virt(i0,1) - do j0 = 1, n_occ_ab_virt(1) - j = occ_virt(j0,1) - call get_mo_two_e_integrals_coulomb_ii_complex(i,j,mo_num,array_coulomb,mo_integrals_map,mo_integrals_map_2) - call get_mo_two_e_integrals_exch_ii_complex(i,j,mo_num,array_exchange,mo_integrals_map,mo_integrals_map_2) - accu = (0.d0,0.d0) - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * array_coulomb(k) - array_exchange(k) + call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) + ! docc ---> virt single excitations + do i0 = 1, n_occ_ab(1,kblock) + i=occ(i0,1,kblock) + i_i = mod(i-1,mo_num_per_kpt)+1 + do j0 = 1, n_occ_ab_virt(1) + j = occ_virt(j0,1) + i_j = mod(j-1,mo_num_per_kpt)+1 + accu = (0.d0,0.d0) + do kocc = 1,kpt_num + ! + array_coulomb(1:mo_num_per_kpt) = big_array_coulomb_integrals_kpts(1:mo_num_per_kpt,kocc,i_i,i_j,kblock) + ! + array_exchange(1:mo_num_per_kpt) = big_array_exchange_integrals_kpts(1:mo_num_per_kpt,kocc,i_i,i_j,kblock) + do k0 = 1, n_occ_ab(1,kocc) + k = occ(k0,1,kocc) + i_k = mod(k-1,mo_num_per_kpt)+1 + accu += 2.d0 * array_coulomb(i_k) - array_exchange(i_k) + enddo + enddo + fock_op_cshell_ref_bitmask_cplx(i_i,i_j,kblock) = accu + mo_one_e_integrals_kpts(i_i,i_j,kblock) + !fock_op_cshell_ref_bitmask_cplx(j,i) = dconjg(accu) + mo_one_e_integrals_complex(j,i) + fock_op_cshell_ref_bitmask_cplx(i_j,i_i,kblock) = dconjg(fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock)) + enddo enddo - fock_op_cshell_ref_bitmask_cplx(i,j) = accu+ mo_one_e_integrals_complex(i,j) - fock_op_cshell_ref_bitmask_cplx(j,i) = dconjg(accu)+ mo_one_e_integrals_complex(j,i) - enddo - enddo - - ! docc ---> docc single excitations - do i0 = 1, n_occ_ab(1) - i=occ(i0,1) - do j0 = 1, n_occ_ab(1) - j = occ(j0,1) - call get_mo_two_e_integrals_coulomb_ii_complex(i,j,mo_num,array_coulomb,mo_integrals_map,mo_integrals_map_2) - call get_mo_two_e_integrals_exch_ii_complex(i,j,mo_num,array_exchange,mo_integrals_map,mo_integrals_map_2) - accu = (0.d0,0.d0) - do k0 = 1, n_occ_ab(1) - k = occ(k0,1) - accu += 2.d0 * array_coulomb(k) - array_exchange(k) + + ! virt ---> virt single excitations + do i0 = 1, n_occ_ab_virt(1) + i=occ_virt(i0,1) + i_i = mod(i-1,mo_num_per_kpt)+1 + do j0 = 1, n_occ_ab_virt(1) + j = occ_virt(j0,1) + i_j = mod(j-1,mo_num_per_kpt)+1 + accu = (0.d0,0.d0) + do kocc = 1,kpt_num + array_coulomb(1:mo_num_per_kpt) = big_array_coulomb_integrals_kpts(1:mo_num_per_kpt,kocc,i_i,i_j,kblock) + array_exchange(1:mo_num_per_kpt) = big_array_exchange_integrals_kpts(1:mo_num_per_kpt,kocc,i_i,i_j,kblock) + do k0 = 1, n_occ_ab(1,kocc) + k = occ(k0,1,kocc) + i_k = mod(k-1,mo_num_per_kpt)+1 + accu += 2.d0 * array_coulomb(i_k) - array_exchange(i_k) + enddo + enddo + fock_op_cshell_ref_bitmask_cplx(i_i,i_j,kblock) = accu + mo_one_e_integrals_kpts(i_i,i_j,kblock) + fock_op_cshell_ref_bitmask_cplx(i_j,i_i,kblock) = dconjg(fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock)) + enddo + enddo + + ! docc ---> docc single excitations + do i0 = 1, n_occ_ab(1,kblock) + i=occ(i0,1,kblock) + i_i = mod(i-1,mo_num_per_kpt)+1 + do j0 = 1, n_occ_ab(1,kblock) + j = occ(j0,1,kblock) + i_j = mod(j-1,mo_num_per_kpt)+1 + accu = (0.d0,0.d0) + do kocc = 1,kpt_num + array_coulomb(1:mo_num_per_kpt) = big_array_coulomb_integrals_kpts(1:mo_num_per_kpt,kocc,i_i,i_j,kblock) + array_exchange(1:mo_num_per_kpt) = big_array_exchange_integrals_kpts(1:mo_num_per_kpt,kocc,i_i,i_j,kblock) + do k0 = 1, n_occ_ab(1,kocc) + k = occ(k0,1,kocc) + i_k = mod(k-1,mo_num_per_kpt)+1 + accu += 2.d0 * array_coulomb(i_k) - array_exchange(i_k) + enddo + enddo + fock_op_cshell_ref_bitmask_cplx(i_i,i_j,kblock) = accu + mo_one_e_integrals_kpts(i_i,i_j,kblock) + fock_op_cshell_ref_bitmask_cplx(i_j,i_i,kblock) = dconjg(fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock)) + enddo enddo - fock_op_cshell_ref_bitmask_cplx(i,j) = accu+ mo_one_e_integrals_complex(i,j) - fock_op_cshell_ref_bitmask_cplx(j,i) = dconjg(accu)+ mo_one_e_integrals_complex(j,i) enddo - enddo deallocate(array_coulomb,array_exchange) END_PROVIDER -subroutine get_single_excitation_from_fock_complex(det_1,det_2,h,p,spin,phase,hij) +subroutine get_single_excitation_from_fock_kpts(det_1,det_2,h,p,spin,phase,hij) use bitmasks implicit none integer,intent(in) :: h,p,spin @@ -411,9 +447,10 @@ subroutine get_single_excitation_from_fock_complex(det_1,det_2,h,p,spin,phase,hi integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) integer :: i0,i complex*16 :: buffer_c(mo_num),buffer_x(mo_num) +! do do i=1, mo_num - buffer_c(i) = big_array_coulomb_integrals_complex(i,h,p) - buffer_x(i) = big_array_exchange_integrals_complex(i,h,p) + buffer_c(i) = big_array_coulomb_integrals_kpts(i,ki,h,p,kp) + buffer_x(i) = big_array_exchange_integrals_kpts(i,ki,h,p,kp) enddo do i = 1, N_int differences(i,1) = xor(det_1(i,1),ref_closed_shell_bitmask(i,1)) diff --git a/src/mo_two_e_ints/integrals_3_index.irp.f b/src/mo_two_e_ints/integrals_3_index.irp.f index 983e2642..438f4102 100644 --- a/src/mo_two_e_ints/integrals_3_index.irp.f +++ b/src/mo_two_e_ints/integrals_3_index.irp.f @@ -55,3 +55,40 @@ END_PROVIDER END_PROVIDER + BEGIN_PROVIDER [complex*16, big_array_coulomb_integrals_kpts, (mo_num_per_kpt,kpt_num,mo_num_per_kpt, mo_num_per_kpt,kpt_num)] +&BEGIN_PROVIDER [complex*16, big_array_exchange_integrals_kpts,(mo_num_per_kpt,kpt_num,mo_num_per_kpt, mo_num_per_kpt,kpt_num)] + implicit none + BEGIN_DOC + ! big_array_coulomb_integrals(j,kj,i,k,ki) = = (ik|jj) + ! big_array_exchange_integrals(j,kj,i,k,ki) = = (ij|jk) + ! for both of these, i and k must be from same kpt for integral to be nonzero + ! TODO: only loop over half, and assign two elements: + ! b_a_coul_int(j,i,k) = b_a_coul_int(j,k,i)* + ! b_a_exch_int(j,i,k) = b_a_exch_int(j,k,i)* + END_DOC + integer :: i,j,k,l + integer :: ki,kj,kk,kl + complex*16 :: get_two_e_integral_kpts + complex*16 :: integral + + do ki = 1,kpt_num + kk=ki + do k = 1, mo_num_per_kpt + do i = 1, mo_num_per_kpt + do kj = 1,kpt_num + kl=kj + do j = 1, mo_num_per_kpt + l = j + integral = get_two_e_integral_kpts(i,j,k,l,ki,kj,kk,kl,mo_integrals_map,mo_integrals_map_2) + big_array_coulomb_integrals_kpts(j,kj,i,k,ki) = integral + l = j + integral = get_two_e_integral_kpts(i,j,l,k,ki,kj,kl,kk,mo_integrals_map,mo_integrals_map_2) + big_array_exchange_integrals_kpts(j,kj,i,k,ki) = integral + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + diff --git a/src/mo_two_e_ints/map_integrals_cplx.irp.f b/src/mo_two_e_ints/map_integrals_cplx.irp.f index 66b9e3c0..1db8da3c 100644 --- a/src/mo_two_e_ints/map_integrals_cplx.irp.f +++ b/src/mo_two_e_ints/map_integrals_cplx.irp.f @@ -119,6 +119,25 @@ complex*16 function get_two_e_integral_complex(i,j,k,l,map,map2) get_two_e_integral_complex = tmp end +complex*16 function get_two_e_integral_kpts(i,j,k,l,ki,kj,kk,kl,map,map2) + use map_module + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + ! TODO: finish this + END_DOC + integer, intent(in) :: i,j,k,l + integer, intent(in) :: ki,kj,kk,kl + type(map_type), intent(inout) :: map,map2 + complex*16 :: get_two_e_integral_complex + complex*16 :: tmp + tmp = get_two_e_integral_complex( i + mo_num_per_kpt*(ki-1), & + j + mo_num_per_kpt*(kj-1), & + k + mo_num_per_kpt*(kk-1), & + l + mo_num_per_kpt*(kl-1), map,map2) + get_two_e_integral_kpts = tmp +end + complex*16 function mo_two_e_integral_complex(i,j,k,l) implicit none BEGIN_DOC @@ -133,6 +152,25 @@ complex*16 function mo_two_e_integral_complex(i,j,k,l) return end +complex*16 function mo_two_e_integral_kpts(i,j,k,l,ki,kj,kk,kl) + implicit none + BEGIN_DOC + ! Returns one integral in the MO basis + END_DOC + integer, intent(in) :: i,j,k,l + integer, intent(in) :: ki,kj,kk,kl + complex*16 :: get_two_e_integral_complex + PROVIDE mo_two_e_integrals_in_map mo_integrals_cache_complex + PROVIDE mo_two_e_integrals_in_map + !DIR$ FORCEINLINE + mo_two_e_integral_kpts = get_two_e_integral_complex( & + i + mo_num_per_kpt*(ki-1), & + j + mo_num_per_kpt*(kj-1), & + k + mo_num_per_kpt*(kk-1), & + l + mo_num_per_kpt*(kl-1),mo_integrals_map,mo_integrals_map_2) + return +end + subroutine get_mo_two_e_integrals_complex(j,k,l,sze,out_val,map,map2) use map_module implicit none