LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-25 22:52:15 +02:00
Code Issues Releases Wiki Activity

Repaired selection

Browse Source
This commit is contained in:
Anthony Scemama 2017-01-30 09:38:04 +01:00
parent b08ced8741
commit 097083db47
6 changed files with 109 additions and 63 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f
View File

@ -1,7 +1,5 @@
BEGIN_PROVIDER [ integer, fragment_count ]
&BEGIN_PROVIDER [ integer, fragment_first ]
fragment_count = (elec_alpha_num-n_core_orb)**2
BEGIN_PROVIDER [ integer, fragment_first ]
implicit none
fragment_first = first_det_of_teeth(1)
END_PROVIDER
@ -111,7 +109,7 @@ subroutine do_carlo(tbc, Ncomb, comb, pt2_detail, computed, sumabove, sum2above,
myVal = 0d0
myVal2 = 0d0
do j=comb_teeth,1,-1
myVal += pt2_detail(1, dets(j)) / weight(dets(j)) * comb_step
myVal += pt2_detail(1, dets(j)) / pt2_weight(dets(j)) * comb_step
sumabove(j) += myVal
sum2above(j) += myVal**2
Nabove(j) += 1
@ -229,8 +227,8 @@ subroutine pt2_collector(b, tbc, comb, Ncomb, computed, pt2_detail, sumabove, su
end do
E0 = sum(pt2_detail(1,:first_det_of_teeth(tooth)-1))
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - cweight(first_det_of_teeth(tooth)-1))
prop = prop / weight(first_det_of_teeth(tooth))
prop = ((1d0 - dfloat(comb_teeth - tooth + 1) * comb_step) - pt2_cweight(first_det_of_teeth(tooth)-1))
prop = prop / pt2_weight(first_det_of_teeth(tooth))
E0 += pt2_detail(1,first_det_of_teeth(tooth)) * prop
avg = E0 + (sumabove(tooth) / Nabove(tooth))
eqt = sqrt(1d0 / (Nabove(tooth)-1) * abs(sum2above(tooth) / Nabove(tooth) - (sumabove(tooth)/Nabove(tooth))**2))
@ -393,7 +391,7 @@ subroutine get_comb(stato, dets)
curs = 1d0 - stato
do j = comb_teeth, 1, -1
dets(j) = pt2_find(curs, cweight)
dets(j) = pt2_find(curs, pt2_cweight)
curs -= comb_step
end do
end subroutine
@ -421,8 +419,8 @@ end subroutine
BEGIN_PROVIDER [ double precision, weight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, cweight, (N_det_generators) ]
BEGIN_PROVIDER [ double precision, pt2_weight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, pt2_cweight, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, comb_workload, (N_det_generators) ]
&BEGIN_PROVIDER [ double precision, comb_step ]
&BEGIN_PROVIDER [ integer, first_det_of_teeth, (comb_teeth+1) ]
@ -432,34 +430,34 @@ end subroutine
double precision :: norm_left, stato
integer, external :: pt2_find
weight(1) = psi_coef_generators(1,1)**2
cweight(1) = psi_coef_generators(1,1)**2
pt2_weight(1) = psi_coef_generators(1,1)**2
pt2_cweight(1) = psi_coef_generators(1,1)**2
do i=2,N_det_generators
weight(i) = psi_coef_generators(i,1)**2
cweight(i) = cweight(i-1) + psi_coef_generators(i,1)**2
pt2_weight(i) = psi_coef_generators(i,1)**2
pt2_cweight(i) = pt2_cweight(i-1) + psi_coef_generators(i,1)**2
end do
weight = weight / cweight(N_det_generators)
cweight = cweight / cweight(N_det_generators)
pt2_weight = pt2_weight / pt2_cweight(N_det_generators)
pt2_cweight = pt2_cweight / pt2_cweight(N_det_generators)
comb_workload = 1d0 / dfloat(N_det_generators)
norm_left = 1d0
comb_step = 1d0/dfloat(comb_teeth)
do i=1,N_det_generators
if(weight(i)/norm_left < comb_step/2d0) then
if(pt2_weight(i)/norm_left < comb_step/2d0) then
first_det_of_comb = i
exit
end if
norm_left -= weight(i)
norm_left -= pt2_weight(i)
end do
comb_step = 1d0 / dfloat(comb_teeth) * (1d0 - cweight(first_det_of_comb-1))
comb_step = 1d0 / dfloat(comb_teeth) * (1d0 - pt2_cweight(first_det_of_comb-1))
stato = 1d0 - comb_step! + 1d-5
do i=comb_teeth, 1, -1
first_det_of_teeth(i) = pt2_find(stato, cweight)
first_det_of_teeth(i) = pt2_find(stato, pt2_cweight)
stato -= comb_step
end do
first_det_of_teeth(comb_teeth+1) = N_det_generators + 1

19
plugins/Full_CI_ZMQ/selection.irp.f
View File

@ -1,5 +1,10 @@
use bitmasks
BEGIN_PROVIDER [ integer, fragment_count ]
implicit none
fragment_count = (elec_alpha_num-n_core_orb)**2
END_PROVIDER
double precision function integral8(i,j,k,l)
implicit none
@ -356,20 +361,6 @@ subroutine select_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,p
integer :: nb_count
do s1=1,2
do i1=N_holes(s1),1,-1 ! Generate low excitations first
! will_compute = (subset == 0)
! nb_count = 0
! if (s1==1) then
! nb_count = N_holes(1)-i1 + N_holes(2)
! else
! nb_count = N_holes(2)-i1
! endif
! maskInd = 12345
! fragment_count = 400
! subset = 3
! nb_count = 100
! if( nb_count >= (fragment_count - mod(maskInd+1, fragment_count) + subset-1) ) then
! will_compute = .true.
! end if
h1 = hole_list(i1,s1)
call apply_hole(psi_det_generators(1,1,i_generator), s1,h1, pmask, ok, N_int)

1
plugins/Psiref_CAS/psi_ref.irp.f
View File

@ -77,6 +77,7 @@ END_PROVIDER
norm_psi_ref(j) += psi_ref_coef(i,j) * psi_ref_coef(i,j)
enddo
inv_norm_psi_ref(j) = 1.d0/(dsqrt(norm_psi_Ref(j)))
print *, inv_norm_psi_ref(j)
enddo
END_PROVIDER

63
plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f
View File

@ -4,44 +4,95 @@
implicit none
integer :: i,j,k,l
integer, allocatable :: idx(:)
integer, allocatable :: holes_part(:,:)
double precision, allocatable :: e_corr(:,:)
double precision, allocatable :: accu(:)
double precision, allocatable :: ihpsi_current(:)
double precision, allocatable :: H_jj(:),H_jj_total(:),S2_jj(:)
integer :: number_of_particles, number_of_holes, n_h,n_p
allocate(e_corr(N_det_non_ref,N_states),ihpsi_current(N_states),accu(N_states),H_jj(N_det_non_ref),idx(0:N_det_non_ref))
allocate(H_jj_total(N_det),S2_jj(N_det))
allocate(holes_part(N_det,2))
accu = 0.d0
do i = 1, N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_interm_norm, N_int, N_det_ref,&
holes_part(i,1) = number_of_holes(psi_non_ref(1,1,i))
holes_part(i,2) = number_of_particles(psi_non_ref(1,1,i))
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,&
size(psi_ref_coef_interm_norm,1), N_states,ihpsi_current)
do j = 1, N_states
e_corr(i,j) = psi_non_ref_coef_interm_norm(i,j) * ihpsi_current(j)
e_corr(i,j) = psi_non_ref_coef(i,j) * ihpsi_current(j) * inv_norm_psi_ref(j)
accu(j) += e_corr(i,j)
enddo
enddo
print *, 'accu = ',accu
double precision :: hjj,diag_h_mat_elem
do i = 1, N_det_non_ref
call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx)
H_jj(i) = 0.d0
n_h = holes_part(i,1)
n_p = holes_part(i,2)
integer :: degree
! do j = 1, N_det_non_ref
! call get_excitation_degree(psi_non_ref(1,1,i),psi_non_ref(1,1,j),degree,N_int)
! if(degree .gt. 2)then
! if(n_h + holes_part(j,1) .gt. 2 .or. n_p + holes_part(j,2) .gt. 2 ) then
! H_jj(i) += e_corr(j,1)
! endif
! endif
! enddo
call filter_not_connected(psi_non_ref,psi_non_ref(1,1,i),N_int,N_det_non_ref,idx)
do j = 1, idx(0)
H_jj(i) += e_corr(idx(j),1)
if(n_h + holes_part(idx(j),1) .gt. 2 .or. n_p + holes_part(idx(j),2) .gt. 2 ) then
H_jj(i) += e_corr(idx(j),1)
endif
enddo
enddo
do i=1,N_Det
H_jj_total(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
call get_s2(psi_det(1,1,i),psi_det(1,1,i),N_int,S2_jj(i))
enddo
do i=1, N_det_non_ref
do i = 1, N_det_non_ref
H_jj_total(idx_non_ref(i)) += H_jj(i)
enddo
print *, 'coef'
call davidson_diag_hjj_sjj(psi_det,CI_eigenvectors_sc2_no_amp,H_jj_total,S2_jj,CI_electronic_energy_sc2_no_amp,size(CI_eigenvectors_sc2_no_amp,1),N_Det,N_states,N_states_diag,N_int,6)
do i = 1, N_det
hjj = diag_h_mat_elem(psi_det(1,1,i),N_int)
! if(hjj<-210.d0)then
! call debug_det(psi_det(1,1,i),N_int)
! print *, CI_eigenvectors_sc2_no_amp((i),1),hjj, H_jj_total(i)
! endif
enddo
print *, 'ref',N_det_ref
do i =1, N_det_ref
call debug_det(psi_det(1,1,idx_ref(i)),N_int)
print *, CI_eigenvectors_sc2_no_amp(idx_ref(i),1), H_jj_total(idx_ref(i))
enddo
print *, 'non ref',N_det_non_ref
do i=1, N_det_non_ref
hjj = diag_h_mat_elem(psi_non_ref(1,1,i),N_int)
! print *, CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1),H_jj_total(idx_non_ref(i)), H_jj(i)
! if(dabs(CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1)).gt.1.d-1)then
! if(hjj<-210.d0)then
! call debug_det(psi_det(1,1,idx_non_ref(i)),N_int)
! write(*,'(10(F16.10,X))') CI_eigenvectors_sc2_no_amp(idx_non_ref(i),1),hjj, H_jj(i),H_jj_total(idx_non_ref(i))
! endif
enddo
! do i = 1, N_det
! print *, CI_eigenvectors_sc2_no_amp(i,1)
! enddo
do i=1,N_states_diag
CI_eigenvectors_s2_sc2_no_amp(i) = S2_jj(i)
enddo
deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj)
deallocate(e_corr,ihpsi_current,accu,H_jj,idx,H_jj_total,s2_jj,holes_part)
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy_sc2_no_amp, (N_states_diag) ]

