diff --git a/src/hartree_fock/hf_energy.irp.f b/src/hartree_fock/hf_energy.irp.f index 5a68164f..ea5e4b1c 100644 --- a/src/hartree_fock/hf_energy.irp.f +++ b/src/hartree_fock/hf_energy.irp.f @@ -13,19 +13,22 @@ END_PROVIDER BEGIN_PROVIDER [ double precision, hf_energy] &BEGIN_PROVIDER [ double precision, hf_two_electron_energy] +&BEGIN_PROVIDER [ double precision, hf_two_electron_energy_jk, (2)] &BEGIN_PROVIDER [ double precision, hf_one_electron_energy] implicit none BEGIN_DOC ! Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components. END_DOC - integer :: i,j,k + integer :: i,j,k,jk hf_energy = nuclear_repulsion hf_two_electron_energy = 0.d0 + hf_two_electron_energy_jk = 0.d0 hf_one_electron_energy = 0.d0 if (is_complex) then - complex*16 :: hf_1e_tmp, hf_2e_tmp + complex*16 :: hf_1e_tmp, hf_2e_tmp, hf_2e_tmp_jk(2) hf_1e_tmp = (0.d0,0.d0) hf_2e_tmp = (0.d0,0.d0) + hf_2e_tmp_jk = (0.d0,0.d0) do k=1,kpt_num do j=1,ao_num_per_kpt do i=1,ao_num_per_kpt @@ -33,9 +36,21 @@ END_PROVIDER +ao_two_e_integral_beta_kpts(i,j,k) * scf_density_matrix_ao_beta_kpts(j,i,k) ) hf_1e_tmp += ao_one_e_integrals_kpts(i,j,k) * (scf_density_matrix_ao_alpha_kpts(j,i,k) & + scf_density_matrix_ao_beta_kpts (j,i,k) ) + do jk=1,2 + hf_2e_tmp_jk(jk) += 0.5d0 * ( ao_two_e_integral_alpha_kpts_jk(i,j,k,jk) * scf_density_matrix_ao_alpha_kpts(j,i,k) & + +ao_two_e_integral_beta_kpts_jk(i,j,k,jk) * scf_density_matrix_ao_beta_kpts(j,i,k) ) + enddo enddo enddo enddo + do jk=1,2 + if (dabs(dimag(hf_2e_tmp_jk(jk))).gt.1.d-10) then + print*,'HF_2e energy (jk) should be real:',jk,irp_here + stop -1 + else + hf_two_electron_energy_jk(jk) = dble(hf_2e_tmp_jk(jk)) + endif + enddo if (dabs(dimag(hf_2e_tmp)).gt.1.d-10) then print*,'HF_2e energy should be real:',irp_here stop -1 diff --git a/src/hartree_fock/print_e_scf.irp.f b/src/hartree_fock/print_e_scf.irp.f index 989c0b9c..8efeba0d 100644 --- a/src/hartree_fock/print_e_scf.irp.f +++ b/src/hartree_fock/print_e_scf.irp.f @@ -15,6 +15,9 @@ subroutine run print*,hf_one_electron_energy print*,hf_two_electron_energy print*,hf_energy + print*,'hf 2e J,K energy' + print*,hf_two_electron_energy_jk(1) + print*,hf_two_electron_energy_jk(2) end diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index e2ada6fc..2f7d4a1a 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -542,27 +542,26 @@ BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_kpts, (ao_num_per_kpt, ao_num_per_kp endif END_PROVIDER - - BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] -&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts_jk, (ao_num_per_kpt, ao_num_per_kpt, kpt_num, 2) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts_jk , (ao_num_per_kpt, ao_num_per_kpt, kpt_num, 2) ] use map_module implicit none BEGIN_DOC - ! Alpha and Beta Fock matrices in AO basis set + ! Alpha and Beta Fock matrices in AO basis set separated into j/k END_DOC !TODO: finish implementing this: see complex qp1 (different mapping) - + integer :: i,j,k,l,k1,r,s integer :: i0,j0,k0,l0 integer*8 :: p,q complex*16 :: integral, c0 - complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:) - complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:) - - ao_two_e_integral_alpha_kpts = (0.d0,0.d0) - ao_two_e_integral_beta_kpts = (0.d0,0.d0) + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:,:) + + ao_two_e_integral_alpha_kpts_jk = (0.d0,0.d0) + ao_two_e_integral_beta_kpts_jk = (0.d0,0.d0) PROVIDE ao_two_e_integrals_in_map scf_density_matrix_ao_alpha_kpts scf_density_matrix_ao_beta_kpts - + integer(omp_lock_kind) :: lck(ao_num) integer(map_size_kind) :: i8 integer :: ii(4), jj(4), kk(4), ll(4), k2 @@ -572,7 +571,267 @@ END_PROVIDER complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) integer(key_kind) :: key1 integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l - + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts_jk, ao_two_e_integral_beta_kpts_jk) + + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_1(ii,jj,kk,ll,key1) + ! i<=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + + if (shiftl(key1,1)==keys(k1)) then !imaginary part (even) + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = i_sign(k2)*values(k1) !for klij and lkji, take complex conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = values(k1) + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_kpts_jk += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts_jk += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & + !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & + !$OMP kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l, & + !$OMP c0,key1)& + !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & + !$OMP SCF_density_matrix_ao_beta_kpts, & + !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts_jk, ao_two_e_integral_beta_kpts_jk) + + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) + allocate(keys(n_elements_max), values(n_elements_max)) + allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2), & + ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num,2)) + ao_two_e_integral_alpha_tmp = (0.d0,0.d0) + ao_two_e_integral_beta_tmp = (0.d0,0.d0) + + !$OMP DO SCHEDULE(static,1) + do i8=0_8,ao_integrals_map_2%map_size + n_elements = n_elements_max + call get_cache_map(ao_integrals_map_2,i8,keys,values,n_elements) + do k1=1,n_elements + ! get original key + ! reverse of 2*key (imag part) and 2*key-1 (real part) + key1 = shiftr(keys(k1)+1,1) + + call two_e_integrals_index_reverse_complex_2(ii,jj,kk,ll,key1) + ! i>=k, j<=l, ik<=jl + ! ijkl, jilk, klij*, lkji* + if (shiftl(key1,1)==keys(k1)) then !imaginary part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = i_sign(k2)*values(k1) ! for klij and lkji, take conjugate + + !G_a(i,k) += D_{ab}(l,j)*() + !G_b(i,k) += D_{ab}(l,j)*() + !G_a(i,l) -= D_a (k,j)*() + !G_b(i,l) -= D_b (k,j)*() + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + else ! real part + do k2=1,4 + if (ii(k2)==0) then + cycle + endif + i = ii(k2) + j = jj(k2) + k = kk(k2) + l = ll(k2) + call get_kpt_idx_ao(i,kpt_i,idx_i) + call get_kpt_idx_ao(j,kpt_j,idx_j) + call get_kpt_idx_ao(k,kpt_k,idx_k) + call get_kpt_idx_ao(l,kpt_l,idx_l) + integral = values(k1) + + if (kpt_l.eq.kpt_j) then + c0 = (scf_density_matrix_ao_alpha_kpts(idx_l,idx_j,kpt_j)+scf_density_matrix_ao_beta_kpts(idx_l,idx_j,kpt_j))*integral + if(kpt_i.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_k,kpt_i,1) += c0 + ao_two_e_integral_beta_tmp (idx_i,idx_k,kpt_i,1) += c0 + endif + + if (kpt_l.eq.kpt_i) then + if(kpt_j.ne.kpt_k) then + print*,'problem in ',irp_here,' ikjl: ',kpt_i,kpt_k,kpt_j,kpt_l + stop 1 + endif + ao_two_e_integral_alpha_tmp(idx_i,idx_l,kpt_i,2) -= SCF_density_matrix_ao_alpha_kpts(idx_k,idx_j,kpt_j) * integral + ao_two_e_integral_beta_tmp (idx_i,idx_l,kpt_i,2) -= scf_density_matrix_ao_beta_kpts (idx_k,idx_j,kpt_j) * integral + endif + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_kpts_jk += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_kpts_jk += ao_two_e_integral_beta_tmp + !$OMP END CRITICAL + deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) + !$OMP END PARALLEL + + +END_PROVIDER + + + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_kpts, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_kpts , (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ] + use map_module + implicit none + BEGIN_DOC + ! Alpha and Beta Fock matrices in AO basis set + END_DOC + !TODO: finish implementing this: see complex qp1 (different mapping) + + integer :: i,j,k,l,k1,r,s + integer :: i0,j0,k0,l0 + integer*8 :: p,q + complex*16 :: integral, c0 + complex*16, allocatable :: ao_two_e_integral_alpha_tmp(:,:,:) + complex*16, allocatable :: ao_two_e_integral_beta_tmp(:,:,:) + + ao_two_e_integral_alpha_kpts = (0.d0,0.d0) + ao_two_e_integral_beta_kpts = (0.d0,0.d0) + PROVIDE ao_two_e_integrals_in_map scf_density_matrix_ao_alpha_kpts scf_density_matrix_ao_beta_kpts + + integer(omp_lock_kind) :: lck(ao_num) + integer(map_size_kind) :: i8 + integer :: ii(4), jj(4), kk(4), ll(4), k2 + integer(cache_map_size_kind) :: n_elements_max, n_elements + integer(key_kind), allocatable :: keys(:) + double precision, allocatable :: values(:) + complex*16, parameter :: i_sign(4) = (/(0.d0,1.d0),(0.d0,1.d0),(0.d0,-1.d0),(0.d0,-1.d0)/) + integer(key_kind) :: key1 + integer :: kpt_i,kpt_j,kpt_k,kpt_l,idx_i,idx_j,idx_k,idx_l + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & @@ -581,14 +840,14 @@ END_PROVIDER !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts, kpt_num, irp_here, & !$OMP SCF_density_matrix_ao_beta_kpts, & !$OMP ao_integrals_map, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) - + call get_cache_map_n_elements_max(ao_integrals_map,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) - + !$OMP DO SCHEDULE(static,1) do i8=0_8,ao_integrals_map%map_size n_elements = n_elements_max @@ -686,7 +945,7 @@ END_PROVIDER deallocate(keys,values,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp) !$OMP END PARALLEL - + !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, & !$OMP n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, & @@ -695,14 +954,14 @@ END_PROVIDER !$OMP SHARED(ao_num_per_kpt,SCF_density_matrix_ao_alpha_kpts,kpt_num, irp_here, & !$OMP SCF_density_matrix_ao_beta_kpts, & !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_kpts, ao_two_e_integral_beta_kpts) - + call get_cache_map_n_elements_max(ao_integrals_map_2,n_elements_max) allocate(keys(n_elements_max), values(n_elements_max)) allocate(ao_two_e_integral_alpha_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num), & ao_two_e_integral_beta_tmp(ao_num_per_kpt,ao_num_per_kpt,kpt_num)) ao_two_e_integral_alpha_tmp = (0.d0,0.d0) ao_two_e_integral_beta_tmp = (0.d0,0.d0) - + !$OMP DO SCHEDULE(static,1) do i8=0_8,ao_integrals_map_2%map_size n_elements = n_elements_max