From df4c9431d086c27f219d34ad4a542bfb8b0d6c01 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Thu, 5 Mar 2020 15:53:45 +0100 Subject: [PATCH] Fixed compilation problems --- src/generators_cas/generators.irp.f | 36 ++++ src/hartree_fock/fock_matrix_hf_cplx.irp.f | 213 --------------------- src/scf_utils/fock_matrix_cplx.irp.f | 213 +++++++++++++++++++++ src/single_ref_method/generators.irp.f | 22 +++ 4 files changed, 271 insertions(+), 213 deletions(-) delete mode 100644 src/hartree_fock/fock_matrix_hf_cplx.irp.f diff --git a/src/generators_cas/generators.irp.f b/src/generators_cas/generators.irp.f index b2f58202..cc87a0af 100644 --- a/src/generators_cas/generators.irp.f +++ b/src/generators_cas/generators.irp.f @@ -82,3 +82,39 @@ BEGIN_PROVIDER [ double precision, select_max, (size_select_max) ] select_max = huge(1.d0) END_PROVIDER + + BEGIN_PROVIDER [ complex*16, psi_coef_generators_complex, (psi_det_size,N_states) ] +&BEGIN_PROVIDER [ complex*16, psi_coef_sorted_gen_complex, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! For Single reference wave functions, the generator is the + ! Hartree-Fock determinant + END_DOC + integer :: i, k, l, m + logical :: good + integer, external :: number_of_holes,number_of_particles + integer, allocatable :: nongen(:) + integer :: inongen + + allocate(nongen(N_det)) + + inongen = 0 + m=0 + do i=1,N_det + good = ( number_of_holes(psi_det_sorted(1,1,i)) ==0).and.(number_of_particles(psi_det_sorted(1,1,i))==0 ) + if (good) then + m = m+1 + psi_coef_generators_complex(m,:) = psi_coef_sorted_complex(i,:) + else + inongen += 1 + nongen(inongen) = i + endif + enddo + ASSERT (m == N_det_generators) + + psi_coef_sorted_gen_complex(:N_det_generators, :) = psi_coef_generators_complex(:N_det_generators, :) + do i=1,inongen + psi_coef_sorted_gen_complex(N_det_generators+i, :) = psi_coef_sorted_complex(nongen(i),:) + end do +END_PROVIDER + diff --git a/src/hartree_fock/fock_matrix_hf_cplx.irp.f b/src/hartree_fock/fock_matrix_hf_cplx.irp.f deleted file mode 100644 index 3dea06fe..00000000 --- a/src/hartree_fock/fock_matrix_hf_cplx.irp.f +++ /dev/null @@ -1,213 +0,0 @@ - - BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_complex, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_complex , (ao_num, ao_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_complex = (0.d0,0.d0) - ao_two_e_integral_beta_complex = (0.d0,0.d0) - PROVIDE ao_two_e_integrals_in_map - - 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 - - !$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 c0,key1)& - !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha_complex, & - !$OMP SCF_density_matrix_ao_beta_complex, & - !$OMP ao_integrals_map, ao_two_e_integral_alpha_complex, ao_two_e_integral_beta_complex) - - 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,ao_num), & - ao_two_e_integral_beta_tmp(ao_num,ao_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 - 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) - 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)*() - - c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral - - ao_two_e_integral_alpha_tmp(i,k) += c0 - ao_two_e_integral_beta_tmp (i,k) += c0 - - ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral - ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral - 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) - integral = values(k1) - - c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral - - ao_two_e_integral_alpha_tmp(i,k) += c0 - ao_two_e_integral_beta_tmp (i,k) += c0 - - ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral - ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral - enddo - endif - enddo - enddo - !$OMP END DO NOWAIT - !$OMP CRITICAL - ao_two_e_integral_alpha_complex += ao_two_e_integral_alpha_tmp - ao_two_e_integral_beta_complex += 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 c0,key1)& - !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha_complex, & - !$OMP SCF_density_matrix_ao_beta_complex, & - !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_complex, ao_two_e_integral_beta_complex) - - 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,ao_num), & - ao_two_e_integral_beta_tmp(ao_num,ao_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 - 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) - 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)*() - - c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral - - ao_two_e_integral_alpha_tmp(i,k) += c0 - ao_two_e_integral_beta_tmp (i,k) += c0 - - ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral - ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral - 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) - integral = values(k1) - - c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral - - ao_two_e_integral_alpha_tmp(i,k) += c0 - ao_two_e_integral_beta_tmp (i,k) += c0 - - ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral - ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral - enddo - endif - enddo - enddo - !$OMP END DO NOWAIT - !$OMP CRITICAL - ao_two_e_integral_alpha_complex += ao_two_e_integral_alpha_tmp - ao_two_e_integral_beta_complex += 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, Fock_matrix_ao_alpha_complex, (ao_num, ao_num) ] -&BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_beta_complex, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Alpha Fock matrix in AO basis set - END_DOC - - integer :: i,j - do j=1,ao_num - do i=1,ao_num - Fock_matrix_ao_alpha_complex(i,j) = ao_one_e_integrals_complex(i,j) + ao_two_e_integral_alpha_complex(i,j) - Fock_matrix_ao_beta_complex (i,j) = ao_one_e_integrals_complex(i,j) + ao_two_e_integral_beta_complex (i,j) - enddo - enddo - -END_PROVIDER diff --git a/src/scf_utils/fock_matrix_cplx.irp.f b/src/scf_utils/fock_matrix_cplx.irp.f index 290f9b9d..577fe5c2 100644 --- a/src/scf_utils/fock_matrix_cplx.irp.f +++ b/src/scf_utils/fock_matrix_cplx.irp.f @@ -146,3 +146,216 @@ BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ] endif END_PROVIDER + + BEGIN_PROVIDER [ complex*16, ao_two_e_integral_alpha_complex, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ complex*16, ao_two_e_integral_beta_complex , (ao_num, ao_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_complex = (0.d0,0.d0) + ao_two_e_integral_beta_complex = (0.d0,0.d0) + PROVIDE ao_two_e_integrals_in_map + + 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 + + !$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 c0,key1)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha_complex, & + !$OMP SCF_density_matrix_ao_beta_complex, & + !$OMP ao_integrals_map, ao_two_e_integral_alpha_complex, ao_two_e_integral_beta_complex) + + 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,ao_num), & + ao_two_e_integral_beta_tmp(ao_num,ao_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 + 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) + 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)*() + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + 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) + integral = values(k1) + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_complex += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_complex += 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 c0,key1)& + !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha_complex, & + !$OMP SCF_density_matrix_ao_beta_complex, & + !$OMP ao_integrals_map_2, ao_two_e_integral_alpha_complex, ao_two_e_integral_beta_complex) + + 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,ao_num), & + ao_two_e_integral_beta_tmp(ao_num,ao_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 + 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) + 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)*() + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + 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) + integral = values(k1) + + c0 = (scf_density_matrix_ao_alpha_complex(l,j)+scf_density_matrix_ao_beta_complex(l,j)) * integral + + ao_two_e_integral_alpha_tmp(i,k) += c0 + ao_two_e_integral_beta_tmp (i,k) += c0 + + ao_two_e_integral_alpha_tmp(i,l) -= SCF_density_matrix_ao_alpha_complex(k,j) * integral + ao_two_e_integral_beta_tmp (i,l) -= scf_density_matrix_ao_beta_complex (k,j) * integral + enddo + endif + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + ao_two_e_integral_alpha_complex += ao_two_e_integral_alpha_tmp + ao_two_e_integral_beta_complex += 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, Fock_matrix_ao_alpha_complex, (ao_num, ao_num) ] +&BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_beta_complex, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Alpha Fock matrix in AO basis set + END_DOC + + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_alpha_complex(i,j) = ao_one_e_integrals_complex(i,j) + ao_two_e_integral_alpha_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) = ao_one_e_integrals_complex(i,j) + ao_two_e_integral_beta_complex (i,j) + enddo + enddo + +END_PROVIDER diff --git a/src/single_ref_method/generators.irp.f b/src/single_ref_method/generators.irp.f index ce71f996..860f357a 100644 --- a/src/single_ref_method/generators.irp.f +++ b/src/single_ref_method/generators.irp.f @@ -25,6 +25,7 @@ END_PROVIDER psi_det_generators(i,2,1) = HF_bitmask(i,2) enddo + ! Search for HF determinant do j=1,N_det call get_excitation_degree(HF_bitmask,psi_det(1,1,j),degree,N_int) if (degree == 0) then @@ -55,4 +56,25 @@ BEGIN_PROVIDER [ integer, size_select_max ] END_PROVIDER +BEGIN_PROVIDER [ complex*16, psi_coef_generators_complex, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! Complex variant of psi_coef_generators + END_DOC + + integer :: i,j,k + integer :: degree + + ! Search for HF determinant + do j=1,N_det + call get_excitation_degree(HF_bitmask,psi_det(1,1,j),degree,N_int) + if (degree == 0) then + k = j + exit + endif + end do + + psi_coef_generators_complex(1,:) = psi_coef_generators_complex(j,:) + +END_PROVIDER