5
plugins/mrsc2_no_amp/sc2_no_amp.irp.f
View File

@ -1,9 +1,14 @@
program pouet
provide ao_bielec_integrals_in_map
call bla
end
subroutine bla
implicit none
integer :: i
do i = 1, 10
call diagonalize_CI_sc2_no_amp
TOUCH psi_coef
enddo
print *, "E+PT2 = ", ci_energy_sc2_no_amp(:)
end

46
src/Determinants/print_wf.irp.f
View File

@ -28,32 +28,32 @@ subroutine routine
if(degree == 0)then
print*,'Reference determinant '
else
call i_H_j(psi_det(1,1,i),psi_det(1,1,1),N_int,hij)
call i_H_j(psi_det(1,1,i),psi_det(1,1,i),N_int,hij)
call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
print*,'phase = ',phase
if(degree == 1)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
if(s1 == 1)then
norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
else
norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
endif
print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map)
double precision :: hmono,hdouble
call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
print*,'hmono = ',hmono
print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij
else
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2
print*,'h2,p2 = ',h2,p2
print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
endif
! if(degree == 1)then
! print*,'s1',s1
! print*,'h1,p1 = ',h1,p1
! if(s1 == 1)then
! norm_mono_a += dabs(psi_coef(i,1)/psi_coef(1,1))
! else
! norm_mono_b += dabs(psi_coef(i,1)/psi_coef(1,1))
! endif
! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,list_act(1),list_act(1),p1,mo_integrals_map)
! double precision :: hmono,hdouble
! call i_H_j_verbose(psi_det(1,1,1),psi_det(1,1,i),N_int,hij,hmono,hdouble)
! print*,'hmono = ',hmono
! print*,'hdouble = ',hdouble
! print*,'hmono+hdouble = ',hmono+hdouble
! print*,'hij = ',hij
! else
! print*,'s1',s1
! print*,'h1,p1 = ',h1,p1
! print*,'s2',s2
! print*,'h2,p2 = ',h2,p2
! print*,'< h | Ka| p > = ',get_mo_bielec_integral(h1,h2,p1,p2,mo_integrals_map)
! endif
print*,'<Ref| H |D_I> = ',hij
endif