From 1c838a30d6fcd36d3778f9f3ec7871b29e245965 Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Thu, 20 Feb 2020 14:56:47 -0600 Subject: [PATCH] working on complex determinants --- src/determinants/psi_energy_mono_elec.irp.f | 21 +- src/determinants/s2.irp.f | 4 + src/determinants/s2_complex.irp.f | 290 ++++++++++++++++++++ src/determinants/utils.irp.f | 22 ++ src/utils_complex/qp2-pbc-diff.txt | 68 +++-- 5 files changed, 386 insertions(+), 19 deletions(-) create mode 100644 src/determinants/s2_complex.irp.f diff --git a/src/determinants/psi_energy_mono_elec.irp.f b/src/determinants/psi_energy_mono_elec.irp.f index 74e69160..1f6d69d5 100644 --- a/src/determinants/psi_energy_mono_elec.irp.f +++ b/src/determinants/psi_energy_mono_elec.irp.f @@ -9,7 +9,26 @@ ! computed using the :c:data:`one_e_dm_mo_alpha` + ! :c:data:`one_e_dm_mo_beta` and :c:data:`mo_one_e_integrals` END_DOC + double precision :: accu psi_energy_h_core = 0.d0 + if (is_complex) then + do i = 1, N_states + do j = 1, mo_num + do k = 1, mo_num + psi_energy_h_core(i) += dble(mo_one_e_integrals_complex(k,j) * & + (one_e_dm_mo_alpha_complex(j,k,i) + one_e_dm_mo_beta_complex(j,k,i))) + enddo + enddo + enddo + do i = 1, N_states + accu = 0.d0 + do j = 1, mo_num + accu += dble(one_e_dm_mo_alpha_complex(j,j,i) + one_e_dm_mo_beta_complex(j,j,i)) + enddo + accu = (elec_alpha_num + elec_beta_num ) / accu + psi_energy_h_core(i) = psi_energy_h_core(i) * accu + enddo + else do i = 1, N_states do j = 1, mo_num do k = 1, mo_num @@ -17,7 +36,6 @@ enddo enddo enddo - double precision :: accu do i = 1, N_states accu = 0.d0 do j = 1, mo_num @@ -26,4 +44,5 @@ accu = (elec_alpha_num + elec_beta_num ) / accu psi_energy_h_core(i) = psi_energy_h_core(i) * accu enddo + endif END_PROVIDER diff --git a/src/determinants/s2.irp.f b/src/determinants/s2.irp.f index 391d0073..6f9560ca 100644 --- a/src/determinants/s2.irp.f +++ b/src/determinants/s2.irp.f @@ -98,7 +98,11 @@ BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] ! array of the averaged values of the S^2 operator on the various states END_DOC integer :: i + if (is_complex) then + call u_0_S2_u_0_complex(s2_values,psi_coef_complex,n_det,psi_det,N_int,N_states,psi_det_size) + else call u_0_S2_u_0(s2_values,psi_coef,n_det,psi_det,N_int,N_states,psi_det_size) + endif END_PROVIDER diff --git a/src/determinants/s2_complex.irp.f b/src/determinants/s2_complex.irp.f new file mode 100644 index 00000000..bb368f55 --- /dev/null +++ b/src/determinants/s2_complex.irp.f @@ -0,0 +1,290 @@ +subroutine u_0_S2_u_0_complex(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8) + print*,irp_here,' not implemented for complex' + stop -1 +! use bitmasks +! implicit none +! BEGIN_DOC +! ! Computes e_0 = / +! ! +! ! n : number of determinants +! ! +! END_DOC +! integer, intent(in) :: n,Nint, N_st, sze_8 +! double precision, intent(out) :: e_0(N_st) +! double precision, intent(in) :: u_0(sze_8,N_st) +! integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) +! +! double precision, allocatable :: v_0(:,:) +! double precision :: u_dot_u,u_dot_v +! integer :: i,j +! allocate (v_0(sze_8,N_st)) +! +! call S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8) +! do i=1,N_st +! e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n) + S_z2_Sz +! enddo +end + + + +subroutine S2_u_0_complex(v_0,u_0,n,keys_tmp,Nint) + print*,irp_here,' not implemented for complex' + stop -1 +! use bitmasks +! implicit none +! BEGIN_DOC +! ! Computes v_0 = S^2|u_0> +! ! +! ! n : number of determinants +! ! +! END_DOC +! integer, intent(in) :: n,Nint +! double precision, intent(out) :: v_0(n) +! double precision, intent(in) :: u_0(n) +! integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) +! call S2_u_0_nstates(v_0,u_0,n,keys_tmp,Nint,1,n) +end + +subroutine S2_u_0_nstates_complex(v_0,u_0,n,keys_tmp,Nint,N_st,sze_8) + print*,irp_here,' not implemented for complex' + stop -1 +! use bitmasks +! implicit none +! BEGIN_DOC +! ! Computes v_0 = S^2|u_0> +! ! +! ! n : number of determinants +! ! +! END_DOC +! integer, intent(in) :: N_st,n,Nint, sze_8 +! double precision, intent(out) :: v_0(sze_8,N_st) +! double precision, intent(in) :: u_0(sze_8,N_st) +! integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) +! double precision :: s2_tmp +! double precision, allocatable :: vt(:,:) +! integer :: i,j,k,l, jj,ii +! integer :: i0, j0 +! +! integer, allocatable :: shortcut(:,:), sort_idx(:,:) +! integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:) +! integer(bit_kind) :: sorted_i(Nint) +! +! integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate +! +! +! ASSERT (Nint > 0) +! ASSERT (Nint == N_int) +! ASSERT (n>0) +! PROVIDE ref_bitmask_energy +! +! allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2)) +! v_0 = 0.d0 +! +! call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint) +! call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint) +! +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE(i,s2_tmp,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i,istate)& +! !$OMP SHARED(n,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version,N_st,sze_8) +! allocate(vt(sze_8,N_st)) +! vt = 0.d0 +! +! do sh=1,shortcut(0,1) +! !$OMP DO SCHEDULE(static,1) +! do sh2=sh,shortcut(0,1) +! exa = 0 +! do ni=1,Nint +! exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1))) +! end do +! if(exa > 2) then +! cycle +! end if +! +! do i=shortcut(sh,1),shortcut(sh+1,1)-1 +! org_i = sort_idx(i,1) +! if(sh==sh2) then +! endi = i-1 +! else +! endi = shortcut(sh2+1,1)-1 +! end if +! do ni=1,Nint +! sorted_i(ni) = sorted(ni,i,1) +! enddo +! +! do j=shortcut(sh2,1),endi +! org_j = sort_idx(j,1) +! ext = exa +! do ni=1,Nint +! ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1))) +! end do +! if(ext <= 4) then +! call get_s2(keys_tmp(1,1,org_i),keys_tmp(1,1,org_j),Nint,s2_tmp) +! do istate=1,N_st +! vt (org_i,istate) = vt (org_i,istate) + s2_tmp*u_0(org_j,istate) +! vt (org_j,istate) = vt (org_j,istate) + s2_tmp*u_0(org_i,istate) +! enddo +! endif +! enddo +! enddo +! enddo +! !$OMP END DO NOWAIT +! enddo +! +! do sh=1,shortcut(0,2) +! !$OMP DO +! do i=shortcut(sh,2),shortcut(sh+1,2)-1 +! org_i = sort_idx(i,2) +! do j=shortcut(sh,2),i-1 +! org_j = sort_idx(j,2) +! ext = 0 +! do ni=1,Nint +! ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2))) +! end do +! if(ext == 4) then +! call get_s2(keys_tmp(1,1,org_i),keys_tmp(1,1,org_j),Nint,s2_tmp) +! do istate=1,N_st +! vt (org_i,istate) = vt (org_i,istate) + s2_tmp*u_0(org_j,istate) +! vt (org_j,istate) = vt (org_j,istate) + s2_tmp*u_0(org_i,istate) +! enddo +! end if +! end do +! end do +! !$OMP END DO NOWAIT +! enddo +! !$OMP BARRIER +! +! do istate=1,N_st +! do i=n,1,-1 +! !$OMP ATOMIC +! v_0(i,istate) = v_0(i,istate) + vt(i,istate) +! enddo +! enddo +! +! deallocate(vt) +! !$OMP END PARALLEL +! +! do i=1,n +! call get_s2(keys_tmp(1,1,i),keys_tmp(1,1,i),Nint,s2_tmp) +! do istate=1,N_st +! v_0(i,istate) += s2_tmp * u_0(i,istate) +! enddo +! enddo +! +! deallocate (shortcut, sort_idx, sorted, version) +end + + + + + + + +subroutine get_uJ_s2_uI_complex(psi_keys_tmp,psi_coefs_tmp,n,nmax_coefs,nmax_keys,s2,nstates) + print*,irp_here,' not implemented for complex' + stop -1 +! implicit none +! use bitmasks +! integer, intent(in) :: n,nmax_coefs,nmax_keys,nstates +! integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax_keys) +! double precision, intent(in) :: psi_coefs_tmp(nmax_coefs,nstates) +! double precision, intent(out) :: s2(nstates,nstates) +! double precision :: s2_tmp,accu +! integer :: i,j,l,jj,ll,kk +! integer, allocatable :: idx(:) +! BEGIN_DOC +! ! returns the matrix elements of S^2 "s2(i,j)" between the "nstates" states +! ! psi_coefs_tmp(:,i) and psi_coefs_tmp(:,j) +! END_DOC +! s2 = 0.d0 +! do ll = 1, nstates +! do jj = 1, nstates +! accu = 0.d0 +! !$OMP PARALLEL DEFAULT(NONE) & +! !$OMP PRIVATE (i,j,kk,idx,s2_tmp) & +! !$OMP SHARED (ll,jj,psi_keys_tmp,psi_coefs_tmp,N_int,n,nstates)& +! !$OMP REDUCTION(+:accu) +! allocate(idx(0:n)) +! !$OMP DO SCHEDULE(dynamic) +! do i = n,1,-1 ! Better OMP scheduling +! call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),N_int,s2_tmp) +! accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(i,jj) +! call filter_connected(psi_keys_tmp,psi_keys_tmp(1,1,i),N_int,i-1,idx) +! do kk=1,idx(0) +! j = idx(kk) +! call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),N_int,s2_tmp) +! accu += psi_coefs_tmp(i,ll) * s2_tmp * psi_coefs_tmp(j,jj) + psi_coefs_tmp(i,jj) * s2_tmp * psi_coefs_tmp(j,ll) +! enddo +! enddo +! !$OMP END DO +! deallocate(idx) +! !$OMP END PARALLEL +! s2(ll,jj) += accu +! enddo +! enddo +! do i = 1, nstates +! do j =i+1,nstates +! accu = 0.5d0 * (s2(i,j) + s2(j,i)) +! s2(i,j) = accu +! s2(j,i) = accu +! enddo +! enddo +end + + +subroutine i_S2_psi_minilist_complex(key,keys,idx_key,N_minilist,coef,Nint,Ndet,Ndet_max,Nstate,i_S2_psi_array) + print*,irp_here,' not implemented for complex' + stop -1 +! 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_S2_psi_array(Nstate) +! +! integer :: i, ii,j, i_in_key, i_in_coef +! double precision :: phase +! integer :: exc(0:2,2,2) +! double precision :: s2ij +! integer :: idx(0:Ndet) +! BEGIN_DOC +!! Computes $\langle i|S^2|\Psi \rangle = \sum_J c_J \langle i|S^2|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) +! i_S2_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 get_s2(keys(1,1,i_in_key),key,Nint,s2ij) +! ! TODO : Cache misses +! i_S2_psi_array(1) = i_S2_psi_array(1) + coef(i_in_coef,1)*s2ij +! enddo +! +! else +! +! do ii=1,idx(0) +! i_in_key = idx(ii) +! i_in_coef = idx_key(idx(ii)) +! !DIR$ FORCEINLINE +! call get_s2(keys(1,1,i_in_key),key,Nint,s2ij) +! do j = 1, Nstate +! i_S2_psi_array(j) = i_S2_psi_array(j) + coef(i_in_coef,j)*s2ij +! enddo +! enddo +! +! endif +! +end diff --git a/src/determinants/utils.irp.f b/src/determinants/utils.irp.f index 3aec16f9..97258216 100644 --- a/src/determinants/utils.irp.f +++ b/src/determinants/utils.irp.f @@ -20,6 +20,28 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] !$OMP END PARALLEL DO END_PROVIDER +BEGIN_PROVIDER [ complex*16, h_matrix_all_dets_complex,(N_det,N_det) ] + use bitmasks + implicit none + BEGIN_DOC + ! |H| matrix on the basis of the Slater determinants defined by psi_det + END_DOC + integer :: i,j,k + complex*16 :: hij + integer :: degree(N_det),idx(0:N_det) + call i_h_j_complex(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + !$OMP PARALLEL DO SCHEDULE(GUIDED) DEFAULT(NONE) PRIVATE(i,j,hij,degree,idx,k) & + !$OMP SHARED (N_det, psi_det, N_int,h_matrix_all_dets_complex) + do i =1,N_det + do j = i, N_det + call i_h_j_complex(psi_det(1,1,i),psi_det(1,1,j),N_int,hij) + H_matrix_all_dets_complex(i,j) = hij + H_matrix_all_dets_complex(j,i) = dconjg(hij) + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + BEGIN_PROVIDER [ double precision, S2_matrix_all_dets,(N_det,N_det) ] use bitmasks diff --git a/src/utils_complex/qp2-pbc-diff.txt b/src/utils_complex/qp2-pbc-diff.txt index 06d74811..884c4341 100644 --- a/src/utils_complex/qp2-pbc-diff.txt +++ b/src/utils_complex/qp2-pbc-diff.txt @@ -2,25 +2,57 @@ ------------------------------------------------------------------------------------- current: -determinants: - TODO - create_excitations - do_single_excitation - use symmetry rules to simplify? - should this be general, or should we only allow singles that conserve momentum? - density_matrix - determinants - ezfio_set_determinants_psi_coef_complex_qp_edit? (need ocaml?) - psi_coef_{max,min}? - save_wavefunction_specified{,_complex} qp_edit save? - psi_energy_mono_elec - diag_h_mat_elem for complex - ... +general: + i_h_j_complex + diag_h_mat_elem if is_complex - DONE - create_excitations - build_singly_excited_wavefunction{_complex} - ... + +determinants: + (done) connected_to_ref.irp.f + (done) create_excitations.irp.f + (****) density_matrix.irp.f + (done) determinants_bitmasks.irp.f + (****) determinants{_complex}.irp.f + mostly done + could separate/combine some providers instead of copying + for psi_{det,coef}_sorted: + use same linked provider for psi_average_norm_contrib_sorted + psi_det_sorted_order + psi_det_sorted + different providers for psi_coef{,_complex}_sorted + need to figure out {,abs_}psi_coef_{min,max} + need to modify ocaml for psi_coef_complex_qp_edit? + save_wavefunction_specified? qp_edit save? (wrong for real?) + (done) energy.irp.f + needs diag_h_mat_elem function to be modified for complex + (????) example.irp.f + (****) EZFIO.cfg + (done) filter_connected.irp.f + (done) fock_diag.irp.f + (****) h_apply.irp.f + (****) h_apply_nozmq.template.f + (****) h_apply.template.f + (****) h_apply_zmq.template.f + (****) occ_pattern.irp.f + mostly done? + might need to change calls to fill_h_apply_buffer_no_selection? + check again after modifying h_apply for complex + (done) prune_wf.irp.f + (done) psi_cas{,_complex}.irp.f + might be able to combine some providers?? + (done) psi_energy_mono_elec.irp.f + (****) ref_bitmask.irp.f + (****) s2{,_complex}.irp.f + (****) single_excitations.irp.f + (****) single_excitation_two_e.irp.f + (****) slater_rules.irp.f + (****) slater_rules_wee_mono.irp.f + (done) sort_dets_ab.irp.f + spindeterminants.ezfio_config + (****) spindeterminants.irp.f + (****) two_e_density_matrix.irp.pouet + (done) utils.irp.f + (****) zmq.irp.f -------------------------------------------------------------------------------------