From e638a640f07ec20383e7a701957893a9b277d6c6 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Tue, 24 Mar 2020 16:43:04 -0500 Subject: [PATCH] problem with 1rdm kpts --- src/cipsi/pt2_stoch_routines.irp.f | 3 +- src/davidson/print_e_components.irp.f | 20 +- src/determinants/density_matrix_cplx.irp.f | 382 ++++++++++++++++++++- src/determinants/single_excitations.irp.f | 155 +++++---- src/determinants/slater_rules.irp.f | 3 +- 5 files changed, 482 insertions(+), 81 deletions(-) diff --git a/src/cipsi/pt2_stoch_routines.irp.f b/src/cipsi/pt2_stoch_routines.irp.f index 635353c5..918f26ed 100644 --- a/src/cipsi/pt2_stoch_routines.irp.f +++ b/src/cipsi/pt2_stoch_routines.irp.f @@ -177,7 +177,8 @@ subroutine ZMQ_pt2(E, pt2,relative_error, error, variance, norm, N_in) if (is_complex) then !todo: psi_selectors isn't linked to psi_selectors_coef anymore; should we provide both? - PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals_complex pt2_w + !PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals_complex pt2_w + PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals_kpts pt2_w PROVIDE psi_selectors pt2_u pt2_J pt2_R else PROVIDE nproc pt2_F mo_two_e_integrals_in_map mo_one_e_integrals pt2_w diff --git a/src/davidson/print_e_components.irp.f b/src/davidson/print_e_components.irp.f index 9c3caf23..fcea369d 100644 --- a/src/davidson/print_e_components.irp.f +++ b/src/davidson/print_e_components.irp.f @@ -6,7 +6,7 @@ subroutine print_energy_components() integer, save :: ifirst = 0 double precision :: Vee, Ven, Vnn, Vecp, T, f complex*16 :: fc - integer :: i,j,k + integer :: i,j,k,kk Vnn = nuclear_repulsion @@ -20,12 +20,18 @@ subroutine print_energy_components() T = 0.d0 if (is_complex) then - do j=1,mo_num - do i=1,mo_num - fc = one_e_dm_mo_alpha_complex(i,j,k) + one_e_dm_mo_beta_complex(i,j,k) - Ven = Ven + dble(fc * mo_integrals_n_e_complex(j,i)) - Vecp = Vecp + dble(fc * mo_pseudo_integrals_complex(j,i)) - T = T + dble(fc * mo_kinetic_integrals_complex(j,i)) + do kk=1,kpt_num + do j=1,mo_num_per_kpt + do i=1,mo_num_per_kpt + !fc = one_e_dm_mo_alpha_complex(i,j,k) + one_e_dm_mo_beta_complex(i,j,k) + !Ven = Ven + dble(fc * mo_integrals_n_e_complex(j,i)) + !Vecp = Vecp + dble(fc * mo_pseudo_integrals_complex(j,i)) + !T = T + dble(fc * mo_kinetic_integrals_complex(j,i)) + fc = one_e_dm_mo_alpha_kpts(i,j,kk,k) + one_e_dm_mo_beta_kpts(i,j,kk,k) + Ven = Ven + dble(fc * mo_integrals_n_e_kpts(j,i,kk)) + Vecp = Vecp + dble(fc * mo_pseudo_integrals_kpts(j,i,kk)) + T = T + dble(fc * mo_kinetic_integrals_kpts(j,i,kk)) + enddo enddo enddo else diff --git a/src/determinants/density_matrix_cplx.irp.f b/src/determinants/density_matrix_cplx.irp.f index 18bae7db..c1b63cf2 100644 --- a/src/determinants/density_matrix_cplx.irp.f +++ b/src/determinants/density_matrix_cplx.irp.f @@ -290,8 +290,8 @@ END_PROVIDER integer :: i,j,k,l complex*16 :: mo_alpha,mo_beta - one_e_dm_ao_alpha = (0.d0,0.d0) - one_e_dm_ao_beta = (0.d0,0.d0) + one_e_dm_ao_alpha_complex = (0.d0,0.d0) + one_e_dm_ao_beta_complex = (0.d0,0.d0) do k = 1, ao_num do l = 1, ao_num do i = 1, mo_num @@ -309,3 +309,381 @@ END_PROVIDER END_PROVIDER +!============================================! +! ! +! kpts ! +! ! +!============================================! + + BEGIN_PROVIDER [ complex*16, one_e_dm_mo_alpha_average_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num) ] +&BEGIN_PROVIDER [ complex*16, one_e_dm_mo_beta_average_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! $\alpha$ and $\beta$ one-body density matrix for each state + END_DOC + integer :: i,k + one_e_dm_mo_alpha_average_kpts = (0.d0,0.d0) + one_e_dm_mo_beta_average_kpts = (0.d0,0.d0) + do i = 1,N_states + do k=1,kpt_num + one_e_dm_mo_alpha_average_kpts(:,:,k) += one_e_dm_mo_alpha_kpts(:,:,k,i) * state_average_weight(i) + one_e_dm_mo_beta_average_kpts(:,:,k) += one_e_dm_mo_beta_kpts(:,:,k,i) * state_average_weight(i) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, one_e_dm_mo_diff_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Difference of the one-body density matrix with respect to the ground state + END_DOC + integer :: i,j, istate,k + + do istate=2,N_states + do k=1,kpt_num + do j=1,mo_num_per_kpt + do i=1,mo_num_per_kpt + one_e_dm_mo_diff_kpts(i,j,k,istate) = & + one_e_dm_mo_alpha_kpts(i,j,k,istate) - one_e_dm_mo_alpha_kpts(i,j,k,1) +& + one_e_dm_mo_beta_kpts (i,j,k,istate) - one_e_dm_mo_beta_kpts (i,j,k,1) + enddo + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, one_e_dm_mo_spin_index_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num,N_states,2) ] + implicit none + integer :: i,j,k,ispin,istate + ispin = 1 + do istate = 1, N_states + do k=1,kpt_num + do j = 1, mo_num_per_kpt + do i = 1, mo_num_per_kpt + one_e_dm_mo_spin_index_kpts(i,j,k,istate,ispin) = one_e_dm_mo_alpha_kpts(i,j,k,istate) + enddo + enddo + enddo + enddo + + ispin = 2 + do istate = 1, N_states + do k=1,kpt_num + do j = 1, mo_num_per_kpt + do i = 1, mo_num_per_kpt + one_e_dm_mo_spin_index_kpts(i,j,k,istate,ispin) = one_e_dm_mo_beta_kpts(i,j,k,istate) + enddo + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, one_e_dm_dagger_mo_spin_index_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num,N_states,2) ] + print*,irp_here,' not implemented for kpts' + stop -1 +! implicit none +! integer :: i,j,ispin,istate +! ispin = 1 +! do istate = 1, N_states +! do j = 1, mo_num +! one_e_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_e_dm_mo_alpha(j,j,istate) +! do i = j+1, mo_num +! one_e_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_e_dm_mo_alpha(i,j,istate) +! one_e_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_e_dm_mo_alpha(i,j,istate) +! enddo +! enddo +! enddo +! +! ispin = 2 +! do istate = 1, N_states +! do j = 1, mo_num +! one_e_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_e_dm_mo_beta(j,j,istate) +! do i = j+1, mo_num +! one_e_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_e_dm_mo_beta(i,j,istate) +! one_e_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_e_dm_mo_beta(i,j,istate) +! enddo +! enddo +! enddo +! +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, one_e_dm_mo_alpha_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num,N_states) ] +&BEGIN_PROVIDER [ complex*16, one_e_dm_mo_beta_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num,N_states) ] + implicit none + BEGIN_DOC + ! $\alpha$ and $\beta$ one-body density matrix for each state + ! $\gamma_{\mu\nu} = \langle\Psi|a_{\nu}^{\dagger}a_{\mu}|\Psi\rangle$ + ! $\gamma_{\mu\nu} = \langle a_{\nu} \Psi|a_{\mu} \Psi\rangle$ + ! $\gamma_{\mu\nu} = \sum_{IJ} c^*_J c_I \langle a_{\nu} I|a_{\mu} J\rangle$ + END_DOC + !todo: implement for kpts + integer :: j,k,l,m,k_a,k_b + integer :: occ(N_int*bit_kind_size,2) + complex*16 :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: ih1,ip1,kh1,kp1,kk,k_shft,ii + integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int) + integer(bit_kind) :: tmp_det_kpts(N_int,2) + integer :: exc(0:2,2),n_occ(2) + complex*16, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:) + integer :: krow, kcol, lrow, lcol + + PROVIDE psi_det psi_coef_complex + + one_e_dm_mo_alpha_kpts = (0.d0,0.d0) + one_e_dm_mo_beta_kpts = (0.d0,0.d0) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,& + !$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2,ih1,ip1,kh1,kp1,kk,& + !$OMP tmp_det_kpts,k_shft,ii)& + !$OMP SHARED(psi_det,psi_coef_complex,N_int,N_states,elec_alpha_num_kpts, & + !$OMP elec_beta_num_kpts,one_e_dm_mo_alpha_kpts,one_e_dm_mo_beta_kpts,N_det,& + !$OMP mo_num_per_kpt,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,& + !$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,& + !$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,& + !$OMP psi_bilinear_matrix_values_complex, psi_bilinear_matrix_transp_values_complex,& + !$OMP N_det_alpha_unique,N_det_beta_unique,irp_here,kpt_num,kpts_bitmask) + allocate(tmp_a(mo_num_per_kpt,mo_num_per_kpt,kpt_num,N_states), tmp_b(mo_num_per_kpt,mo_num_per_kpt,kpt_num,N_states) ) + tmp_a = (0.d0,0.d0) + !$OMP DO SCHEDULE(dynamic,64) + do k_a=1,N_det + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) + + ! Diagonal part + ! ------------- + + do kk=1,kpt_num + k_shft = (kk-1)*mo_num_per_kpt + do ii=1,N_int + tmp_det_kpts(ii,1) = iand(tmp_det(ii,1),kpts_bitmask(ii,kk)) + tmp_det_kpts(ii,2) = iand(tmp_det(ii,2),kpts_bitmask(ii,kk)) + enddo + call bitstring_to_list_ab(tmp_det_kpts, occ, n_occ, N_int) + do m=1,N_states + ck = cdabs(psi_bilinear_matrix_values_complex(k_a,m)*psi_bilinear_matrix_values_complex(k_a,m)) + !do l=1,elec_alpha_num_kpts(kk) + do l=1,n_occ(1) + j = occ(l,1) - k_shft + tmp_a(j,j,kk,m) += ck + enddo + enddo + enddo + + if (k_a == N_det) cycle + l = k_a+1 + lrow = psi_bilinear_matrix_rows(l) + lcol = psi_bilinear_matrix_columns(l) + ! Fix beta determinant, loop over alphas + do while ( lcol == kcol ) + tmp_det2(:) = psi_det_alpha_unique(:, lrow) + call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int) + if (degree == 1) then + exc = 0 + call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int) + call decode_exc_spin(exc,h1,p1,h2,p2) + ! h1 occ in k + ! p1 occ in l + ih1 = mod(h1-1,mo_num_per_kpt)+1 + ip1 = mod(p1-1,mo_num_per_kpt)+1 + kh1 = (h1-1)/mo_num_per_kpt + 1 + kp1 = (p1-1)/mo_num_per_kpt + 1 + if (kh1.ne.kp1) then + print *,'problem in: ',irp_here,'a' + print *,' h1 = ',h1 + print *,' p1 = ',p1 + print *,'ih1 = ',ih1 + print *,'ip1 = ',ip1 + print *,'kh1 = ',kh1 + print *,'kp1 = ',kp1 + !call debug_det(tmp_det,N_int) + !call debug_spindet(tmp_det2,N_int) + !call print_spindet(tmp_det2,N_int) + !stop -2 + endif + do m=1,N_states + ckl = dconjg(psi_bilinear_matrix_values_complex(k_a,m))*psi_bilinear_matrix_values_complex(l,m) * phase + tmp_a(ih1,ip1,kh1,m) += dconjg(ckl) + tmp_a(ip1,ih1,kh1,m) += ckl + enddo + endif + l = l+1 + if (l>N_det) exit + lrow = psi_bilinear_matrix_rows(l) + lcol = psi_bilinear_matrix_columns(l) + enddo + + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + one_e_dm_mo_alpha_kpts(:,:,:,:) = one_e_dm_mo_alpha_kpts(:,:,:,:) + tmp_a(:,:,:,:) + !$OMP END CRITICAL + deallocate(tmp_a) + + tmp_b = (0.d0,0.d0) + !$OMP DO SCHEDULE(dynamic,64) + do k_b=1,N_det + krow = psi_bilinear_matrix_transp_rows(k_b) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_transp_columns(k_b) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) + + ! Diagonal part + ! ------------- + + do kk=1,kpt_num + k_shft = (kk-1)*mo_num_per_kpt + do ii=1,N_int + tmp_det_kpts(ii,1) = iand(tmp_det(ii,1),kpts_bitmask(ii,kk)) + tmp_det_kpts(ii,2) = iand(tmp_det(ii,2),kpts_bitmask(ii,kk)) + enddo + call bitstring_to_list_ab(tmp_det_kpts, occ, n_occ, N_int) + do m=1,N_states + ck = cdabs(psi_bilinear_matrix_transp_values_complex(k_b,m)*psi_bilinear_matrix_transp_values_complex(k_b,m)) + do l=1,n_occ(2) + j = occ(l,2) - k_shft + tmp_b(j,j,kk,m) += ck + enddo + enddo + enddo + + if (k_b == N_det) cycle + l = k_b+1 + lrow = psi_bilinear_matrix_transp_rows(l) + lcol = psi_bilinear_matrix_transp_columns(l) + ! Fix beta determinant, loop over alphas + do while ( lrow == krow ) + tmp_det2(:) = psi_det_beta_unique(:, lcol) + call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int) + if (degree == 1) then + exc = 0 + call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int) + call decode_exc_spin(exc,h1,p1,h2,p2) + ih1 = mod(h1-1,mo_num_per_kpt)+1 + ip1 = mod(p1-1,mo_num_per_kpt)+1 + kh1 = (h1-1)/mo_num_per_kpt + 1 + kp1 = (p1-1)/mo_num_per_kpt + 1 + if (kh1.ne.kp1) then + print *,'problem in: ',irp_here,'b' + print *,' h1 = ',h1 + print *,' p1 = ',p1 + print *,'ih1 = ',ih1 + print *,'ip1 = ',ip1 + print *,'kh1 = ',kh1 + print *,'kp1 = ',kp1 + !stop -3 + endif + do m=1,N_states + ckl = dconjg(psi_bilinear_matrix_transp_values_complex(k_b,m))*psi_bilinear_matrix_transp_values_complex(l,m) * phase + tmp_b(ih1,ip1,kh1,m) += dconjg(ckl) + tmp_b(ip1,ih1,kh1,m) += ckl + enddo + endif + l = l+1 + if (l>N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l) + lcol = psi_bilinear_matrix_transp_columns(l) + enddo + + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_e_dm_mo_beta_kpts(:,:,:,:) = one_e_dm_mo_beta_kpts(:,:,:,:) + tmp_b(:,:,:,:) + !$OMP END CRITICAL + + deallocate(tmp_b) + !$OMP END PARALLEL + +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, one_e_dm_mo_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! One-body density matrix + END_DOC + one_e_dm_mo_kpts = one_e_dm_mo_alpha_average_kpts + one_e_dm_mo_beta_average_kpts +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, one_e_spin_density_mo_kpts, (mo_num_per_kpt,mo_num_per_kpt,kpt_num) ] + implicit none + BEGIN_DOC + ! $\rho(\alpha) - \rho(\beta)$ + END_DOC + one_e_spin_density_mo_kpts = one_e_dm_mo_alpha_average_kpts - one_e_dm_mo_beta_average_kpts +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, one_e_spin_density_ao_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ] + BEGIN_DOC + ! One body spin density matrix on the |AO| basis : $\rho_{AO}(\alpha) - \rho_{AO}(\beta)$ + ! todo: verify that this is correct for complex + ! equivalent to using mo_to_ao_no_overlap? + END_DOC + implicit none + integer :: i,j,k,l,kk + complex*16 :: dm_mo + + one_e_spin_density_ao_kpts = (0.d0,0.d0) + do kk=1,kpt_num + do k = 1, ao_num_per_kpt + do l = 1, ao_num_per_kpt + do i = 1, mo_num_per_kpt + do j = 1, mo_num_per_kpt + dm_mo = one_e_spin_density_mo_kpts(j,i,kk) + ! if(dabs(dm_mo).le.1.d-10)cycle + one_e_spin_density_ao_kpts(l,k,kk) += dconjg(mo_coef_kpts(k,i,kk)) * mo_coef_kpts(l,j,kk) * dm_mo + + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, one_e_dm_ao_alpha_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ] +&BEGIN_PROVIDER [ complex*16, one_e_dm_ao_beta_kpts, (ao_num_per_kpt,ao_num_per_kpt,kpt_num) ] + BEGIN_DOC + ! One body density matrix on the |AO| basis : $\rho_{AO}(\alpha), \rho_{AO}(\beta)$. + END_DOC + implicit none + integer :: i,j,k,l,kk + complex*16 :: mo_alpha,mo_beta + + one_e_dm_ao_alpha_kpts = (0.d0,0.d0) + one_e_dm_ao_beta_kpts = (0.d0,0.d0) + do kk=1,kpt_num + do k = 1, ao_num_per_kpt + do l = 1, ao_num_per_kpt + do i = 1, mo_num_per_kpt + do j = 1, mo_num_per_kpt + mo_alpha = one_e_dm_mo_alpha_average_kpts(j,i,kk) + mo_beta = one_e_dm_mo_beta_average_kpts(j,i,kk) + ! if(dabs(dm_mo).le.1.d-10)cycle + one_e_dm_ao_alpha_kpts(l,k,kk) += dconjg(mo_coef_kpts(k,i,kk)) * mo_coef_kpts(l,j,kk) * mo_alpha + one_e_dm_ao_beta_kpts(l,k,kk) += dconjg(mo_coef_kpts(k,i,kk)) * mo_coef_kpts(l,j,kk) * mo_beta + enddo + enddo + enddo + enddo + enddo + +END_PROVIDER + + diff --git a/src/determinants/single_excitations.irp.f b/src/determinants/single_excitations.irp.f index 6e394572..044c7d06 100644 --- a/src/determinants/single_excitations.irp.f +++ b/src/determinants/single_excitations.irp.f @@ -336,6 +336,7 @@ 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,kblock,kvirt + integer :: i_i, i_j, i_k, kocc integer :: n_occ_ab(2,kpt_num) integer :: occ(N_int*bit_kind_size,2,kpt_num) integer :: n_occ_ab_virt(2) @@ -343,7 +344,7 @@ BEGIN_PROVIDER [complex*16, fock_op_cshell_ref_bitmask_kpts, (mo_num_per_kpt, mo 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(:,:) + complex*16, allocatable :: array_coulomb(:),array_exchange(:) do kblock = 1,kpt_num call bitstring_to_list_ab(ref_closed_shell_bitmask_kpts(1,1,kblock), & @@ -378,9 +379,9 @@ BEGIN_PROVIDER [complex*16, fock_op_cshell_ref_bitmask_kpts, (mo_num_per_kpt, mo 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_kpts(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)) + fock_op_cshell_ref_bitmask_kpts(i_j,i_i,kblock) = dconjg(fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock)) enddo enddo @@ -401,8 +402,8 @@ BEGIN_PROVIDER [complex*16, fock_op_cshell_ref_bitmask_kpts, (mo_num_per_kpt, mo 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)) + fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock) = accu + mo_one_e_integrals_kpts(i_i,i_j,kblock) + fock_op_cshell_ref_bitmask_kpts(i_j,i_i,kblock) = dconjg(fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock)) enddo enddo @@ -423,8 +424,8 @@ BEGIN_PROVIDER [complex*16, fock_op_cshell_ref_bitmask_kpts, (mo_num_per_kpt, mo 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)) + fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock) = accu + mo_one_e_integrals_kpts(i_i,i_j,kblock) + fock_op_cshell_ref_bitmask_kpts(i_j,i_i,kblock) = dconjg(fock_op_cshell_ref_bitmask_kpts(i_i,i_j,kblock)) enddo enddo enddo @@ -432,69 +433,83 @@ BEGIN_PROVIDER [complex*16, fock_op_cshell_ref_bitmask_kpts, (mo_num_per_kpt, mo END_PROVIDER -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 - double precision, intent(in) :: phase - integer(bit_kind), intent(in) :: det_1(N_int,2), det_2(N_int,2) - complex*16, intent(out) :: hij - integer(bit_kind) :: differences(N_int,2) - integer(bit_kind) :: hole(N_int,2) - integer(bit_kind) :: partcl(N_int,2) - integer :: occ_hole(N_int*bit_kind_size,2) - integer :: occ_partcl(N_int*bit_kind_size,2) - 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_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)) - differences(i,2) = xor(det_1(i,2),ref_closed_shell_bitmask(i,2)) - hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) - hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) - partcl(i,1) = iand(differences(i,1),det_1(i,1)) - partcl(i,2) = iand(differences(i,2),det_1(i,2)) - enddo - call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) - call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) - hij = fock_op_cshell_ref_bitmask_cplx(h,p) - ! holes :: direct terms - do i0 = 1, n_occ_ab_hole(1) - i = occ_hole(i0,1) - hij -= buffer_c(i) - enddo - do i0 = 1, n_occ_ab_hole(2) - i = occ_hole(i0,2) - hij -= buffer_c(i) - enddo - - ! holes :: exchange terms - do i0 = 1, n_occ_ab_hole(spin) - i = occ_hole(i0,spin) - hij += buffer_x(i) - enddo - - ! particles :: direct terms - do i0 = 1, n_occ_ab_partcl(1) - i = occ_partcl(i0,1) - hij += buffer_c(i) - enddo - do i0 = 1, n_occ_ab_partcl(2) - i = occ_partcl(i0,2) - hij += buffer_c(i) - enddo - - ! particles :: exchange terms - do i0 = 1, n_occ_ab_partcl(spin) - i = occ_partcl(i0,spin) - hij -= buffer_x(i) - enddo - hij = hij * phase +subroutine get_single_excitation_from_fock_kpts(det_1,det_2,ih,ip,spin,phase,hij) + use bitmasks + !called by i_h_j{,_s2,_single_spin}_complex + ! ih, ip are indices in total mo list (not per kpt) + implicit none + integer,intent(in) :: ih,ip,spin + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: det_1(N_int,2), det_2(N_int,2) + complex*16, intent(out) :: hij + integer(bit_kind) :: differences(N_int,2) + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: partcl(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_partcl(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i,h,p + integer :: ki,khp + complex*16 :: buffer_c(mo_num_per_kpt),buffer_x(mo_num_per_kpt) + khp = (ip-1)/mo_num_per_kpt+1 + p = mod(ip-1,mo_num_per_kpt)+1 + h = mod(ih-1,mo_num_per_kpt)+1 + !todo: omp kpts + do ki=1,kpt_num + do i=1, mo_num_per_kpt + ! + buffer_c(i) = big_array_coulomb_integrals_kpts(i,ki,h,p,khp) + ! + buffer_x(i) = big_array_exchange_integrals_kpts(i,ki,h,p,khp) + enddo + do i = 1, N_int + !holes in ref, not in det1 + !part in det1, not in ref + differences(i,1) = iand(xor(det_1(i,1),ref_closed_shell_bitmask(i,1)),kpts_bitmask(i,ki)) + differences(i,2) = iand(xor(det_1(i,2),ref_closed_shell_bitmask(i,2)),kpts_bitmask(i,ki)) + !differences(i,1) = xor(det_1(i,1),ref_closed_shell_bitmask_kpts(i,1,ki)) + !differences(i,2) = xor(det_1(i,2),ref_closed_shell_bitmask_kpts(i,2,ki)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask_kpts(i,1,ki)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask_kpts(i,2,ki)) + partcl(i,1) = iand(differences(i,1),det_1(i,1)) + partcl(i,2) = iand(differences(i,2),det_1(i,2)) + enddo + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + hij = fock_op_cshell_ref_bitmask_kpts(h,p,khp) + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) - (ki-1)*mo_num_per_kpt + hij -= buffer_c(i) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) - (ki-1)*mo_num_per_kpt + hij -= buffer_c(i) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) - (ki-1)*mo_num_per_kpt + hij += buffer_x(i) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) - (ki-1)*mo_num_per_kpt + hij += buffer_c(i) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) - (ki-1)*mo_num_per_kpt + hij += buffer_c(i) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) - (ki-1)*mo_num_per_kpt + hij -= buffer_x(i) + enddo + enddo + hij = hij * phase end diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 723f3194..ddd469b1 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -2491,7 +2491,8 @@ subroutine i_H_j_complex(key_i,key_j,Nint,hij) p = exc(1,2,2) spin = 2 endif - call get_single_excitation_from_fock_complex(key_i,key_j,m,p,spin,phase,hij) + !call get_single_excitation_from_fock_complex(key_i,key_j,m,p,spin,phase,hij) + call get_single_excitation_from_fock_kpts(key_i,key_j,m,p,spin,phase,hij) case (0) hij = dcmplx(diag_H_mat_elem(key_i,Nint),0.d0)