From bcc23bf47f930bd58ce0c79ade56a2e2bb3141fc Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Fri, 24 Jan 2020 07:42:37 -0600 Subject: [PATCH] finished complex mapping, starting comples hartree fock --- src/ao_two_e_ints/map_integrals.irp.f | 165 ++++++++---------- src/mo_basis/mos.irp.f | 29 +++ src/mo_one_e_ints/ao_to_mo.irp.f | 68 ++++++++ src/scf_utils/fock_matrix.irp.f | 58 +++++- .../export_integrals_ao_periodic.irp.f | 8 + 5 files changed, 236 insertions(+), 92 deletions(-) diff --git a/src/ao_two_e_ints/map_integrals.irp.f b/src/ao_two_e_ints/map_integrals.irp.f index e9d2740c..4c30c4df 100644 --- a/src/ao_two_e_ints/map_integrals.irp.f +++ b/src/ao_two_e_ints/map_integrals.irp.f @@ -258,6 +258,7 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] complex(integral_kind) :: integral integer(key_kind) :: p,q,r,s,ik,jl logical :: ilek, jlel, iklejl + complex*16 :: get_ao_two_e_integral_periodic_simple !$OMP PARALLEL DO PRIVATE (ilek,jlel,p,q,r,s, ik,jl,iklejl, & @@ -267,36 +268,8 @@ BEGIN_PROVIDER [ complex*16, ao_integrals_cache_periodic, (0:64*64*64*64) ] do j=ao_integrals_cache_min,ao_integrals_cache_max do i=ao_integrals_cache_min,ao_integrals_cache_max !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx1) - ilek = (i.le.k) - jlel = (j.le.l) - idx1 = 2*idx1 - 1 - if (ilek.eqv.jlel) then !map1 - !TODO: merge these calls using map_get_2 - call map_get(ao_integrals_map,idx1,tmp_re) - call map_get(ao_integrals_map,idx1+1,tmp_im) - if (ilek) then - integral = dcmplx(tmp_re,tmp_im) - else - integral = dcmplx(tmp_re,-tmp_im) - endif - else !map2 - !TODO: merge these calls using map_get_2 - call map_get(ao_integrals_map_2,idx1,tmp_re) - call map_get(ao_integrals_map_2,idx1+1,tmp_im) - p = min(i,k) - r = max(i,k) - ik = p+shiftr(r*r-r,1) - q = min(j,l) - s = max(j,l) - jl = q+shiftr(s*s-s,1) - iklejl = (ik.le.jl) - if (ilek.eqv.iklejl) then - integral = dcmplx(tmp_re,tmp_im) - else - integral = dcmplx(tmp_re,-tmp_im) - endif - endif + integral = get_ao_two_e_integral_periodic_simple(i,j,k,l,& + ao_integrals_map,ao_integrals_map_2) ii = l-ao_integrals_cache_min ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) @@ -324,35 +297,61 @@ subroutine ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign) integer(key_kind), intent(out) :: idx logical, intent(out) :: use_map1 double precision, intent(out) :: sign - integer(key_kind) :: p,q,r,s,ik,jl - logical :: ilek, jlel, iklejl, ikeqjl + integer(key_kind) :: p,q,r,s,ik,jl,ij,kl + logical :: iltk, jltl, ikltjl, ieqk, jeql, ikeqjl, ijeqkl ! i.le.k, j.le.l, tri(i,k).le.tri(j,l) !DIR$ FORCEINLINE call two_e_integrals_index_periodic(i,j,k,l,idx,ik,jl) - ilek = (i.le.k) - jlel = (j.le.l) + p = min(i,j) + r = max(i,j) + ij = p+shiftr(r*r-r,1) + q = min(k,l) + s = max(k,l) + kl = q+shiftr(s*s-s,1) + + idx = 2*idx-1 - ikeqjl = (ik.eq.jl) - if (ilek.eqv.jlel) then !map1 + + if (ij==kl) then !real, map1 + sign=0.d0 use_map1=.True. + else + iltk = (i.lt.k) + jltl = (j.lt.l) + ieqk = (i.eq.k) + jeql = (j.eq.l) + ikltjl = (ik.lt.jl) + ikeqjl = (ik.eq.jl) if (ikeqjl) then - sign=0.d0 - else if (ilek) then - sign=1.d0 - else - sign=-1.d0 - endif - else !map2 - use_map1=.False. - if (ikeqjl) then - sign=0.d0 - else - iklejl = (ik.le.jl) - if (ilek.eqv.iklejl) then - sign=1.d0 + if (iltk) then + sign=1.d0 + use_map1=.False. else sign=-1.d0 + use_map1=.False. endif + else if (ieqk) then + if (jltl) then + sign=1.d0 + use_map1=.True. + else + sign=-1.d0 + use_map1=.True. + endif + else if (jeql) then + if (iltk) then + sign=1.d0 + use_map1=.True. + else + sign=-1.d0 + use_map1=.True. + endif + else if (iltk.eqv.ikltjl) then + sign=1.d0 + use_map1=.False. + else + sign=-1.d0 + use_map1=.False. endif endif end @@ -364,48 +363,33 @@ complex*16 function get_ao_two_e_integral_periodic_simple(i,j,k,l,map,map2) resu ! Gets one AO bi-electronic integral from the AO map END_DOC integer, intent(in) :: i,j,k,l - integer(key_kind) :: idx1,idx2 + integer(key_kind) :: idx1,idx2,idx real(integral_kind) :: tmp_re, tmp_im integer(key_kind) :: idx_re,idx_im type(map_type), intent(inout) :: map,map2 integer :: ii complex(integral_kind) :: tmp integer(key_kind) :: p,q,r,s,ik,jl - logical :: ilek, jlel, iklejl + logical :: ilek, jlel, iklejl,use_map1 + double precision :: sign ! a.le.c, b.le.d, tri(a,c).le.tri(b,d) - PROVIDE ao_two_e_integrals_in_map - !DIR$ FORCEINLINE - call two_e_integrals_index(i,j,k,l,idx1) - ilek = (i.le.k) - jlel = (j.le.l) - idx1 = idx1*2-1 - if (ilek.eqv.jlel) then !map1 - !TODO: merge these calls using map_get_2 - call map_get(map,idx1,tmp_re) - call map_get(map,idx1+1,tmp_im) - if (ilek) then - tmp = dcmplx(tmp_re,tmp_im) - else - tmp = dcmplx(tmp_re,-tmp_im) - endif - else !map2 - !TODO: merge these calls using map_get_2 - call map_get(map2,idx1,tmp_re) - call map_get(map2,idx1+1,tmp_im) - p = min(i,k) - r = max(i,k) - ik = p+shiftr(r*r-r,1) - q = min(j,l) - s = max(j,l) - jl = q+shiftr(s*s-s,1) - iklejl = (ik.le.jl) - if (ilek.eqv.iklejl) then - tmp = dcmplx(tmp_re,tmp_im) - else - tmp = dcmplx(tmp_re,-tmp_im) - endif - endif - result = tmp + PROVIDE ao_two_e_integrals_in_map + call ao_two_e_integral_periodic_map_idx_sign(i,j,k,l,use_map1,idx,sign) + if (use_map1) then + call map_get(map,idx,tmp_re) + if (sign/=0.d0) then + call map_get(map,idx+1,tmp_im) + tmp_im *= sign + else + tmp_im=0.d0 + endif + else + call map_get(map2,idx,tmp_re) + call map_get(map2,idx+1,tmp_im) + tmp_im *= sign + endif + tmp = dcmplx(tmp_re,tmp_im) + result = tmp end @@ -428,11 +412,12 @@ complex*16 function get_ao_two_e_integral_periodic(i,j,k,l,map,map2) result(resu ! a.le.c, b.le.d, tri(a,c).le.tri(b,d) PROVIDE ao_two_e_integrals_in_map ao_integrals_cache_periodic ao_integrals_cache_min !DIR$ FORCEINLINE - if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then - tmp = (0.d0,0.d0) - else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then - tmp = (0.d0,0.d0) - else +! if (ao_overlap_abs(i,k)*ao_overlap_abs(j,l) < ao_integrals_threshold ) then +! tmp = (0.d0,0.d0) +! else if (ao_two_e_integral_schwartz(i,k)*ao_two_e_integral_schwartz(j,l) < ao_integrals_threshold) then +! tmp = (0.d0,0.d0) +! else + if (.True.) then ii = l-ao_integrals_cache_min ii = ior(ii, k-ao_integrals_cache_min) ii = ior(ii, j-ao_integrals_cache_min) diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 73d33901..23713f34 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -280,6 +280,35 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) deallocate(T) end +subroutine ao_to_mo_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 + ! + ! Ct.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,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('C','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 mix_mo_jk(j,k) implicit none diff --git a/src/mo_one_e_ints/ao_to_mo.irp.f b/src/mo_one_e_ints/ao_to_mo.irp.f index a0d8caaa..5279608f 100644 --- a/src/mo_one_e_ints/ao_to_mo.irp.f +++ b/src/mo_one_e_ints/ao_to_mo.irp.f @@ -50,6 +50,61 @@ subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao) deallocate(T) end +subroutine mo_to_ao_complex(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the AO basis + ! + ! (S.C).A_mo.(S.C)t + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_mo(LDA_mo,mo_num) + complex*16, intent(out) :: A_ao(LDA_ao,ao_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(mo_num,ao_num) ) + + call zgemm('N','C', mo_num, ao_num, mo_num, & + (1.d0,0.d0), A_mo,size(A_mo,1), & + S_mo_coef_complex, size(S_mo_coef_complex,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('N','N', ao_num, ao_num, mo_num, & + (1.d0,0.d0), S_mo_coef_complex, size(S_mo_coef_complex,1), & + T, size(T,1), & + (0.d0,0.d0), A_ao, size(A_ao,1)) + + deallocate(T) +end + +subroutine mo_to_ao_no_overlap_complex(A_mo,LDA_mo,A_ao,LDA_ao) + implicit none + BEGIN_DOC + ! Transform A from the MO basis to the S^-1 AO basis + ! Useful for density matrix + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_mo(LDA_mo,mo_num) + complex*16, intent(out) :: A_ao(LDA_ao,ao_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(mo_num,ao_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call zgemm('N','C', mo_num, ao_num, mo_num, & + (1.d0,0.d0), A_mo,size(A_mo,1), & + mo_coef_complex, size(mo_coef_complex,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('N','N', ao_num, ao_num, mo_num, & + (1.d0,0.d0), mo_coef_complex,size(mo_coef_complex,1), & + T, size(T,1), & + (0.d0,0.d0), A_ao, size(A_ao,1)) + + deallocate(T) +end + + BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num, mo_num) ] implicit none BEGIN_DOC @@ -63,4 +118,17 @@ BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num, mo_num) ] END_PROVIDER +BEGIN_PROVIDER [ complex*16, S_mo_coef_complex, (ao_num, mo_num) ] + implicit none + BEGIN_DOC + ! Product S.C where S is the overlap matrix in the AO basis and C the mo_coef matrix. + END_DOC + + call zgemm('N','N',ao_num, mo_num, ao_num, (1.d0,0.d0), & + ao_overlap_complex, size(ao_overlap_complex,1), & + mo_coef_complex, size(mo_coef_complex,1), & + (0.d0,0.d0), & + S_mo_coef_complex, size(S_mo_coef_complex,1)) + +END_PROVIDER diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index fc9eaadd..a49d51e5 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -101,18 +101,45 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_num,mo_num) ] BEGIN_DOC ! Fock matrix on the MO basis END_DOC - call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & + if (is_periodic) then + print*,'error',irp_here + stop -1 + else + call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), & Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1)) + endif END_PROVIDER +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo_complex(Fock_matrix_ao_alpha_complex,size(Fock_matrix_ao_alpha_complex,1), & + Fock_matrix_mo_alpha_complex,size(Fock_matrix_mo_alpha_complex,1)) +END_PROVIDER BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_num,mo_num) ] implicit none BEGIN_DOC ! Fock matrix on the MO basis END_DOC - call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & + if (is_periodic) then + print*,'error',irp_here + stop -1 + else + call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), & Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1)) + endif +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo_complex(Fock_matrix_ao_beta_complex,size(Fock_matrix_ao_beta_complex,1), & + Fock_matrix_mo_beta_complex,size(Fock_matrix_mo_beta_complex,1)) END_PROVIDER @@ -143,6 +170,33 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] END_PROVIDER +BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if(frozen_orb_scf)then + call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & + Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) + else + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_complex(i,j) = Fock_matrix_ao_alpha_complex(i,j) + enddo + enddo + else + call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & + Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) + endif + endif +END_PROVIDER + + BEGIN_PROVIDER [ double precision, SCF_energy ] implicit none BEGIN_DOC diff --git a/src/utils_periodic/export_integrals_ao_periodic.irp.f b/src/utils_periodic/export_integrals_ao_periodic.irp.f index 8d60ff49..1f190749 100644 --- a/src/utils_periodic/export_integrals_ao_periodic.irp.f +++ b/src/utils_periodic/export_integrals_ao_periodic.irp.f @@ -153,6 +153,14 @@ provide ao_two_e_integrals_in_map call map_get(ao_integrals_map_2,idx_tmp+1,tmp6) print*,tmp3,tmp4 print*,tmp5,tmp6 + integer*8 :: ii + ii = l-ao_integrals_cache_min + ii = ior( shiftl(ii,6), k-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), j-ao_integrals_cache_min) + ii = ior( shiftl(ii,6), i-ao_integrals_cache_min) + print*,'cache(pbc)=', ao_integrals_cache_periodic(ii) + print*,'cache(old)=', ao_integrals_cache(ii) + print* ! if (use_map1) then ! n_integrals_1 += 1 ! buffer_i_1(n_integrals_1-1)=idx_tmp