From b1e14142c655e6400a975273f506c62702996d13 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Wed, 5 Feb 2020 17:50:17 -0600 Subject: [PATCH] working on complex MO 2e ints --- src/ao_two_e_ints/map_integrals_complex.irp.f | 8 +- src/mo_basis/mos_complex.irp.f | 32 +- .../four_idx_novvvv_complex.irp.f | 247 ++++ src/mo_two_e_ints/map_integrals_complex.irp.f | 6 +- src/mo_two_e_ints/mo_bi_integrals.irp.f | 204 ++- .../mo_bi_integrals_complex.irp.f | 1163 +++++++++++++++++ 6 files changed, 1644 insertions(+), 16 deletions(-) create mode 100644 src/mo_two_e_ints/four_idx_novvvv_complex.irp.f create mode 100644 src/mo_two_e_ints/mo_bi_integrals_complex.irp.f diff --git a/src/ao_two_e_ints/map_integrals_complex.irp.f b/src/ao_two_e_ints/map_integrals_complex.irp.f index f087e1b8..611bc4cb 100644 --- a/src/ao_two_e_ints/map_integrals_complex.irp.f +++ b/src/ao_two_e_ints/map_integrals_complex.irp.f @@ -376,7 +376,9 @@ subroutine get_ao_two_e_integrals_periodic(j,k,l,sze,out_val) end -!subroutine get_ao_two_e_integrals_non_zero_periodic(j,k,l,sze,out_val,out_val_index,non_zero_int) +subroutine get_ao_two_e_integrals_non_zero_periodic(j,k,l,sze,out_val,out_val_index,non_zero_int) + print*,'not implemented for periodic',irp_here + stop -1 ! use map_module ! implicit none ! BEGIN_DOC @@ -418,8 +420,8 @@ end ! out_val_index(non_zero_int) = i ! out_val(non_zero_int) = tmp ! enddo -! -!end + +end !subroutine get_ao_two_e_integrals_non_zero_jl_periodic(j,l,thresh,sze_max,sze,out_val,out_val_index,non_zero_int) diff --git a/src/mo_basis/mos_complex.irp.f b/src/mo_basis/mos_complex.irp.f index 7a4361b7..75d3e169 100644 --- a/src/mo_basis/mos_complex.irp.f +++ b/src/mo_basis/mos_complex.irp.f @@ -166,7 +166,7 @@ subroutine ao_to_mo_complex(A_ao,LDA_ao,A_mo,LDA_mo) ! Transform A from the AO basis to the MO basis ! where A is complex in the AO basis ! - ! Ct.A_ao.C + ! C^\dagger.A_ao.C END_DOC integer, intent(in) :: LDA_ao,LDA_mo complex*16, intent(in) :: A_ao(LDA_ao,ao_num) @@ -189,6 +189,36 @@ subroutine ao_to_mo_complex(A_ao,LDA_ao,A_mo,LDA_mo) deallocate(T) end +subroutine ao_to_mo_noconjg_complex(A_ao,LDA_ao,A_mo,LDA_mo) + implicit none + BEGIN_DOC + ! Transform A from the AO basis to the MO basis + ! where A is complex in the AO basis + ! + ! C^T.A_ao.C + ! needed for 4idx tranform in four_idx_novvvv + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_ao(LDA_ao,ao_num) + complex*16, intent(out) :: A_mo(LDA_mo,mo_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num,mo_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call zgemm('N','N', ao_num, mo_num, ao_num, & + (1.d0,0.d0), A_ao,LDA_ao, & + mo_coef_complex, size(mo_coef_complex,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('T','N', mo_num, mo_num, ao_num, & + (1.d0,0.d0), mo_coef_complex,size(mo_coef_complex,1), & + T, ao_num, & + (0.d0,0.d0), A_mo, size(A_mo,1)) + + deallocate(T) +end + subroutine ao_ortho_cano_to_ao_complex(A_ao,LDA_ao,A,LDA) implicit none diff --git a/src/mo_two_e_ints/four_idx_novvvv_complex.irp.f b/src/mo_two_e_ints/four_idx_novvvv_complex.irp.f new file mode 100644 index 00000000..e02de3b7 --- /dev/null +++ b/src/mo_two_e_ints/four_idx_novvvv_complex.irp.f @@ -0,0 +1,247 @@ +BEGIN_PROVIDER [ complex*16, mo_coef_novirt_complex, (ao_num,n_core_inact_act_orb) ] + implicit none + BEGIN_DOC + ! MO coefficients without virtual MOs + END_DOC + integer :: j,jj + + do j=1,n_core_inact_act_orb + jj = list_core_inact_act(j) + mo_coef_novirt_complex(:,j) = mo_coef_complex(:,jj) + enddo + +END_PROVIDER + +subroutine ao_to_mo_novirt_complex(A_ao,LDA_ao,A_mo,LDA_mo) + implicit none + BEGIN_DOC + ! Transform A from the |AO| basis to the |MO| basis excluding virtuals + ! + ! $C^\dagger.A_{ao}.C$ + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_ao(LDA_ao,ao_num) + complex*16, intent(out) :: A_mo(LDA_mo,n_core_inact_act_orb) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num,n_core_inact_act_orb) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call zgemm('N','N', ao_num, n_core_inact_act_orb, ao_num, & + (1.d0,0.d0), A_ao,LDA_ao, & + mo_coef_novirt_complex, size(mo_coef_novirt_complex,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('C','N', n_core_inact_act_orb, n_core_inact_act_orb, ao_num,& + (1.d0,0.d0), mo_coef_novirt_complex,size(mo_coef_novirt_complex,1), & + T, ao_num, & + (0.d0,0.d0), A_mo, size(A_mo,1)) + + deallocate(T) +end + +subroutine ao_to_mo_novirt_conjg_complex(A_ao,LDA_ao,A_mo,LDA_mo) + implicit none + BEGIN_DOC + ! Transform A from the |AO| basis to the |MO| basis excluding virtuals + ! + ! $C^\dagger.A_{ao}.C^*$ + ! half-transformed ints as handled by four_idx_novvvv need to use this + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_ao(LDA_ao,ao_num) + complex*16, intent(out) :: A_mo(LDA_mo,n_core_inact_act_orb) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num,n_core_inact_act_orb) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call zgemm('N','N', ao_num, n_core_inact_act_orb, ao_num, & + (1.d0,0.d0), A_ao,LDA_ao, & + dconjg(mo_coef_novirt_complex), size(mo_coef_novirt_complex,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('C','N', n_core_inact_act_orb, n_core_inact_act_orb, ao_num,& + (1.d0,0.d0), mo_coef_novirt_complex,size(mo_coef_novirt_complex,1), & + T, ao_num, & + (0.d0,0.d0), A_mo, size(A_mo,1)) + + deallocate(T) +end + + +subroutine four_idx_novvvv_complex + use map_module + implicit none + BEGIN_DOC + ! Retransform MO integrals for next CAS-SCF step + END_DOC + integer :: i,j,k,l,n_integrals1,n_integrals2 + logical :: use_map1 + complex*16, allocatable :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:) + complex*16, external :: get_ao_two_e_integral_periodic + integer(key_kind), allocatable :: idx1(:),idx2(:) + complex(integral_kind), allocatable :: values1(:),values2(:) + double precision :: sign_tmp + integer(key_kind) :: idx_tmp + + integer :: p,q,r,s + allocate( T(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) , & + T2(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) ) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(mo_num,ao_num,T,n_core_inact_act_orb, & + !$OMP mo_integrals_threshold,mo_integrals_map, & + !$OMP mo_integrals_map_2,ao_integrals_map_2, & + !$OMP list_core_inact_act,T2,ao_integrals_map) & + !$OMP PRIVATE(i,j,k,l,p,q,r,s,idx1,idx2,values1,values2,n_integrals1, & + !$OMP n_integrals2,use_map1,idx_tmp,sign_tmp, & + !$OMP f,f2,d) + allocate(f(ao_num,ao_num,ao_num), f2(ao_num,ao_num,ao_num), d(mo_num,mo_num), & + idx1(2*mo_num*mo_num), values1(2*mo_num*mo_num), & + idx2(2*mo_num*mo_num), values2(2*mo_num*mo_num) ) + + ! + !$OMP DO + do s=1,ao_num + do r=1,ao_num + do q=1,ao_num + do p=1,r + f (p,q,r) = get_ao_two_e_integral_periodic(p,q,r,s,ao_integrals_map,ao_integrals_map_2) + f (r,q,p) = get_ao_two_e_integral_periodic(r,q,p,s,ao_integrals_map,ao_integrals_map_2) + enddo + enddo + enddo + do r=1,ao_num + do q=1,ao_num + do p=1,ao_num + f2(p,q,r) = f(p,r,q) + enddo + enddo + enddo + ! f (p,q,r) = + ! f2(p,q,r) = + + do r=1,ao_num + call ao_to_mo_novirt_conjg_complex(f (1,1,r),size(f ,1),T (1,1,r,s),size(T,1)) + call ao_to_mo_novirt_complex(f2(1,1,r),size(f2,1),T2(1,1,r,s),size(T,1)) + enddo + ! T (i,j,p,q) = + ! T2(i,j,p,q) = + + enddo + !$OMP END DO + + !$OMP DO + do j=1,n_core_inact_act_orb + do i=1,n_core_inact_act_orb + do s=1,ao_num + do r=1,ao_num + f (r,s,1) = T (i,j,r,s) + f2(r,s,1) = T2(i,j,r,s) + enddo + enddo + call ao_to_mo_noconjg_complex(f ,size(f ,1),d,size(d,1)) + n_integrals1 = 0 + n_integrals2 = 0 + do l=1,mo_num + do k=1,mo_num + call ao_two_e_integral_periodic_map_idx_sign(list_core_inact_act(i),list_core_inact_act(j),k,l,use_map1,idx_tmp,sign_tmp) + if (use_map1) then + n_integrals1+=1 + values1(n_integrals1) = dble(d(k,l)) + idx1(n_integrals1) = idx_tmp + if (sign_tmp /= 0.d0) then ! should always be true, but might change in the future + n_integrals1+=1 + values1(n_integrals1) = sign_tmp*dimag(d(k,l)) + idx1(n_integrals1) = idx_tmp+1 + endif + else + n_integrals2+=1 + values2(n_integrals2) = dble(d(k,l)) + idx2(n_integrals2) = idx_tmp + if (sign_tmp /= 0.d0) then + n_integrals2+=1 + values2(n_integrals2) = sign_tmp*dimag(d(k,l)) + idx2(n_integrals2) = idx_tmp+1 + endif + endif + enddo + enddo + call map_append(mo_integrals_map, idx1, values1, n_integrals1) + call map_append(mo_integrals_map_2, idx2, values2, n_integrals2) + + call ao_to_mo(f2,size(f2,1),d,size(d,1)) + n_integrals1 = 0 + n_integrals2 = 0 + do l=1,mo_num + do k=1,mo_num + call ao_two_e_integral_periodic_map_idx_sign(list_core_inact_act(i),k,list_core_inact_act(j),l,use_map1,idx_tmp,sign_tmp) + if (use_map1) then + n_integrals1+=1 + values1(n_integrals1) = dble(d(k,l)) + idx1(n_integrals1) = idx_tmp + if (sign_tmp /= 0.d0) then ! should always be true, but might change in the future + n_integrals1+=1 + values1(n_integrals1) = sign_tmp*dimag(d(k,l)) + idx1(n_integrals1) = idx_tmp+1 + endif + else + n_integrals2+=1 + values2(n_integrals2) = dble(d(k,l)) + idx2(n_integrals2) = idx_tmp + if (sign_tmp /= 0.d0) then + n_integrals2+=1 + values2(n_integrals2) = sign_tmp*dimag(d(k,l)) + idx2(n_integrals2) = idx_tmp+1 + endif + endif + enddo + enddo + call map_append(mo_integrals_map, idx1, values1, n_integrals1) + call map_append(mo_integrals_map_2, idx2, values2, n_integrals2) + enddo + enddo + !$OMP END DO + deallocate(f,f2,d,idx1,idx2,values1,values2) + + !$OMP END PARALLEL + + deallocate(T,T2) + + + call map_sort(mo_integrals_map) + call map_unique(mo_integrals_map) + call map_shrink(mo_integrals_map,real(mo_integrals_threshold,integral_kind)) + + call map_sort(mo_integrals_map_2) + call map_unique(mo_integrals_map_2) + call map_shrink(mo_integrals_map_2,real(mo_integrals_threshold,integral_kind)) + +end + +subroutine four_idx_novvvv2_complex + use bitmasks + implicit none + integer :: i + integer(bit_kind) :: mask_ijkl(N_int,4) + + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = full_ijkl_bitmask_4(i,1) + mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,4) = full_ijkl_bitmask_4(i,1) + enddo + call add_integrals_to_map_complex(mask_ijkl) + + print*, '' + do i = 1,N_int + mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1) + mask_ijkl(i,3) = virt_bitmask(i,1) + mask_ijkl(i,4) = virt_bitmask(i,1) + enddo + call add_integrals_to_map_complex(mask_ijkl) + +end diff --git a/src/mo_two_e_ints/map_integrals_complex.irp.f b/src/mo_two_e_ints/map_integrals_complex.irp.f index 1864d600..20970a15 100644 --- a/src/mo_two_e_ints/map_integrals_complex.irp.f +++ b/src/mo_two_e_ints/map_integrals_complex.irp.f @@ -306,7 +306,7 @@ subroutine get_mo_two_e_integrals_coulomb_ii_periodic(k,l,sze,out_val,map,map2) integer :: i integer(key_kind) :: hash(sze),hash_re(sze),hash_im(sze) real(integral_kind) :: tmp_re(sze),tmp_im(sze) - complex*16 :: out_re(sze),out_im(sze) + double precision :: out_re(sze),out_im(sze) double precision :: sign PROVIDE mo_two_e_integrals_in_map @@ -400,10 +400,10 @@ subroutine get_mo_two_e_integrals_exch_ii_periodic(k,l,sze,out_val,map,map2) integer, intent(in) :: k,l, sze double precision, intent(out) :: out_val(sze) type(map_type), intent(inout) :: map,map2 - integer :: i + integer :: i,klmin,klmax integer(key_kind) :: hash(sze),hash_re(sze),hash_im(sze) real(integral_kind) :: tmp_re(sze),tmp_im(sze) - complex*16 :: out_re(sze),out_im(sze) + double precision :: out_re(sze),out_im(sze) double precision :: sign,sign2(sze) PROVIDE mo_two_e_integrals_in_map diff --git a/src/mo_two_e_ints/mo_bi_integrals.irp.f b/src/mo_two_e_ints/mo_bi_integrals.irp.f index 11687602..bb998c26 100644 --- a/src/mo_two_e_ints/mo_bi_integrals.irp.f +++ b/src/mo_two_e_ints/mo_bi_integrals.irp.f @@ -56,11 +56,11 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ] if(no_vvvv_integrals)then print*,'not implemented for periodic',irp_here stop -1 - call four_idx_novvvv_periodic + call four_idx_novvvv_complex else print*,'not implemented for periodic',irp_here stop -1 - call add_integrals_to_map_periodic(full_ijkl_bitmask_4) + call add_integrals_to_map_complex(full_ijkl_bitmask_4) endif call wall_time(wall_2) @@ -981,13 +981,94 @@ end ! mo_two_e_integrals_jj_exchange_from_ao(i,j) = J_ij ! mo_two_e_integrals_jj_anti_from_ao(i,j) = J_ij - K_ij END_DOC - + integer :: i,j,p,q,r,s double precision :: c - real(integral_kind) :: integral integer :: n, pp - real(integral_kind), allocatable :: int_value(:) integer, allocatable :: int_idx(:) + if (is_periodic) then + complex(integral_kind) :: integral2 + complex(integral_kind), allocatable :: int_value2(:) + complex*16 :: cz + + complex*16, allocatable :: iqrs2(:,:), iqsr2(:,:), iqis2(:), iqri2(:) + PROVIDE ao_two_e_integrals_in_map mo_coef_complex + mo_two_e_integral_jj_from_ao = 0.d0 + mo_two_e_integrals_jj_exchange_from_ao = 0.d0 + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs2, iqsr2 + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i,j,p,q,r,s,integral2,c,n,pp,int_value2,int_idx, & + !$OMP iqrs2, iqsr2,iqri2,iqis2,cz) & + !$OMP SHARED(mo_num,mo_coef_transp_complex,mo_coef_transp_complex_conjg,ao_num, & + !$OMP ao_integrals_threshold) & + !$OMP REDUCTION(+:mo_two_e_integral_jj_from_ao,mo_two_e_integrals_jj_exchange_from_ao) + + allocate( int_value2(ao_num), int_idx(ao_num), & + iqrs2(mo_num,ao_num), iqis2(mo_num), iqri2(mo_num), & + iqsr2(mo_num,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=1,ao_num + + do j=1,ao_num + do i=1,mo_num + iqrs2(i,j) = (0.d0,0.d0) + iqsr2(i,j) = (0.d0,0.d0) + enddo + enddo + + + do r=1,ao_num + call get_ao_two_e_integrals_non_zero_periodic(q,r,s,ao_num,int_value2,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral2 = int_value2(pp) + if (cdabs(integral2) > ao_integrals_threshold) then + do i=1,mo_num + iqrs2(i,r) += mo_coef_transp_complex_conjg(i,p) * integral2 + enddo + endif + enddo + call get_ao_two_e_integrals_non_zero_periodic(q,s,r,ao_num,int_value2,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral2 = int_value2(pp) + if (cdabs(integral2) > ao_integrals_threshold) then + do i=1,mo_num + iqsr2(i,r) += mo_coef_transp_complex_conjg(i,p) * integral2 + enddo + endif + enddo + enddo + iqis2 = (0.d0,0.d0) + iqri2 = (0.d0,0.d0) + do r=1,ao_num + do i=1,mo_num + iqis2(i) += mo_coef_transp_complex(i,r) * iqrs2(i,r) + iqri2(i) += mo_coef_transp_complex(i,r) * iqsr2(i,r) + enddo + enddo + do i=1,mo_num + do j=1,mo_num + cz = mo_coef_transp_complex_conjg(j,q)*mo_coef_transp_complex(j,s) + mo_two_e_integral_jj_from_ao(j,i) += dble(cz * iqis2(i)) + mo_two_e_integrals_jj_exchange_from_ao(j,i) += dble(cz * iqri2(i)) + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + deallocate(iqrs2,iqsr2,int_value2,int_idx) + !$OMP END PARALLEL + + + else + real(integral_kind) :: integral + real(integral_kind), allocatable :: int_value(:) double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) @@ -1092,7 +1173,7 @@ end !$OMP END DO NOWAIT deallocate(iqrs,iqsr,int_value,int_idx) !$OMP END PARALLEL - + endif mo_two_e_integrals_jj_anti_from_ao = mo_two_e_integral_jj_from_ao - mo_two_e_integrals_jj_exchange_from_ao @@ -1112,11 +1193,100 @@ END_PROVIDER integer :: i,j,p,q,r,s integer :: i0,j0 double precision :: c - real(integral_kind) :: integral integer :: n, pp - real(integral_kind), allocatable :: int_value(:) integer, allocatable :: int_idx(:) + if (is_periodic) then + complex*16 :: cz + complex(integral_kind) :: integral2 + complex(integral_kind), allocatable :: int_value2(:) + complex*16, allocatable :: iqrs2(:,:), iqsr2(:,:), iqis2(:), iqri2(:) + + PROVIDE ao_two_e_integrals_in_map mo_coef_complex + + mo_two_e_integrals_vv_from_ao = 0.d0 + mo_two_e_integrals_vv_exchange_from_ao = 0.d0 + + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs2, iqsr2 + + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE (i0,j0,i,j,p,q,r,s,integral2,c,n,pp,int_value2,int_idx, & + !$OMP iqrs2, iqsr2,iqri2,iqis2,cz) & + !$OMP SHARED(n_virt_orb,mo_num,list_virt,mo_coef_transp_complex,ao_num, & + !$OMP mo_coef_transp_complex_conjg, & + !$OMP ao_integrals_threshold,do_direct_integrals) & + !$OMP REDUCTION(+:mo_two_e_integrals_vv_from_ao,mo_two_e_integrals_vv_exchange_from_ao) + + allocate( int_value2(ao_num), int_idx(ao_num), & + iqrs2(mo_num,ao_num), iqis2(mo_num), iqri2(mo_num),& + iqsr2(mo_num,ao_num) ) + + !$OMP DO SCHEDULE (guided) + do s=1,ao_num + do q=1,ao_num + + do j=1,ao_num + do i0=1,n_virt_orb + i = list_virt(i0) + iqrs2(i,j) = (0.d0,0.d0) + iqsr2(i,j) = (0.d0,0.d0) + enddo + enddo + + + do r=1,ao_num + call get_ao_two_e_integrals_non_zero_periodic(q,r,s,ao_num,int_value2,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral2 = int_value2(pp) + if (cdabs(integral2) > ao_integrals_threshold) then + do i0=1,n_virt_orb + i =list_virt(i0) + iqrs2(i,r) += mo_coef_transp_complex_conjg(i,p) * integral2 + enddo + endif + enddo + call get_ao_two_e_integrals_non_zero_periodic(q,s,r,ao_num,int_value2,int_idx,n) + do pp=1,n + p = int_idx(pp) + integral2 = int_value2(pp) + if (cdabs(integral2) > ao_integrals_threshold) then + do i0=1,n_virt_orb + i = list_virt(i0) + iqsr2(i,r) += mo_coef_transp_complex_conjg(i,p) * integral2 + enddo + endif + enddo + enddo + + iqis2 = (0.d0,0.d0) + iqri2 = (0.d0,0.d0) + do r=1,ao_num + do i0=1,n_virt_orb + i = list_virt(i0) + iqis2(i) += mo_coef_transp_complex(i,r) * iqrs2(i,r) + iqri2(i) += mo_coef_transp_complex(i,r) * iqsr2(i,r) + enddo + enddo + do i0=1,n_virt_orb + i= list_virt(i0) + do j0=1,n_virt_orb + j = list_virt(j0) + cz = mo_coef_transp_complex_conjg(j,q)*mo_coef_transp_complex(j,s) + mo_two_e_integrals_vv_from_ao(j,i) += dble(cz * iqis2(i)) + mo_two_e_integrals_vv_exchange_from_ao(j,i) += dble(cz * iqri2(i)) + enddo + enddo + + enddo + enddo + !$OMP END DO NOWAIT + deallocate(iqrs2,iqsr2,iqis2,iqri2,int_value2,int_idx) + !$OMP END PARALLEL + else + real(integral_kind) :: integral + real(integral_kind), allocatable :: int_value(:) double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) if (.not.do_direct_integrals) then @@ -1228,6 +1398,7 @@ END_PROVIDER !$OMP END DO NOWAIT deallocate(iqrs,iqsr,int_value,int_idx) !$OMP END PARALLEL + endif mo_two_e_integrals_vv_anti_from_ao = mo_two_e_integrals_vv_from_ao - mo_two_e_integrals_vv_exchange_from_ao ! print*, '**********' @@ -1257,7 +1428,18 @@ END_PROVIDER PROVIDE mo_two_e_integrals_in_map mo_two_e_integrals_jj = 0.d0 mo_two_e_integrals_jj_exchange = 0.d0 - + if (is_periodic) then + complex*16 :: get_two_e_integral_periodic + do j=1,mo_num + do i=1,mo_num + mo_two_e_integrals_jj(i,j) = dble(get_two_e_integral_periodic(i,j,i,j,& + mo_integrals_map,mo_integrals_map_2)) + mo_two_e_integrals_jj_exchange(i,j) = dble(get_two_e_integral_periodic(i,j,j,i,& + mo_integrals_map,mo_integrals_map_2)) + mo_two_e_integrals_jj_anti(i,j) = mo_two_e_integrals_jj(i,j) - mo_two_e_integrals_jj_exchange(i,j) + enddo + enddo + else do j=1,mo_num do i=1,mo_num mo_two_e_integrals_jj(i,j) = get_two_e_integral(i,j,i,j,mo_integrals_map) @@ -1265,6 +1447,7 @@ END_PROVIDER mo_two_e_integrals_jj_anti(i,j) = mo_two_e_integrals_jj(i,j) - mo_two_e_integrals_jj_exchange(i,j) enddo enddo + endif END_PROVIDER @@ -1275,6 +1458,9 @@ subroutine clear_mo_map ! Frees the memory of the MO map END_DOC call map_deinit(mo_integrals_map) + if (is_periodic) then + call map_deinit(mo_integrals_map_2) + endif FREE mo_integrals_map mo_two_e_integrals_jj mo_two_e_integrals_jj_anti FREE mo_two_e_integrals_jj_exchange mo_two_e_integrals_in_map end diff --git a/src/mo_two_e_ints/mo_bi_integrals_complex.irp.f b/src/mo_two_e_ints/mo_bi_integrals_complex.irp.f new file mode 100644 index 00000000..8d1dd9f0 --- /dev/null +++ b/src/mo_two_e_ints/mo_bi_integrals_complex.irp.f @@ -0,0 +1,1163 @@ + + +subroutine add_integrals_to_map_complex(mask_ijkl) + use map_module + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: two_e_tmp_0_idx(:) + real(integral_kind), allocatable :: two_e_tmp_0(:,:) + double precision, allocatable :: two_e_tmp_1(:) + double precision, allocatable :: two_e_tmp_2(:,:) + double precision, allocatable :: two_e_tmp_3(:,:,:) + !DIR$ ATTRIBUTES ALIGN : 64 :: two_e_tmp_1, two_e_tmp_2, two_e_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + double precision, external :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + print*,'not implemented for complex',irp_here + stop -1 +! PROVIDE ao_two_e_integrals_in_map mo_coef +! +! !Get list of MOs for i,j,k and l +! !------------------------------- +! +! allocate(list_ijkl(mo_num,4)) +! call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) +! call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) +! call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) +! call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijkl(i,1)) +! enddo +! if(j==0)then +! return +! endif +! +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijkl(i,2)) +! enddo +! if(j==0)then +! return +! endif +! +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijkl(i,3)) +! enddo +! if(j==0)then +! return +! endif +! +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijkl(i,4)) +! enddo +! if(j==0)then +! return +! endif +! +! size_buffer = min(ao_num*ao_num*ao_num,16000000) +! print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+& +! ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' +! +! double precision :: accu_bis +! accu_bis = 0.d0 +! call wall_time(wall_1) +! +! !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & +! !$OMP two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,& +! !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & +! !$OMP wall_0,thread_num,accu_bis) & +! !$OMP DEFAULT(NONE) & +! !$OMP SHARED(size_buffer,ao_num,mo_num,n_i,n_j,n_k,n_l, & +! !$OMP mo_coef_transp, & +! !$OMP mo_coef_transp_is_built, list_ijkl, & +! !$OMP mo_coef_is_built, wall_1, & +! !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) +! n_integrals = 0 +! wall_0 = wall_1 +! allocate(two_e_tmp_3(mo_num, n_j, n_k), & +! two_e_tmp_1(mo_num), & +! two_e_tmp_0(ao_num,ao_num), & +! two_e_tmp_0_idx(ao_num), & +! two_e_tmp_2(mo_num, n_j), & +! buffer_i(size_buffer), & +! buffer_value(size_buffer) ) +! +! thread_num = 0 +! !$ thread_num = omp_get_thread_num() +! !$OMP DO SCHEDULE(guided) +! do l1 = 1,ao_num +! two_e_tmp_3 = 0.d0 +! do k1 = 1,ao_num +! two_e_tmp_2 = 0.d0 +! do j1 = 1,ao_num +! call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) +! ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) +! enddo +! do j1 = 1,ao_num +! kmax = 0 +! do i1 = 1,ao_num +! c = two_e_tmp_0(i1,j1) +! if (c == 0.d0) then +! cycle +! endif +! kmax += 1 +! two_e_tmp_0(kmax,j1) = c +! two_e_tmp_0_idx(kmax) = i1 +! enddo +! +! if (kmax==0) then +! cycle +! endif +! +! two_e_tmp_1 = 0.d0 +! ii1=1 +! do ii1 = 1,kmax-4,4 +! i1 = two_e_tmp_0_idx(ii1) +! i2 = two_e_tmp_0_idx(ii1+1) +! i3 = two_e_tmp_0_idx(ii1+2) +! i4 = two_e_tmp_0_idx(ii1+3) +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_1(i) = two_e_tmp_1(i) + & +! mo_coef_transp(i,i1) * two_e_tmp_0(ii1,j1) + & +! mo_coef_transp(i,i2) * two_e_tmp_0(ii1+1,j1) + & +! mo_coef_transp(i,i3) * two_e_tmp_0(ii1+2,j1) + & +! mo_coef_transp(i,i4) * two_e_tmp_0(ii1+3,j1) +! enddo ! i +! enddo ! ii1 +! +! i2 = ii1 +! do ii1 = i2,kmax +! i1 = two_e_tmp_0_idx(ii1) +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_1(i) = two_e_tmp_1(i) + mo_coef_transp(i,i1) * two_e_tmp_0(ii1,j1) +! enddo ! i +! enddo ! ii1 +! c = 0.d0 +! +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! c = max(c,abs(two_e_tmp_1(i))) +! if (c>mo_integrals_threshold) exit +! enddo +! if ( c < mo_integrals_threshold ) then +! cycle +! endif +! +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! c = mo_coef_transp(j,j1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_2(i,j0) = two_e_tmp_2(i,j0) + c * two_e_tmp_1(i) +! enddo ! i +! enddo ! j +! enddo !j1 +! if ( maxval(abs(two_e_tmp_2)) < mo_integrals_threshold ) then +! cycle +! endif +! +! +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! c = mo_coef_transp(k,k1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! do i = list_ijkl(1,1), k +! two_e_tmp_3(i,j0,k0) = two_e_tmp_3(i,j0,k0) + c* two_e_tmp_2(i,j0) +! enddo!i +! enddo !j +! +! enddo !k +! enddo !k1 +! +! +! +! do l0 = 1,n_l +! l = list_ijkl(l0,4) +! c = mo_coef_transp(l,l1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! j1 = shiftr((l*l-l),1) +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! if (j > l) then +! exit +! endif +! j1 += 1 +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! i1 = shiftr((k*k-k),1) +! if (i1<=j1) then +! continue +! else +! exit +! endif +! two_e_tmp_1 = 0.d0 +! do i0 = 1, n_i +! i = list_ijkl(i0,1) +! if (i>k) then +! exit +! endif +! two_e_tmp_1(i) = c*two_e_tmp_3(i,j0,k0) +! ! i1+=1 +! enddo +! +! do i0 = 1, n_i +! i = list_ijkl(i0,1) +! if(i> min(k,j1-i1+list_ijkl(1,1)-1))then +! exit +! endif +! if (abs(two_e_tmp_1(i)) < mo_integrals_threshold) then +! cycle +! endif +! n_integrals += 1 +! buffer_value(n_integrals) = two_e_tmp_1(i) +! !DIR$ FORCEINLINE +! call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) +! if (n_integrals == size_buffer) then +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! n_integrals = 0 +! endif +! enddo +! enddo +! enddo +! enddo +! +! call wall_time(wall_2) +! if (thread_num == 0) then +! if (wall_2 - wall_0 > 1.d0) then +! wall_0 = wall_2 +! print*, 100.*float(l1)/float(ao_num), '% in ', & +! wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' +! endif +! endif +! enddo +! !$OMP END DO NOWAIT +! deallocate (two_e_tmp_1,two_e_tmp_2,two_e_tmp_3) +! +! integer :: index_needed +! +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! deallocate(buffer_i, buffer_value) +! !$OMP END PARALLEL +! call map_merge(mo_integrals_map) +! +! call wall_time(wall_2) +! call cpu_time(cpu_2) +! integer*8 :: get_mo_map_size, mo_map_size +! mo_map_size = get_mo_map_size() +! +! deallocate(list_ijkl) + + +end + + +subroutine add_integrals_to_map_three_indices_complex(mask_ijk) + use map_module + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijk(N_int,3) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k + integer :: m + integer, allocatable :: two_e_tmp_0_idx(:) + real(integral_kind), allocatable :: two_e_tmp_0(:,:) + double precision, allocatable :: two_e_tmp_1(:) + double precision, allocatable :: two_e_tmp_2(:,:) + double precision, allocatable :: two_e_tmp_3(:,:,:) + !DIR$ ATTRIBUTES ALIGN : 64 :: two_e_tmp_1, two_e_tmp_2, two_e_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + double precision :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + print*,'not implemented for complex',irp_here + stop -1 +! PROVIDE ao_two_e_integrals_in_map mo_coef +! +! !Get list of MOs for i,j,k and l +! !------------------------------- +! +! allocate(list_ijkl(mo_num,4)) +! call bitstring_to_list( mask_ijk(1,1), list_ijkl(1,1), n_i, N_int ) +! call bitstring_to_list( mask_ijk(1,2), list_ijkl(1,2), n_j, N_int ) +! call bitstring_to_list( mask_ijk(1,3), list_ijkl(1,3), n_k, N_int ) +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijk(i,1)) +! enddo +! if(j==0)then +! return +! endif +! +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijk(i,2)) +! enddo +! if(j==0)then +! return +! endif +! +! j = 0 +! do i = 1, N_int +! j += popcnt(mask_ijk(i,3)) +! enddo +! if(j==0)then +! return +! endif +! +! size_buffer = min(ao_num*ao_num*ao_num,16000000) +! print*, 'Providing the molecular integrals ' +! print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+& +! ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' +! +! call wall_time(wall_1) +! call cpu_time(cpu_1) +! double precision :: accu_bis +! accu_bis = 0.d0 +! !$OMP PARALLEL PRIVATE(m,l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & +! !$OMP two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,& +! !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & +! !$OMP wall_0,thread_num,accu_bis) & +! !$OMP DEFAULT(NONE) & +! !$OMP SHARED(size_buffer,ao_num,mo_num,n_i,n_j,n_k, & +! !$OMP mo_coef_transp, & +! !$OMP mo_coef_transp_is_built, list_ijkl, & +! !$OMP mo_coef_is_built, wall_1, & +! !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) +! n_integrals = 0 +! wall_0 = wall_1 +! allocate(two_e_tmp_3(mo_num, n_j, n_k), & +! two_e_tmp_1(mo_num), & +! two_e_tmp_0(ao_num,ao_num), & +! two_e_tmp_0_idx(ao_num), & +! two_e_tmp_2(mo_num, n_j), & +! buffer_i(size_buffer), & +! buffer_value(size_buffer) ) +! +! thread_num = 0 +! !$ thread_num = omp_get_thread_num() +! !$OMP DO SCHEDULE(guided) +! do l1 = 1,ao_num +! two_e_tmp_3 = 0.d0 +! do k1 = 1,ao_num +! two_e_tmp_2 = 0.d0 +! do j1 = 1,ao_num +! call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) +! enddo +! do j1 = 1,ao_num +! kmax = 0 +! do i1 = 1,ao_num +! c = two_e_tmp_0(i1,j1) +! if (c == 0.d0) then +! cycle +! endif +! kmax += 1 +! two_e_tmp_0(kmax,j1) = c +! two_e_tmp_0_idx(kmax) = i1 +! enddo +! +! if (kmax==0) then +! cycle +! endif +! +! two_e_tmp_1 = 0.d0 +! ii1=1 +! do ii1 = 1,kmax-4,4 +! i1 = two_e_tmp_0_idx(ii1) +! i2 = two_e_tmp_0_idx(ii1+1) +! i3 = two_e_tmp_0_idx(ii1+2) +! i4 = two_e_tmp_0_idx(ii1+3) +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_1(i) = two_e_tmp_1(i) + & +! mo_coef_transp(i,i1) * two_e_tmp_0(ii1,j1) + & +! mo_coef_transp(i,i2) * two_e_tmp_0(ii1+1,j1) + & +! mo_coef_transp(i,i3) * two_e_tmp_0(ii1+2,j1) + & +! mo_coef_transp(i,i4) * two_e_tmp_0(ii1+3,j1) +! enddo ! i +! enddo ! ii1 +! +! i2 = ii1 +! do ii1 = i2,kmax +! i1 = two_e_tmp_0_idx(ii1) +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_1(i) = two_e_tmp_1(i) + mo_coef_transp(i,i1) * two_e_tmp_0(ii1,j1) +! enddo ! i +! enddo ! ii1 +! c = 0.d0 +! +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! c = max(c,abs(two_e_tmp_1(i))) +! if (c>mo_integrals_threshold) exit +! enddo +! if ( c < mo_integrals_threshold ) then +! cycle +! endif +! +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! c = mo_coef_transp(j,j1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_2(i,j0) = two_e_tmp_2(i,j0) + c * two_e_tmp_1(i) +! enddo ! i +! enddo ! j +! enddo !j1 +! if ( maxval(abs(two_e_tmp_2)) < mo_integrals_threshold ) then +! cycle +! endif +! +! +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! c = mo_coef_transp(k,k1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! do i = list_ijkl(1,1), k +! two_e_tmp_3(i,j0,k0) = two_e_tmp_3(i,j0,k0) + c* two_e_tmp_2(i,j0) +! enddo!i +! enddo !j +! +! enddo !k +! enddo !k1 +! +! +! +! do l0 = 1,n_j +! l = list_ijkl(l0,2) +! c = mo_coef_transp(l,l1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! i1 = shiftr((k*k-k),1) +! two_e_tmp_1 = 0.d0 +! j0 = l0 +! j = list_ijkl(j0,2) +! do i0 = 1, n_i +! i = list_ijkl(i0,1) +! if (i>k) then +! exit +! endif +! two_e_tmp_1(i) = c*two_e_tmp_3(i,j0,k0) +! enddo +! +! do i0 = 1, n_i +! i = list_ijkl(i0,1) +! if (i>k) then !min(k,j1-i1) +! exit +! endif +! if (abs(two_e_tmp_1(i)) < mo_integrals_threshold) then +! cycle +! endif +! n_integrals += 1 +! buffer_value(n_integrals) = two_e_tmp_1(i) +! if(i==k .and. j==l .and. i.ne.j)then +! buffer_value(n_integrals) = buffer_value(n_integrals) *0.5d0 +! endif +! !DIR$ FORCEINLINE +! call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) +! if (n_integrals == size_buffer) then +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! n_integrals = 0 +! endif +! enddo +! enddo +! enddo +! +! do l0 = 1,n_j +! l = list_ijkl(l0,2) +! c = mo_coef_transp(l,l1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! i1 = shiftr((k*k-k),1) +! two_e_tmp_1 = 0.d0 +! j0 = k0 +! j = list_ijkl(k0,2) +! i0 = l0 +! i = list_ijkl(i0,2) +! if (k==l) then +! cycle +! endif +! two_e_tmp_1(i) = c*two_e_tmp_3(i,j0,k0) +! +! n_integrals += 1 +! buffer_value(n_integrals) = two_e_tmp_1(i) +! !DIR$ FORCEINLINE +! call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) +! if (n_integrals == size_buffer) then +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! n_integrals = 0 +! endif +! enddo +! enddo +! +! call wall_time(wall_2) +! if (thread_num == 0) then +! if (wall_2 - wall_0 > 1.d0) then +! wall_0 = wall_2 +! print*, 100.*float(l1)/float(ao_num), '% in ', & +! wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' +! endif +! endif +! enddo +! !$OMP END DO NOWAIT +! deallocate (two_e_tmp_1,two_e_tmp_2,two_e_tmp_3) +! +! integer :: index_needed +! +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! deallocate(buffer_i, buffer_value) +! !$OMP END PARALLEL +! call map_merge(mo_integrals_map) +! +! call wall_time(wall_2) +! call cpu_time(cpu_2) +! integer*8 :: get_mo_map_size, mo_map_size +! mo_map_size = get_mo_map_size() +! +! deallocate(list_ijkl) +! +! +! print*,'Molecular integrals provided:' +! print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' +! print*,' Number of MO integrals: ', mo_map_size +! print*,' cpu time :',cpu_2 - cpu_1, 's' +! print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + +end + + +subroutine add_integrals_to_map_no_exit_34_complex(mask_ijkl) + use map_module + use bitmasks + implicit none + + BEGIN_DOC + ! Adds integrals to tha MO map according to some bitmask + END_DOC + + integer(bit_kind), intent(in) :: mask_ijkl(N_int,4) + + integer :: i,j,k,l + integer :: i0,j0,k0,l0 + double precision :: c, cpu_1, cpu_2, wall_1, wall_2, wall_0 + + integer, allocatable :: list_ijkl(:,:) + integer :: n_i, n_j, n_k, n_l + integer, allocatable :: two_e_tmp_0_idx(:) + real(integral_kind), allocatable :: two_e_tmp_0(:,:) + double precision, allocatable :: two_e_tmp_1(:) + double precision, allocatable :: two_e_tmp_2(:,:) + double precision, allocatable :: two_e_tmp_3(:,:,:) + !DIR$ ATTRIBUTES ALIGN : 64 :: two_e_tmp_1, two_e_tmp_2, two_e_tmp_3 + + integer :: n_integrals + integer :: size_buffer + integer(key_kind),allocatable :: buffer_i(:) + real(integral_kind),allocatable :: buffer_value(:) + double precision :: map_mb + + integer :: i1,j1,k1,l1, ii1, kmax, thread_num + integer :: i2,i3,i4 + double precision,parameter :: thr_coef = 1.d-10 + + print*,'not implemented for complex',irp_here + stop -1 +! PROVIDE ao_two_e_integrals_in_map mo_coef +! +! !Get list of MOs for i,j,k and l +! !------------------------------- +! +! allocate(list_ijkl(mo_num,4)) +! call bitstring_to_list( mask_ijkl(1,1), list_ijkl(1,1), n_i, N_int ) +! call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int ) +! call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int ) +! call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int ) +! +! size_buffer = min(ao_num*ao_num*ao_num,16000000) +! print*, 'Providing the molecular integrals ' +! print*, 'Buffers : ', 8.*(mo_num*(n_j)*(n_k+1) + mo_num+& +! ao_num+ao_num*ao_num+ size_buffer*3)/(1024*1024), 'MB / core' +! +! call wall_time(wall_1) +! call cpu_time(cpu_1) +! +! !$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, & +! !$OMP two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,& +! !$OMP buffer_i,buffer_value,n_integrals,wall_2,i0,j0,k0,l0, & +! !$OMP wall_0,thread_num) & +! !$OMP DEFAULT(NONE) & +! !$OMP SHARED(size_buffer,ao_num,mo_num,n_i,n_j,n_k,n_l, & +! !$OMP mo_coef_transp, & +! !$OMP mo_coef_transp_is_built, list_ijkl, & +! !$OMP mo_coef_is_built, wall_1, & +! !$OMP mo_coef,mo_integrals_threshold,mo_integrals_map) +! n_integrals = 0 +! wall_0 = wall_1 +! allocate(two_e_tmp_3(mo_num, n_j, n_k), & +! two_e_tmp_1(mo_num), & +! two_e_tmp_0(ao_num,ao_num), & +! two_e_tmp_0_idx(ao_num), & +! two_e_tmp_2(mo_num, n_j), & +! buffer_i(size_buffer), & +! buffer_value(size_buffer) ) +! +! thread_num = 0 +! !$ thread_num = omp_get_thread_num() +! !$OMP DO SCHEDULE(guided) +! do l1 = 1,ao_num +! !IRP_IF COARRAY +! ! if (mod(l1-this_image(),num_images()) /= 0 ) then +! ! cycle +! ! endif +! !IRP_ENDIF +! two_e_tmp_3 = 0.d0 +! do k1 = 1,ao_num +! two_e_tmp_2 = 0.d0 +! do j1 = 1,ao_num +! call get_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) +! ! call compute_ao_two_e_integrals(j1,k1,l1,ao_num,two_e_tmp_0(1,j1)) +! enddo +! do j1 = 1,ao_num +! kmax = 0 +! do i1 = 1,ao_num +! c = two_e_tmp_0(i1,j1) +! if (c == 0.d0) then +! cycle +! endif +! kmax += 1 +! two_e_tmp_0(kmax,j1) = c +! two_e_tmp_0_idx(kmax) = i1 +! enddo +! +! if (kmax==0) then +! cycle +! endif +! +! two_e_tmp_1 = 0.d0 +! ii1=1 +! do ii1 = 1,kmax-4,4 +! i1 = two_e_tmp_0_idx(ii1) +! i2 = two_e_tmp_0_idx(ii1+1) +! i3 = two_e_tmp_0_idx(ii1+2) +! i4 = two_e_tmp_0_idx(ii1+3) +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_1(i) = two_e_tmp_1(i) + & +! mo_coef_transp(i,i1) * two_e_tmp_0(ii1,j1) + & +! mo_coef_transp(i,i2) * two_e_tmp_0(ii1+1,j1) + & +! mo_coef_transp(i,i3) * two_e_tmp_0(ii1+2,j1) + & +! mo_coef_transp(i,i4) * two_e_tmp_0(ii1+3,j1) +! enddo ! i +! enddo ! ii1 +! +! i2 = ii1 +! do ii1 = i2,kmax +! i1 = two_e_tmp_0_idx(ii1) +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_1(i) = two_e_tmp_1(i) + mo_coef_transp(i,i1) * two_e_tmp_0(ii1,j1) +! enddo ! i +! enddo ! ii1 +! c = 0.d0 +! +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! c = max(c,abs(two_e_tmp_1(i))) +! if (c>mo_integrals_threshold) exit +! enddo +! if ( c < mo_integrals_threshold ) then +! cycle +! endif +! +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! c = mo_coef_transp(j,j1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! do i = list_ijkl(1,1), list_ijkl(n_i,1) +! two_e_tmp_2(i,j0) = two_e_tmp_2(i,j0) + c * two_e_tmp_1(i) +! enddo ! i +! enddo ! j +! enddo !j1 +! if ( maxval(abs(two_e_tmp_2)) < mo_integrals_threshold ) then +! cycle +! endif +! +! +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! c = mo_coef_transp(k,k1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! do i = list_ijkl(1,1), k +! two_e_tmp_3(i,j0,k0) = two_e_tmp_3(i,j0,k0) + c* two_e_tmp_2(i,j0) +! enddo!i +! enddo !j +! +! enddo !k +! enddo !k1 +! +! +! +! do l0 = 1,n_l +! l = list_ijkl(l0,4) +! c = mo_coef_transp(l,l1) +! if (abs(c) < thr_coef) then +! cycle +! endif +! j1 = shiftr((l*l-l),1) +! do j0 = 1, n_j +! j = list_ijkl(j0,2) +! if (j > l) then +! exit +! endif +! j1 += 1 +! do k0 = 1, n_k +! k = list_ijkl(k0,3) +! i1 = shiftr((k*k-k),1) +! two_e_tmp_1 = 0.d0 +! do i0 = 1, n_i +! i = list_ijkl(i0,1) +! if (i>k) then +! exit +! endif +! two_e_tmp_1(i) = c*two_e_tmp_3(i,j0,k0) +! enddo +! +! do i0 = 1, n_i +! i = list_ijkl(i0,1) +! if(i> k)then +! exit +! endif +! +! if (abs(two_e_tmp_1(i)) < mo_integrals_threshold) then +! cycle +! endif +! n_integrals += 1 +! buffer_value(n_integrals) = two_e_tmp_1(i) +! !DIR$ FORCEINLINE +! call mo_two_e_integrals_index(i,j,k,l,buffer_i(n_integrals)) +! if (n_integrals == size_buffer) then +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! n_integrals = 0 +! endif +! enddo +! enddo +! enddo +! enddo +! +! call wall_time(wall_2) +! if (thread_num == 0) then +! if (wall_2 - wall_0 > 1.d0) then +! wall_0 = wall_2 +! print*, 100.*float(l1)/float(ao_num), '% in ', & +! wall_2-wall_1, 's', map_mb(mo_integrals_map) ,'MB' +! endif +! endif +! enddo +! !$OMP END DO NOWAIT +! deallocate (two_e_tmp_1,two_e_tmp_2,two_e_tmp_3) +! +! call insert_into_mo_integrals_map(n_integrals,buffer_i,buffer_value,& +! real(mo_integrals_threshold,integral_kind)) +! deallocate(buffer_i, buffer_value) +! !$OMP END PARALLEL +! !IRP_IF COARRAY +! ! print*, 'Communicating the map' +! ! call communicate_mo_integrals() +! !IRP_ENDIF +! call map_merge(mo_integrals_map) +! +! call wall_time(wall_2) +! call cpu_time(cpu_2) +! integer*8 :: get_mo_map_size, mo_map_size +! mo_map_size = get_mo_map_size() +! +! deallocate(list_ijkl) +! +! +! print*,'Molecular integrals provided:' +! print*,' Size of MO map ', map_mb(mo_integrals_map) ,'MB' +! print*,' Number of MO integrals: ', mo_map_size +! print*,' cpu time :',cpu_2 - cpu_1, 's' +! print*,' wall time :',wall_2 - wall_1, 's ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')' + + +end + + + +! BEGIN_PROVIDER [ double precision, mo_two_e_integral_jj_from_ao, (mo_num,mo_num) ] +!&BEGIN_PROVIDER [ double precision, mo_two_e_integrals_jj_exchange_from_ao, (mo_num,mo_num) ] +!&BEGIN_PROVIDER [ double precision, mo_two_e_integrals_jj_anti_from_ao, (mo_num,mo_num) ] +! implicit none +! BEGIN_DOC +! ! mo_two_e_integral_jj_from_ao(i,j) = J_ij +! ! mo_two_e_integrals_jj_exchange_from_ao(i,j) = J_ij +! ! mo_two_e_integrals_jj_anti_from_ao(i,j) = J_ij - K_ij +! END_DOC +! +! integer :: i,j,p,q,r,s +! double precision :: c +! real(integral_kind) :: integral +! integer :: n, pp +! real(integral_kind), allocatable :: int_value(:) +! integer, allocatable :: int_idx(:) +! +! double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) +! +! if (.not.do_direct_integrals) then +! PROVIDE ao_two_e_integrals_in_map mo_coef +! endif +! +! mo_two_e_integral_jj_from_ao = 0.d0 +! mo_two_e_integrals_jj_exchange_from_ao = 0.d0 +! +! !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr +! +! +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE (i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, & +! !$OMP iqrs, iqsr,iqri,iqis) & +! !$OMP SHARED(mo_num,mo_coef_transp,ao_num, & +! !$OMP ao_integrals_threshold,do_direct_integrals) & +! !$OMP REDUCTION(+:mo_two_e_integral_jj_from_ao,mo_two_e_integrals_jj_exchange_from_ao) +! +! allocate( int_value(ao_num), int_idx(ao_num), & +! iqrs(mo_num,ao_num), iqis(mo_num), iqri(mo_num), & +! iqsr(mo_num,ao_num) ) +! +! !$OMP DO SCHEDULE (guided) +! do s=1,ao_num +! do q=1,ao_num +! +! do j=1,ao_num +! do i=1,mo_num +! iqrs(i,j) = 0.d0 +! iqsr(i,j) = 0.d0 +! enddo +! enddo +! +! if (do_direct_integrals) then +! double precision :: ao_two_e_integral +! do r=1,ao_num +! call compute_ao_two_e_integrals(q,r,s,ao_num,int_value) +! do p=1,ao_num +! integral = int_value(p) +! if (abs(integral) > ao_integrals_threshold) then +! do i=1,mo_num +! iqrs(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! call compute_ao_two_e_integrals(q,s,r,ao_num,int_value) +! do p=1,ao_num +! integral = int_value(p) +! if (abs(integral) > ao_integrals_threshold) then +! do i=1,mo_num +! iqsr(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! enddo +! +! else +! +! do r=1,ao_num +! call get_ao_two_e_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n) +! do pp=1,n +! p = int_idx(pp) +! integral = int_value(pp) +! if (abs(integral) > ao_integrals_threshold) then +! do i=1,mo_num +! iqrs(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! call get_ao_two_e_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n) +! do pp=1,n +! p = int_idx(pp) +! integral = int_value(pp) +! if (abs(integral) > ao_integrals_threshold) then +! do i=1,mo_num +! iqsr(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! enddo +! endif +! iqis = 0.d0 +! iqri = 0.d0 +! do r=1,ao_num +! do i=1,mo_num +! iqis(i) += mo_coef_transp(i,r) * iqrs(i,r) +! iqri(i) += mo_coef_transp(i,r) * iqsr(i,r) +! enddo +! enddo +! do i=1,mo_num +! do j=1,mo_num +! c = mo_coef_transp(j,q)*mo_coef_transp(j,s) +! mo_two_e_integral_jj_from_ao(j,i) += c * iqis(i) +! mo_two_e_integrals_jj_exchange_from_ao(j,i) += c * iqri(i) +! enddo +! enddo +! +! enddo +! enddo +! !$OMP END DO NOWAIT +! deallocate(iqrs,iqsr,int_value,int_idx) +! !$OMP END PARALLEL +! +! mo_two_e_integrals_jj_anti_from_ao = mo_two_e_integral_jj_from_ao - mo_two_e_integrals_jj_exchange_from_ao +! +! +!END_PROVIDER +! +! BEGIN_PROVIDER [ double precision, mo_two_e_integrals_vv_from_ao, (mo_num,mo_num) ] +!&BEGIN_PROVIDER [ double precision, mo_two_e_integrals_vv_exchange_from_ao, (mo_num,mo_num) ] +!&BEGIN_PROVIDER [ double precision, mo_two_e_integrals_vv_anti_from_ao, (mo_num,mo_num) ] +! implicit none +! BEGIN_DOC +! ! mo_two_e_integrals_vv_from_ao(i,j) = J_ij +! ! mo_two_e_integrals_vv_exchange_from_ao(i,j) = J_ij +! ! mo_two_e_integrals_vv_anti_from_ao(i,j) = J_ij - K_ij +! ! but only for the virtual orbitals +! END_DOC +! +! integer :: i,j,p,q,r,s +! integer :: i0,j0 +! double precision :: c +! real(integral_kind) :: integral +! integer :: n, pp +! real(integral_kind), allocatable :: int_value(:) +! integer, allocatable :: int_idx(:) +! +! double precision, allocatable :: iqrs(:,:), iqsr(:,:), iqis(:), iqri(:) +! +! if (.not.do_direct_integrals) then +! PROVIDE ao_two_e_integrals_in_map mo_coef +! endif +! +! mo_two_e_integrals_vv_from_ao = 0.d0 +! mo_two_e_integrals_vv_exchange_from_ao = 0.d0 +! +! !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: iqrs, iqsr +! +! +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE (i0,j0,i,j,p,q,r,s,integral,c,n,pp,int_value,int_idx, & +! !$OMP iqrs, iqsr,iqri,iqis) & +! !$OMP SHARED(n_virt_orb,mo_num,list_virt,mo_coef_transp,ao_num, & +! !$OMP ao_integrals_threshold,do_direct_integrals) & +! !$OMP REDUCTION(+:mo_two_e_integrals_vv_from_ao,mo_two_e_integrals_vv_exchange_from_ao) +! +! allocate( int_value(ao_num), int_idx(ao_num), & +! iqrs(mo_num,ao_num), iqis(mo_num), iqri(mo_num),& +! iqsr(mo_num,ao_num) ) +! +! !$OMP DO SCHEDULE (guided) +! do s=1,ao_num +! do q=1,ao_num +! +! do j=1,ao_num +! do i0=1,n_virt_orb +! i = list_virt(i0) +! iqrs(i,j) = 0.d0 +! iqsr(i,j) = 0.d0 +! enddo +! enddo +! +! if (do_direct_integrals) then +! double precision :: ao_two_e_integral +! do r=1,ao_num +! call compute_ao_two_e_integrals(q,r,s,ao_num,int_value) +! do p=1,ao_num +! integral = int_value(p) +! if (abs(integral) > ao_integrals_threshold) then +! do i0=1,n_virt_orb +! i = list_virt(i0) +! iqrs(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! call compute_ao_two_e_integrals(q,s,r,ao_num,int_value) +! do p=1,ao_num +! integral = int_value(p) +! if (abs(integral) > ao_integrals_threshold) then +! do i0=1,n_virt_orb +! i =list_virt(i0) +! iqsr(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! enddo +! +! else +! +! do r=1,ao_num +! call get_ao_two_e_integrals_non_zero(q,r,s,ao_num,int_value,int_idx,n) +! do pp=1,n +! p = int_idx(pp) +! integral = int_value(pp) +! if (abs(integral) > ao_integrals_threshold) then +! do i0=1,n_virt_orb +! i =list_virt(i0) +! iqrs(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! call get_ao_two_e_integrals_non_zero(q,s,r,ao_num,int_value,int_idx,n) +! do pp=1,n +! p = int_idx(pp) +! integral = int_value(pp) +! if (abs(integral) > ao_integrals_threshold) then +! do i0=1,n_virt_orb +! i = list_virt(i0) +! iqsr(i,r) += mo_coef_transp(i,p) * integral +! enddo +! endif +! enddo +! enddo +! endif +! iqis = 0.d0 +! iqri = 0.d0 +! do r=1,ao_num +! do i0=1,n_virt_orb +! i = list_virt(i0) +! iqis(i) += mo_coef_transp(i,r) * iqrs(i,r) +! iqri(i) += mo_coef_transp(i,r) * iqsr(i,r) +! enddo +! enddo +! do i0=1,n_virt_orb +! i= list_virt(i0) +! do j0=1,n_virt_orb +! j = list_virt(j0) +! c = mo_coef_transp(j,q)*mo_coef_transp(j,s) +! mo_two_e_integrals_vv_from_ao(j,i) += c * iqis(i) +! mo_two_e_integrals_vv_exchange_from_ao(j,i) += c * iqri(i) +! enddo +! enddo +! +! enddo +! enddo +! !$OMP END DO NOWAIT +! deallocate(iqrs,iqsr,int_value,int_idx) +! !$OMP END PARALLEL +! +! mo_two_e_integrals_vv_anti_from_ao = mo_two_e_integrals_vv_from_ao - mo_two_e_integrals_vv_exchange_from_ao +! ! print*, '**********' +! ! do i0 =1, n_virt_orb +! ! i = list_virt(i0) +! ! print*, mo_two_e_integrals_vv_from_ao(i,i) +! ! enddo +! ! print*, '**********' +! +! +!END_PROVIDER +! +! +! BEGIN_PROVIDER [ double precision, mo_two_e_integrals_jj, (mo_num,mo_num) ] +!&BEGIN_PROVIDER [ double precision, mo_two_e_integrals_jj_exchange, (mo_num,mo_num) ] +!&BEGIN_PROVIDER [ double precision, mo_two_e_integrals_jj_anti, (mo_num,mo_num) ] +! implicit none +! BEGIN_DOC +! ! mo_two_e_integrals_jj(i,j) = J_ij +! ! mo_two_e_integrals_jj_exchange(i,j) = K_ij +! ! mo_two_e_integrals_jj_anti(i,j) = J_ij - K_ij +! END_DOC +! +! integer :: i,j +! double precision :: get_two_e_integral +! +! PROVIDE mo_two_e_integrals_in_map +! mo_two_e_integrals_jj = 0.d0 +! mo_two_e_integrals_jj_exchange = 0.d0 +! +! do j=1,mo_num +! do i=1,mo_num +! mo_two_e_integrals_jj(i,j) = get_two_e_integral(i,j,i,j,mo_integrals_map) +! mo_two_e_integrals_jj_exchange(i,j) = get_two_e_integral(i,j,j,i,mo_integrals_map) +! mo_two_e_integrals_jj_anti(i,j) = mo_two_e_integrals_jj(i,j) - mo_two_e_integrals_jj_exchange(i,j) +! enddo +! enddo +! +!END_PROVIDER +! +! +!subroutine clear_mo_map +! implicit none +! BEGIN_DOC +! ! Frees the memory of the MO map +! END_DOC +! call map_deinit(mo_integrals_map) +! FREE mo_integrals_map mo_two_e_integrals_jj mo_two_e_integrals_jj_anti +! FREE mo_two_e_integrals_jj_exchange mo_two_e_integrals_in_map +!end +!