diff --git a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f index 98ef0b49..f34242ab 100644 --- a/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f +++ b/plugins/Full_CI_ZMQ/pt2_stoch_routines.irp.f @@ -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 diff --git a/plugins/Full_CI_ZMQ/selection.irp.f b/plugins/Full_CI_ZMQ/selection.irp.f index 32c635ec..86e9e9f2 100644 --- a/plugins/Full_CI_ZMQ/selection.irp.f +++ b/plugins/Full_CI_ZMQ/selection.irp.f @@ -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) diff --git a/plugins/Psiref_CAS/psi_ref.irp.f b/plugins/Psiref_CAS/psi_ref.irp.f index ab9e6943..87439764 100644 --- a/plugins/Psiref_CAS/psi_ref.irp.f +++ b/plugins/Psiref_CAS/psi_ref.irp.f @@ -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 diff --git a/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f b/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f index b8b021e8..e4555d8c 100644 --- a/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f +++ b/plugins/mrsc2_no_amp/mrsc2_no_amp.irp.f @@ -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) ] diff --git a/plugins/mrsc2_no_amp/sc2_no_amp.irp.f b/plugins/mrsc2_no_amp/sc2_no_amp.irp.f index 622d7236..f557783b 100644 --- a/plugins/mrsc2_no_amp/sc2_no_amp.irp.f +++ b/plugins/mrsc2_no_amp/sc2_no_amp.irp.f @@ -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 diff --git a/src/Determinants/print_wf.irp.f b/src/Determinants/print_wf.irp.f index af109e2d..737e4d3e 100644 --- a/src/Determinants/print_wf.irp.f +++ b/src/Determinants/print_wf.irp.f @@ -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*,' = ',hij endif