LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 03:53:29 +01:00
Code Releases Activity

working on complex determinants

Browse Source
This commit is contained in:
Kevin Gasperich 2020-02-20 14:56:47 -06:00
parent 5c66e4b99f
commit 1c838a30d6
5 changed files with 386 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

21
src/determinants/psi_energy_mono_elec.irp.f
View File

@ -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

4
src/determinants/s2.irp.f
View File

@ -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

290
src/determinants/s2_complex.irp.f Normal file
View File

@ -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 = <u_0|S2|u_0>/<u_0|u_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

22
src/determinants/utils.irp.f
View File

@ -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

68
src/utils_complex/qp2-pbc-diff.txt
View File

@ -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
-------------------------------------------------------------------------------------