diff --git a/src/determinants/density_matrix.irp.f b/src/determinants/density_matrix.irp.f index e69a1803..4cc20cc1 100644 --- a/src/determinants/density_matrix.irp.f +++ b/src/determinants/density_matrix.irp.f @@ -248,6 +248,11 @@ BEGIN_PROVIDER [ double precision, one_e_spin_density_mo, (mo_num,mo_num) ] END_PROVIDER subroutine set_natural_mos + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif implicit none BEGIN_DOC ! Set natural orbitals, obtained by diagonalization of the one-body density matrix @@ -274,6 +279,11 @@ subroutine set_natural_mos end subroutine save_natural_mos + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif implicit none BEGIN_DOC ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in @@ -292,11 +302,19 @@ BEGIN_PROVIDER [ double precision, c0_weight, (N_states) ] if (N_states > 1) then integer :: i double precision :: c + if (is_complex) then + do i=1,N_states + c0_weight(i) = 1.d-31 + c = maxval(cdabs(psi_coef_complex(:,i) * psi_coef_complex(:,i))) + c0_weight(i) = 1.d0/(c+1.d-20) + enddo + else do i=1,N_states c0_weight(i) = 1.d-31 c = maxval(psi_coef(:,i) * psi_coef(:,i)) c0_weight(i) = 1.d0/(c+1.d-20) enddo + endif c = 1.d0/minval(c0_weight(:)) do i=1,N_states c0_weight(i) = c0_weight(i) * c @@ -398,8 +416,23 @@ subroutine get_occupation_from_dets(istate,occupation) ASSERT (istate <= N_states) occupation = 0.d0 - double precision, external :: u_dot_u + + if (is_complex) then + double precision, external :: u_dot_u_complex + norm_2 = 1.d0/u_dot_u_complex(psi_coef_complex(1,istate),N_det) + do i=1,N_det + c = cdabs(psi_coef_complex(i,istate)*psi_coef_complex(i,istate))*norm_2 + call bitstring_to_list_ab(psi_det(1,1,i), list, n_elements, N_int) + do ispin=1,2 + do j=1,n_elements(ispin) + ASSERT ( list(j,ispin) < mo_num ) + occupation( list(j,ispin) ) += c + enddo + enddo + enddo + else + double precision, external :: u_dot_u norm_2 = 1.d0/u_dot_u(psi_coef(1,istate),N_det) do i=1,N_det @@ -412,5 +445,6 @@ subroutine get_occupation_from_dets(istate,occupation) enddo enddo enddo + endif end diff --git a/src/determinants/density_matrix_complex.irp.f b/src/determinants/density_matrix_complex.irp.f new file mode 100644 index 00000000..e6c81d33 --- /dev/null +++ b/src/determinants/density_matrix_complex.irp.f @@ -0,0 +1,311 @@ + BEGIN_PROVIDER [ complex*16, one_e_dm_mo_alpha_average_complex, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ complex*16, one_e_dm_mo_beta_average_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! $\alpha$ and $\beta$ one-body density matrix for each state + END_DOC + integer :: i + one_e_dm_mo_alpha_average_complex = (0.d0,0.d0) + one_e_dm_mo_beta_average_complex = (0.d0,0.d0) + do i = 1,N_states + one_e_dm_mo_alpha_average_complex(:,:) += one_e_dm_mo_alpha_complex(:,:,i) * state_average_weight(i) + one_e_dm_mo_beta_average_complex(:,:) += one_e_dm_mo_beta_complex(:,:,i) * state_average_weight(i) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, one_e_dm_mo_diff_complex, (mo_num,mo_num,2:N_states) ] + implicit none + BEGIN_DOC + ! Difference of the one-body density matrix with respect to the ground state + END_DOC + integer :: i,j, istate + + do istate=2,N_states + do j=1,mo_num + do i=1,mo_num + one_e_dm_mo_diff_complex(i,j,istate) = & + one_e_dm_mo_alpha_complex(i,j,istate) - one_e_dm_mo_alpha_complex(i,j,1) +& + one_e_dm_mo_beta_complex (i,j,istate) - one_e_dm_mo_beta_complex (i,j,1) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, one_e_dm_mo_spin_index_complex, (mo_num,mo_num,N_states,2) ] + implicit none + integer :: i,j,ispin,istate + ispin = 1 + do istate = 1, N_states + do j = 1, mo_num + do i = 1, mo_num + one_e_dm_mo_spin_index_complex(i,j,istate,ispin) = one_e_dm_mo_alpha_complex(i,j,istate) + enddo + enddo + enddo + + ispin = 2 + do istate = 1, N_states + do j = 1, mo_num + do i = 1, mo_num + one_e_dm_mo_spin_index_complex(i,j,istate,ispin) = one_e_dm_mo_beta_complex(i,j,istate) + enddo + enddo + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, one_e_dm_dagger_mo_spin_index_complex, (mo_num,mo_num,N_states,2) ] + print*,irp_here,' not implemented for complex' + stop -1 +! implicit none +! integer :: i,j,ispin,istate +! ispin = 1 +! do istate = 1, N_states +! do j = 1, mo_num +! one_e_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_e_dm_mo_alpha(j,j,istate) +! do i = j+1, mo_num +! one_e_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_e_dm_mo_alpha(i,j,istate) +! one_e_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_e_dm_mo_alpha(i,j,istate) +! enddo +! enddo +! enddo +! +! ispin = 2 +! do istate = 1, N_states +! do j = 1, mo_num +! one_e_dm_dagger_mo_spin_index(j,j,istate,ispin) = 1 - one_e_dm_mo_beta(j,j,istate) +! do i = j+1, mo_num +! one_e_dm_dagger_mo_spin_index(i,j,istate,ispin) = -one_e_dm_mo_beta(i,j,istate) +! one_e_dm_dagger_mo_spin_index(j,i,istate,ispin) = -one_e_dm_mo_beta(i,j,istate) +! enddo +! enddo +! enddo +! +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, one_e_dm_mo_alpha_complex, (mo_num,mo_num,N_states) ] +&BEGIN_PROVIDER [ complex*16, one_e_dm_mo_beta_complex, (mo_num,mo_num,N_states) ] + implicit none + BEGIN_DOC + ! $\alpha$ and $\beta$ one-body density matrix for each state + ! $\gamma_{\mu\nu} = \langle\Psi|a_{\nu}^{\dagger}a_{\mu}|\Psi\rangle$ + ! $\gamma_{\mu\nu} = \langle a_{\nu} \Psi|a_{\mu} \Psi\rangle$ + ! $\gamma_{\mu\nu} = \sum_{IJ} c^*_J c_I \langle a_{\nu} I|a_{\mu} J\rangle$ + END_DOC + + integer :: j,k,l,m,k_a,k_b + integer :: occ(N_int*bit_kind_size,2) + complex*16 :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int) + integer :: exc(0:2,2),n_occ(2) + complex*16, allocatable :: tmp_a(:,:,:), tmp_b(:,:,:) + integer :: krow, kcol, lrow, lcol + + PROVIDE psi_det + + one_e_dm_mo_alpha_complex = (0.d0,0.d0) + one_e_dm_mo_beta_complex = (0.d0,0.d0) + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,& + !$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)& + !$OMP SHARED(psi_det,psi_coef_complex,N_int,N_states,elec_alpha_num, & + !$OMP elec_beta_num,one_e_dm_mo_alpha_complex,one_e_dm_mo_beta_complex,N_det,& + !$OMP mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,& + !$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,& + !$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,& + !$OMP psi_bilinear_matrix_values_complex, psi_bilinear_matrix_transp_values_complex,& + !$OMP N_det_alpha_unique,N_det_beta_unique,irp_here) + allocate(tmp_a(mo_num,mo_num,N_states), tmp_b(mo_num,mo_num,N_states) ) + tmp_a = (0.d0,0.d0) + !$OMP DO SCHEDULE(dynamic,64) + do k_a=1,N_det + krow = psi_bilinear_matrix_rows(k_a) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_columns(k_a) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) + + ! Diagonal part + ! ------------- + + call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) + do m=1,N_states + ck = cdabs(psi_bilinear_matrix_values_complex(k_a,m)*psi_bilinear_matrix_values_complex(k_a,m)) + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j,m) += ck + enddo + enddo + + if (k_a == N_det) cycle + l = k_a+1 + lrow = psi_bilinear_matrix_rows(l) + lcol = psi_bilinear_matrix_columns(l) + ! Fix beta determinant, loop over alphas + do while ( lcol == kcol ) + tmp_det2(:) = psi_det_alpha_unique(:, lrow) + call get_excitation_degree_spin(tmp_det(1,1),tmp_det2,degree,N_int) + if (degree == 1) then + exc = 0 + call get_single_excitation_spin(tmp_det(1,1),tmp_det2,exc,phase,N_int) + call decode_exc_spin(exc,h1,p1,h2,p2) + ! h1 occ in k + ! p1 occ in l + do m=1,N_states + ckl = dconjg(psi_bilinear_matrix_values_complex(k_a,m))*psi_bilinear_matrix_values_complex(l,m) * phase + tmp_a(h1,p1,m) += dconjg(ckl) + tmp_a(p1,h1,m) += ckl + enddo + endif + l = l+1 + if (l>N_det) exit + lrow = psi_bilinear_matrix_rows(l) + lcol = psi_bilinear_matrix_columns(l) + enddo + + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + one_e_dm_mo_alpha_complex(:,:,:) = one_e_dm_mo_alpha_complex(:,:,:) + tmp_a(:,:,:) + !$OMP END CRITICAL + deallocate(tmp_a) + + tmp_b = (0.d0,0.d0) + !$OMP DO SCHEDULE(dynamic,64) + do k_b=1,N_det + krow = psi_bilinear_matrix_transp_rows(k_b) + ASSERT (krow <= N_det_alpha_unique) + + kcol = psi_bilinear_matrix_transp_columns(k_b) + ASSERT (kcol <= N_det_beta_unique) + + tmp_det(1:N_int,1) = psi_det_alpha_unique(1:N_int,krow) + tmp_det(1:N_int,2) = psi_det_beta_unique (1:N_int,kcol) + + ! Diagonal part + ! ------------- + + call bitstring_to_list_ab(tmp_det, occ, n_occ, N_int) + do m=1,N_states + ck = cdabs(psi_bilinear_matrix_transp_values_complex(k_b,m)*psi_bilinear_matrix_transp_values_complex(k_b,m)) + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j,m) += ck + enddo + enddo + + if (k_b == N_det) cycle + l = k_b+1 + lrow = psi_bilinear_matrix_transp_rows(l) + lcol = psi_bilinear_matrix_transp_columns(l) + ! Fix beta determinant, loop over alphas + do while ( lrow == krow ) + tmp_det2(:) = psi_det_beta_unique(:, lcol) + call get_excitation_degree_spin(tmp_det(1,2),tmp_det2,degree,N_int) + if (degree == 1) then + exc = 0 + call get_single_excitation_spin(tmp_det(1,2),tmp_det2,exc,phase,N_int) + call decode_exc_spin(exc,h1,p1,h2,p2) + do m=1,N_states + ckl = dconjg(psi_bilinear_matrix_transp_values_complex(k_b,m))*psi_bilinear_matrix_transp_values_complex(l,m) * phase + tmp_b(h1,p1,m) += dconjg(ckl) + tmp_b(p1,h1,m) += ckl + enddo + endif + l = l+1 + if (l>N_det) exit + lrow = psi_bilinear_matrix_transp_rows(l) + lcol = psi_bilinear_matrix_transp_columns(l) + enddo + + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_e_dm_mo_beta_complex(:,:,:) = one_e_dm_mo_beta_complex(:,:,:) + tmp_b(:,:,:) + !$OMP END CRITICAL + + deallocate(tmp_b) + !$OMP END PARALLEL + +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, one_e_dm_mo_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! One-body density matrix + END_DOC + one_e_dm_mo_complex = one_e_dm_mo_alpha_average_complex + one_e_dm_mo_beta_average_complex +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, one_e_spin_density_mo_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! $\rho(\alpha) - \rho(\beta)$ + END_DOC + one_e_spin_density_mo_complex = one_e_dm_mo_alpha_average_complex - one_e_dm_mo_beta_average_complex +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, one_e_spin_density_ao_complex, (ao_num,ao_num) ] + BEGIN_DOC + ! One body spin density matrix on the |AO| basis : $\rho_{AO}(\alpha) - \rho_{AO}(\beta)$ + ! todo: verify that this is correct for complex + ! equivalent to using mo_to_ao_no_overlap? + END_DOC + implicit none + integer :: i,j,k,l + complex*16 :: dm_mo + + one_e_spin_density_ao_complex = (0.d0,0.d0) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, mo_num + do j = 1, mo_num + dm_mo = one_e_spin_density_mo_complex(j,i) + ! if(dabs(dm_mo).le.1.d-10)cycle + one_e_spin_density_ao_complex(l,k) += dconjg(mo_coef_complex(k,i)) * mo_coef_complex(l,j) * dm_mo + + enddo + enddo + enddo + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, one_e_dm_ao_alpha_complex, (ao_num,ao_num) ] +&BEGIN_PROVIDER [ complex*16, one_e_dm_ao_beta_complex, (ao_num,ao_num) ] + BEGIN_DOC + ! One body density matrix on the |AO| basis : $\rho_{AO}(\alpha), \rho_{AO}(\beta)$. + END_DOC + implicit none + integer :: i,j,k,l + complex*16 :: mo_alpha,mo_beta + + one_e_dm_ao_alpha = (0.d0,0.d0) + one_e_dm_ao_beta = (0.d0,0.d0) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, mo_num + do j = 1, mo_num + mo_alpha = one_e_dm_mo_alpha_average_complex(j,i) + mo_beta = one_e_dm_mo_beta_average_complex(j,i) + ! if(dabs(dm_mo).le.1.d-10)cycle + one_e_dm_ao_alpha_complex(l,k) += dconjg(mo_coef_complex(k,i)) * mo_coef_complex(l,j) * mo_alpha + one_e_dm_ao_beta_complex(l,k) += dconjg(mo_coef_complex(k,i)) * mo_coef_complex(l,j) * mo_beta + enddo + enddo + enddo + enddo + +END_PROVIDER + + diff --git a/src/determinants/s2_complex.irp.f b/src/determinants/s2_complex.irp.f index bb368f55..e2116db8 100644 --- a/src/determinants/s2_complex.irp.f +++ b/src/determinants/s2_complex.irp.f @@ -1,4 +1,5 @@ subroutine u_0_S2_u_0_complex(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8) + !todo: modify/implement for complex print*,irp_here,' not implemented for complex' stop -1 ! use bitmasks @@ -28,6 +29,7 @@ end subroutine S2_u_0_complex(v_0,u_0,n,keys_tmp,Nint) + !todo: modify/implement for complex print*,irp_here,' not implemented for complex' stop -1 ! use bitmasks @@ -46,6 +48,7 @@ subroutine S2_u_0_complex(v_0,u_0,n,keys_tmp,Nint) end subroutine S2_u_0_nstates_complex(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8) + !todo: modify/implement for complex print*,irp_here,' not implemented for complex' stop -1 ! use bitmasks @@ -180,6 +183,7 @@ end subroutine get_uJ_s2_uI_complex(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates) + !todo: modify/implement for complex print*,irp_here,' not implemented for complex' stop -1 ! implicit none @@ -232,6 +236,7 @@ end subroutine i_S2_psi_minilist_complex(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_S2_psi_array) + !todo: modify/implement for complex print*,irp_here,' not implemented for complex' stop -1 ! use bitmasks diff --git a/src/determinants/single_excitations.irp.f b/src/determinants/single_excitations.irp.f index ccfeaa2e..65c8ac7f 100644 --- a/src/determinants/single_excitations.irp.f +++ b/src/determinants/single_excitations.irp.f @@ -159,3 +159,148 @@ subroutine get_single_excitation_from_fock(det_1,det_2,h,p,spin,phase,hij) end + + +BEGIN_PROVIDER [complex*16, fock_operator_closed_shell_ref_bitmask_complex, (mo_num, mo_num) ] + implicit none + integer :: i0,j0,i,j,k0,k + integer :: n_occ_ab(2) + integer :: occ(N_int*bit_kind_size,2) + integer :: n_occ_ab_virt(2) + integer :: occ_virt(N_int*bit_kind_size,2) + integer(bit_kind) :: key_test(N_int) + integer(bit_kind) :: key_virt(N_int,2) + complex*16 :: accu + + call bitstring_to_list_ab(ref_closed_shell_bitmask, occ, n_occ_ab, N_int) + do i = 1, N_int + key_virt(i,1) = full_ijkl_bitmask(i) + key_virt(i,2) = full_ijkl_bitmask(i) + key_virt(i,1) = xor(key_virt(i,1),ref_closed_shell_bitmask(i,1)) + key_virt(i,2) = xor(key_virt(i,2),ref_closed_shell_bitmask(i,2)) + enddo + complex*16, allocatable :: array_coulomb(:),array_exchange(:) + allocate (array_coulomb(mo_num),array_exchange(mo_num)) + call bitstring_to_list_ab(key_virt, occ_virt, n_occ_ab_virt, N_int) + ! docc ---> virt single excitations + do i0 = 1, n_occ_ab(1) + i=occ(i0,1) + do j0 = 1, n_occ_ab_virt(1) + j = occ_virt(j0,1) + ! + call get_mo_two_e_integrals_coulomb_ii_complex(i,j,mo_num,array_coulomb,mo_integrals_map,mo_integrals_map_2) + ! + call get_mo_two_e_integrals_exch_ii_complex(i,j,mo_num,array_exchange,mo_integrals_map,mo_integrals_map_2) + accu = (0.d0,0.d0) + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * array_coulomb(k) - array_exchange(k) + enddo + fock_operator_closed_shell_ref_bitmask_complex(i,j) = accu + mo_one_e_integrals_complex(i,j) + !fock_operator_closed_shell_ref_bitmask_complex(j,i) = dconjg(accu) + mo_one_e_integrals_complex(j,i) + fock_operator_closed_shell_ref_bitmask_complex(j,i) = dconjg(fock_operator_closed_shell_ref_bitmask_complex(i,j)) + enddo + enddo + + ! virt ---> virt single excitations + do i0 = 1, n_occ_ab_virt(1) + i=occ_virt(i0,1) + do j0 = 1, n_occ_ab_virt(1) + j = occ_virt(j0,1) + call get_mo_two_e_integrals_coulomb_ii_complex(i,j,mo_num,array_coulomb,mo_integrals_map,mo_integrals_map_2) + call get_mo_two_e_integrals_exch_ii_complex(i,j,mo_num,array_exchange,mo_integrals_map,mo_integrals_map_2) + accu = (0.d0,0.d0) + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * array_coulomb(k) - array_exchange(k) + enddo + fock_operator_closed_shell_ref_bitmask_complex(i,j) = accu+ mo_one_e_integrals_complex(i,j) + fock_operator_closed_shell_ref_bitmask_complex(j,i) = dconjg(accu)+ mo_one_e_integrals_complex(j,i) + enddo + enddo + + ! docc ---> docc single excitations + do i0 = 1, n_occ_ab(1) + i=occ(i0,1) + do j0 = 1, n_occ_ab(1) + j = occ(j0,1) + call get_mo_two_e_integrals_coulomb_ii_complex(i,j,mo_num,array_coulomb,mo_integrals_map,mo_integrals_map_2) + call get_mo_two_e_integrals_exch_ii_complex(i,j,mo_num,array_exchange,mo_integrals_map,mo_integrals_map_2) + accu = (0.d0,0.d0) + do k0 = 1, n_occ_ab(1) + k = occ(k0,1) + accu += 2.d0 * array_coulomb(k) - array_exchange(k) + enddo + fock_operator_closed_shell_ref_bitmask_complex(i,j) = accu+ mo_one_e_integrals_complex(i,j) + fock_operator_closed_shell_ref_bitmask_complex(j,i) = dconjg(accu)+ mo_one_e_integrals_complex(j,i) + enddo + enddo + deallocate(array_coulomb,array_exchange) + +END_PROVIDER + +subroutine get_single_excitation_from_fock_complex(det_1,det_2,h,p,spin,phase,hij) + use bitmasks + implicit none + integer,intent(in) :: h,p,spin + double precision, intent(in) :: phase + integer(bit_kind), intent(in) :: det_1(N_int,2), det_2(N_int,2) + complex*16, intent(out) :: hij + integer(bit_kind) :: differences(N_int,2) + integer(bit_kind) :: hole(N_int,2) + integer(bit_kind) :: partcl(N_int,2) + integer :: occ_hole(N_int*bit_kind_size,2) + integer :: occ_partcl(N_int*bit_kind_size,2) + integer :: n_occ_ab_hole(2),n_occ_ab_partcl(2) + integer :: i0,i + complex*16 :: buffer_c(mo_num),buffer_x(mo_num) + do i=1, mo_num + buffer_c(i) = big_array_coulomb_integrals_complex(i,h,p) + buffer_x(i) = big_array_exchange_integrals_complex(i,h,p) + enddo + do i = 1, N_int + differences(i,1) = xor(det_1(i,1),ref_closed_shell_bitmask(i,1)) + differences(i,2) = xor(det_1(i,2),ref_closed_shell_bitmask(i,2)) + hole(i,1) = iand(differences(i,1),ref_closed_shell_bitmask(i,1)) + hole(i,2) = iand(differences(i,2),ref_closed_shell_bitmask(i,2)) + partcl(i,1) = iand(differences(i,1),det_1(i,1)) + partcl(i,2) = iand(differences(i,2),det_1(i,2)) + enddo + call bitstring_to_list_ab(hole, occ_hole, n_occ_ab_hole, N_int) + call bitstring_to_list_ab(partcl, occ_partcl, n_occ_ab_partcl, N_int) + hij = fock_operator_closed_shell_ref_bitmask_complex(h,p) + ! holes :: direct terms + do i0 = 1, n_occ_ab_hole(1) + i = occ_hole(i0,1) + hij -= buffer_c(i) + enddo + do i0 = 1, n_occ_ab_hole(2) + i = occ_hole(i0,2) + hij -= buffer_c(i) + enddo + + ! holes :: exchange terms + do i0 = 1, n_occ_ab_hole(spin) + i = occ_hole(i0,spin) + hij += buffer_x(i) + enddo + + ! particles :: direct terms + do i0 = 1, n_occ_ab_partcl(1) + i = occ_partcl(i0,1) + hij += buffer_c(i) + enddo + do i0 = 1, n_occ_ab_partcl(2) + i = occ_partcl(i0,2) + hij += buffer_c(i) + enddo + + ! particles :: exchange terms + do i0 = 1, n_occ_ab_partcl(spin) + i = occ_partcl(i0,spin) + hij -= buffer_x(i) + enddo + hij = hij * phase + +end + diff --git a/src/determinants/slater_rules.irp.f b/src/determinants/slater_rules.irp.f index 52dfc143..7f0ccdbc 100644 --- a/src/determinants/slater_rules.irp.f +++ b/src/determinants/slater_rules.irp.f @@ -1581,9 +1581,12 @@ subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) end - - double precision function diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif use bitmasks implicit none BEGIN_DOC @@ -2292,3 +2295,592 @@ subroutine connected_to_hf(key_i,yes_no) yes_no = .True. endif end + + +!==============================================================================! +! ! +! Complex ! +! ! +!==============================================================================! + + +subroutine i_H_j_s2_complex(key_i,key_j,Nint,hij,s2) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC + ! Returns $\langle i|H|j \rangle$ and $\langle i|S^2|j \rangle$ + ! where $i$ and $j$ are determinants. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij, s2 + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_two_e_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase + integer :: n_occ_ab(2) + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + s2 = 0d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + ! Single alpha, single beta + if (exc(0,1,1) == 1) then + if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then + s2 = -phase + endif + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_two_e_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + ! Double alpha + else if (exc(0,1,1) == 2) then + hij = phase*(get_two_e_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_two_e_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + ! Double beta + else if (exc(0,1,2) == 2) then + hij = phase*(get_two_e_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_two_e_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_single_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + ! Single alpha + if (exc(0,1,1) == 1) then + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 + ! Single beta + else + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + endif + call get_single_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) + + case (0) + double precision, external :: diag_S_mat_elem + s2 = diag_S_mat_elem(key_i,Nint) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + + +subroutine i_H_j_complex(key_i,key_j,Nint,hij) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC + ! Returns $\langle i|H|j \rangle$ where $i$ and $j$ are determinants. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_two_e_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase + integer :: n_occ_ab(2) + PROVIDE mo_two_e_integrals_in_map mo_integrals_map big_array_exchange_integrals + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + + hij = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + integer :: spin + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Single alpha, single beta + if(exc(1,1,1) == exc(1,2,2) )then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) ==exc(1,1,2))then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_two_e_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_two_e_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_two_e_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_two_e_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_two_e_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_single_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + if (exc(0,1,1) == 1) then + ! Single alpha + m = exc(1,1,1) + p = exc(1,2,1) + spin = 1 + else + ! Single beta + m = exc(1,1,2) + p = exc(1,2,2) + spin = 2 + endif + call get_single_excitation_from_fock(key_i,key_j,p,m,spin,phase,hij) + + case (0) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + + + + +subroutine i_H_j_verbose_complex(key_i,key_j,Nint,hij,hmono,hdouble,phase) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC + ! Returns $\langle i|H|j \rangle$ where $i$ and $j$ are determinants. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij,hmono,hdouble,phase + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_two_e_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem + integer :: n_occ_ab(2) + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_two_e_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + hmono = 0.d0 + hdouble = 0.d0 + !DIR$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Single alpha, single beta + hij = phase*get_two_e_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_two_e_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_two_e_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_two_e_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_two_e_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_single_excitation(key_i,key_j,exc,phase,Nint) + !DIR$ FORCEINLINE + call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Single alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_two_e_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_two_e_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_two_e_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hdouble = hdouble + mipi(occ(k,1)) - miip(occ(k,1)) + enddo + do k = 1, elec_beta_num + hdouble = hdouble + mipi(occ(k,2)) + enddo + + else + ! Single beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_two_e_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_two_e_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_two_e_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hdouble = hdouble + mipi(occ(k,1)) + enddo + do k = 1, elec_beta_num + hdouble = hdouble + mipi(occ(k,2)) - miip(occ(k,2)) + enddo + + endif + hmono = mo_one_e_integrals(m,p) + hij = phase*(hdouble + hmono) + + case (0) + phase = 1.d0 + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + +subroutine i_H_psi_complex(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC +! Computes $\langle i|H|Psi \rangle = \sum_J c_J \langle i | H | J \rangle$. +! +! Uses filter_connected_i_H_psi0 to get all the $|J \rangle$ to which $|i \rangle$ +! is connected. +! The i_H_psi_minilist is much faster but requires to build the +! minilists. + END_DOC + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer, allocatable :: idx(:) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) + + i_H_psi_array = 0.d0 + + call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) + if (Nstate == 1) then + + do ii=1,idx(0) + i = idx(ii) + !DIR$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + i_H_psi_array(1) = i_H_psi_array(1) + coef(i,1)*hij + enddo + + else + + do ii=1,idx(0) + i = idx(ii) + !DIR$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo + + endif + +end + + +subroutine i_H_psi_minilist_complex(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate,idx_key(Ndet), N_minilist + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + + integer :: i, ii,j, i_in_key, i_in_coef + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer, allocatable :: idx(:) + BEGIN_DOC +! Computes $\langle i|H|\Psi \rangle = \sum_J c_J \langle i|H|J\rangle$. +! +! Uses filter_connected_i_H_psi0 to get all the $|J \rangle$ to which $|i \rangle$ +! is connected. The $|J\rangle$ are searched in short pre-computed lists. + END_DOC + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + allocate(idx(0:Ndet)) + i_H_psi_array = 0.d0 + + call filter_connected_i_H_psi0(keys,key,Nint,N_minilist,idx) + if (Nstate == 1) then + + do ii=1,idx(0) + i_in_key = idx(ii) + i_in_coef = idx_key(idx(ii)) + !DIR$ FORCEINLINE + call i_H_j(keys(1,1,i_in_key),key,Nint,hij) + ! TODO : Cache misses + i_H_psi_array(1) = i_H_psi_array(1) + coef(i_in_coef,1)*hij + enddo + + else + + do ii=1,idx(0) + i_in_key = idx(ii) + i_in_coef = idx_key(idx(ii)) + !DIR$ FORCEINLINE + call i_H_j(keys(1,1,i_in_key),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i_in_coef,j)*hij + enddo + enddo + + endif + +end + + + +subroutine i_H_j_single_spin_complex(key_i,key_j,Nint,spin,hij) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC + ! Returns $\langle i|H|j \rangle$ where $i$ and $j$ are determinants differing by + ! a single excitation. + END_DOC + integer, intent(in) :: Nint, spin + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2) + double precision :: phase + + PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + + call get_single_excitation_spin(key_i(1,spin),key_j(1,spin),exc,phase,Nint) + call get_single_excitation_from_fock(key_i,key_j,exc(1,1),exc(1,2),spin,phase,hij) +end + +subroutine i_H_j_double_spin_complex(key_i,key_j,Nint,hij) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC + ! Returns $\langle i|H|j \rangle$ where $i$ and $j$ are determinants differing by + ! a same-spin double excitation. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint), key_j(Nint) + double precision, intent(out) :: hij + + integer :: exc(0:2,2) + double precision :: phase + double precision, external :: get_two_e_integral + + PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + call get_double_excitation_spin(key_i,key_j,exc,phase,Nint) + hij = phase*(get_two_e_integral( & + exc(1,1), & + exc(2,1), & + exc(1,2), & + exc(2,2), mo_integrals_map) - & + get_two_e_integral( & + exc(1,1), & + exc(2,1), & + exc(2,2), & + exc(1,2), mo_integrals_map) ) +end + +subroutine i_H_j_double_alpha_beta_complex(key_i,key_j,Nint,hij) + !todo: modify/implement for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif + use bitmasks + implicit none + BEGIN_DOC + ! Returns $\langle i|H|j \rangle$ where $i$ and $j$ are determinants differing by + ! an opposite-spin double excitation. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + double precision :: phase, phase2 + double precision, external :: get_two_e_integral + + PROVIDE big_array_exchange_integrals mo_two_e_integrals_in_map + + call get_single_excitation_spin(key_i(1,1),key_j(1,1),exc(0,1,1),phase,Nint) + call get_single_excitation_spin(key_i(1,2),key_j(1,2),exc(0,1,2),phase2,Nint) + phase = phase*phase2 + if (exc(1,1,1) == exc(1,2,2)) then + hij = phase * big_array_exchange_integrals(exc(1,1,1),exc(1,1,2),exc(1,2,1)) + else if (exc(1,2,1) == exc(1,1,2)) then + hij = phase * big_array_exchange_integrals(exc(1,2,1),exc(1,1,1),exc(1,2,2)) + else + hij = phase*get_two_e_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + endif +end diff --git a/src/determinants/spindeterminants.ezfio_config b/src/determinants/spindeterminants.ezfio_config index 39ccb82b..bd4b80ce 100644 --- a/src/determinants/spindeterminants.ezfio_config +++ b/src/determinants/spindeterminants.ezfio_config @@ -10,6 +10,7 @@ spindeterminants psi_coef_matrix_rows integer (spindeterminants_n_det) psi_coef_matrix_columns integer (spindeterminants_n_det) psi_coef_matrix_values double precision (spindeterminants_n_det,spindeterminants_n_states) + psi_coef_matrix_values_complex double precision (2,spindeterminants_n_det,spindeterminants_n_states) n_svd_coefs integer psi_svd_alpha double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs,spindeterminants_n_states) psi_svd_beta double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs,spindeterminants_n_states) diff --git a/src/determinants/spindeterminants.irp.f b/src/determinants/spindeterminants.irp.f index 716c81ee..974ec614 100644 --- a/src/determinants/spindeterminants.irp.f +++ b/src/determinants/spindeterminants.irp.f @@ -307,8 +307,12 @@ integer function get_index_in_psi_det_beta_unique(key,Nint) end - subroutine write_spindeterminants + !todo: modify for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif use bitmasks implicit none integer(8), allocatable :: tmpdet(:,:) @@ -349,8 +353,12 @@ subroutine write_spindeterminants enddo call ezfio_set_spindeterminants_psi_det_beta(psi_det_beta_unique) deallocate(tmpdet) - + + if (is_complex) then + call ezfio_set_spindeterminants_psi_coef_matrix_values_complex(psi_bilinear_matrix_values_complex) + else call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_bilinear_matrix_values) + endif call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_bilinear_matrix_rows) call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_bilinear_matrix_columns) @@ -370,6 +378,18 @@ end det_alpha_norm = 0.d0 det_beta_norm = 0.d0 + if (is_complex) then + do k=1,N_det + i = psi_bilinear_matrix_rows(k) + j = psi_bilinear_matrix_columns(k) + f = 0.d0 + do l=1,N_states + f += cdabs(psi_bilinear_matrix_values_complex(k,l)*psi_bilinear_matrix_values_complex(k,l)) * state_average_weight(l) + enddo + det_alpha_norm(i) += f + det_beta_norm(j) += f + enddo + else do k=1,N_det i = psi_bilinear_matrix_rows(k) j = psi_bilinear_matrix_columns(k) @@ -380,6 +400,7 @@ end det_alpha_norm(i) += f det_beta_norm(j) += f enddo + endif det_alpha_norm = det_alpha_norm det_beta_norm = det_beta_norm @@ -392,8 +413,35 @@ END_PROVIDER ! ! !==============================================================================! - BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ] -&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows , (N_det) ] +BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_values, (N_det,N_states) ] + use bitmasks + PROVIDE psi_bilinear_matrix_rows + do k=1,N_det + do l=1,N_states + psi_bilinear_matrix_values(k,l) = psi_coef(k,l) + enddo + enddo + do l=1,N_states + call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det) + enddo +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, psi_bilinear_matrix_values_complex, (N_det,N_states) ] + use bitmasks + PROVIDE psi_bilinear_matrix_rows + do k=1,N_det + do l=1,N_states + psi_bilinear_matrix_values_complex(k,l) = psi_coef_complex(k,l) + enddo + enddo + do l=1,N_states + call cdset_order(psi_bilinear_matrix_values_complex(1,l),psi_bilinear_matrix_order,N_det) + enddo +END_PROVIDER + + + + BEGIN_PROVIDER [ integer, psi_bilinear_matrix_rows , (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns, (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_order , (N_det) ] use bitmasks @@ -408,10 +456,13 @@ END_PROVIDER END_DOC integer :: i,j,k, l integer(bit_kind) :: tmp_det(N_int,2) - integer, external :: get_index_in_psi_det_sorted_bit +! integer, external :: get_index_in_psi_det_sorted_bit - - PROVIDE psi_coef_sorted_bit + if (is_complex) then + PROVIDE psi_coef_sorted_bit_complex + else + PROVIDE psi_coef_sorted_bit + endif integer*8, allocatable :: to_sort(:) integer, external :: get_index_in_psi_det_alpha_unique @@ -427,9 +478,6 @@ END_PROVIDER ASSERT (j>0) ASSERT (j<=N_det_beta_unique) - do l=1,N_states - psi_bilinear_matrix_values(k,l) = psi_coef(k,l) - enddo psi_bilinear_matrix_rows(k) = i psi_bilinear_matrix_columns(k) = j to_sort(k) = int(N_det_alpha_unique,8) * int(j-1,8) + int(i,8) @@ -445,11 +493,6 @@ END_PROVIDER !$OMP SINGLE call iset_order(psi_bilinear_matrix_columns,psi_bilinear_matrix_order,N_det) !$OMP END SINGLE - !$OMP DO - do l=1,N_states - call dset_order(psi_bilinear_matrix_values(1,l),psi_bilinear_matrix_order,N_det) - enddo - !$OMP END DO !$OMP END PARALLEL deallocate(to_sort) ASSERT (minval(psi_bilinear_matrix_rows) == 1) @@ -514,8 +557,71 @@ BEGIN_PROVIDER [ integer, psi_bilinear_matrix_columns_loc, (N_det_beta_unique+1) END_PROVIDER - BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ] -&BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows , (N_det) ] +BEGIN_PROVIDER [ double precision, psi_bilinear_matrix_transp_values, (N_det,N_states) ] + use bitmasks + implicit none + BEGIN_DOC + ! Transpose of :c:data:`psi_bilinear_matrix` + ! + ! $D_\beta^\dagger.C^\dagger.D_\alpha$ + ! + ! Rows are $\alpha$ determinants and columns are $\beta$, but the matrix is stored in row major + ! format. + END_DOC + integer :: i,j,k,l + + PROVIDE psi_bilinear_matrix_transp_rows + + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l) + do l=1,N_states + !$OMP DO + do k=1,N_det + psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l) + enddo + !$OMP ENDDO NOWAIT + enddo + !$OMP END PARALLEL + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l) + do l=1,N_states + call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det) + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, psi_bilinear_matrix_transp_values_complex, (N_det,N_states) ] + use bitmasks + implicit none + BEGIN_DOC + ! Transpose of :c:data:`psi_bilinear_matrix` + ! + ! $D_\beta^\dagger.C^\dagger.D_\alpha$ + ! + ! Rows are $\alpha$ determinants and columns are $\beta$, but the matrix is stored in row major + ! format. + END_DOC + integer :: i,j,k,l + + PROVIDE psi_bilinear_matrix_transp_rows + + !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l) + do l=1,N_states + !$OMP DO + do k=1,N_det + psi_bilinear_matrix_transp_values_complex (k,l) = psi_bilinear_matrix_values_complex (k,l) + enddo + !$OMP ENDDO NOWAIT + enddo + !$OMP END PARALLEL + !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l) + do l=1,N_states + call cdset_order(psi_bilinear_matrix_transp_values_complex(1,l),psi_bilinear_matrix_transp_order,N_det) + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + + BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_rows , (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_columns, (N_det) ] &BEGIN_PROVIDER [ integer, psi_bilinear_matrix_transp_order , (N_det) ] use bitmasks @@ -530,18 +636,15 @@ END_PROVIDER END_DOC integer :: i,j,k,l - PROVIDE psi_coef_sorted_bit + if (is_complex) then + PROVIDE psi_coef_sorted_bit_complex + else + PROVIDE psi_coef_sorted_bit + endif integer*8, allocatable :: to_sort(:) allocate(to_sort(N_det)) !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,l) - do l=1,N_states - !$OMP DO - do k=1,N_det - psi_bilinear_matrix_transp_values (k,l) = psi_bilinear_matrix_values (k,l) - enddo - !$OMP ENDDO NOWAIT - enddo !$OMP DO do k=1,N_det psi_bilinear_matrix_transp_columns(k) = psi_bilinear_matrix_columns(k) @@ -563,11 +666,6 @@ END_PROVIDER call i8radix_sort(to_sort, psi_bilinear_matrix_transp_order, N_det,-1) call iset_order(psi_bilinear_matrix_transp_rows,psi_bilinear_matrix_transp_order,N_det) call iset_order(psi_bilinear_matrix_transp_columns,psi_bilinear_matrix_transp_order,N_det) - !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l) - do l=1,N_states - call dset_order(psi_bilinear_matrix_transp_values(1,l),psi_bilinear_matrix_transp_order,N_det) - enddo - !$OMP END PARALLEL DO deallocate(to_sort) ASSERT (minval(psi_bilinear_matrix_transp_columns) == 1) ASSERT (minval(psi_bilinear_matrix_transp_rows) == 1) @@ -641,7 +739,30 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_de enddo END_PROVIDER +BEGIN_PROVIDER [ complex*16, psi_bilinear_matrix_complex, (N_det_alpha_unique,N_det_beta_unique,N_states) ] + implicit none + BEGIN_DOC + ! Coefficient matrix if the wave function is expressed in a bilinear form : + ! + ! $D_\alpha^\dagger.C.D_\beta$ + END_DOC + integer :: i,j,k,istate + psi_bilinear_matrix_complex = (0.d0,0.d0) + do k=1,N_det + i = psi_bilinear_matrix_rows(k) + j = psi_bilinear_matrix_columns(k) + do istate=1,N_states + psi_bilinear_matrix_complex(i,j,istate) = psi_bilinear_matrix_values_complex(k,istate) + enddo + enddo +END_PROVIDER + subroutine create_wf_of_psi_bilinear_matrix(truncate) + !todo: modify for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif use bitmasks implicit none BEGIN_DOC @@ -713,6 +834,11 @@ subroutine create_wf_of_psi_bilinear_matrix(truncate) end subroutine generate_all_alpha_beta_det_products + !todo: modify for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif implicit none BEGIN_DOC ! Creates a wave function from all possible $\alpha \times \beta$ determinants @@ -856,6 +982,11 @@ end subroutine copy_psi_bilinear_to_psi(psi, isize) + !todo: modify for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif implicit none BEGIN_DOC ! Overwrites :c:data:`psi_det` and :c:data:`psi_coef` with the wave function @@ -1292,6 +1423,11 @@ END_TEMPLATE subroutine wf_of_psi_bilinear_matrix(truncate) + !todo: modify for complex + if (is_complex) then + print*,irp_here,' not implemented for complex' + stop -1 + endif use bitmasks implicit none BEGIN_DOC diff --git a/src/utils/sort.irp.f b/src/utils/sort.irp.f index ce609411..cf5e0038 100644 --- a/src/utils/sort.irp.f +++ b/src/utils/sort.irp.f @@ -409,6 +409,7 @@ BEGIN_TEMPLATE SUBST [ X, type ] ; real ;; d ; double precision ;; + cd; complex*16 ;; i ; integer ;; i8; integer*8 ;; i2; integer*2 ;; diff --git a/src/utils_complex/qp2-pbc-diff.txt b/src/utils_complex/qp2-pbc-diff.txt index 884c4341..ad07aeb4 100644 --- a/src/utils_complex/qp2-pbc-diff.txt +++ b/src/utils_complex/qp2-pbc-diff.txt @@ -3,14 +3,14 @@ current: general: - i_h_j_complex - diag_h_mat_elem if is_complex - + check for dependence on psi_det_sorted, clean up providers determinants: (done) connected_to_ref.irp.f (done) create_excitations.irp.f - (****) density_matrix.irp.f + (done?)density_matrix{,_complex}.irp.f + no one_e_dm_dagger_mo_spin_index_complex + need to test for complex (done) determinants_bitmasks.irp.f (****) determinants{_complex}.irp.f mostly done @@ -41,15 +41,24 @@ determinants: (done) psi_cas{,_complex}.irp.f might be able to combine some providers?? (done) psi_energy_mono_elec.irp.f - (****) ref_bitmask.irp.f + (done) ref_bitmask.irp.f (****) s2{,_complex}.irp.f - (****) single_excitations.irp.f + made copies of needed functions for complex + still need to do implementation + (done) single_excitations.irp.f (****) single_excitation_two_e.irp.f (****) slater_rules.irp.f + made copies of needed functions for complex + still need to do implementation (****) slater_rules_wee_mono.irp.f (done) sort_dets_ab.irp.f spindeterminants.ezfio_config - (****) spindeterminants.irp.f + need svd complex? + (done?) spindeterminants.irp.f + separated psi_bilinear_matrix_values from psi_bilinear_matrix_{rows,columns,order} + new provider for psi_bilinear_matrix_values_complex + same for bilinear matrix transp (no conjugate) + done except for specific functions that are commented with todo (****) two_e_density_matrix.irp.pouet (done) utils.irp.f (****) zmq.irp.f