10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-17 18:55:26 +02:00

Merge branch 'master' into develop

This commit is contained in:
Anthony Scemama 2016-04-04 11:03:31 +02:00
commit a871bc451b
14 changed files with 555 additions and 125 deletions

View File

@ -7,11 +7,11 @@ interface: ezfio
type: Threshold
doc: Threshold on the convergence of the MRCC energy
interface: ezfio,provider,ocaml
default: 1.e-7
default: 1.e-5
[n_it_mrcc_max]
type: Strictly_positive_int
doc: Maximum number of MRCC iterations
interface: ezfio,provider,ocaml
default: 20
default: 10

View File

@ -1,13 +1,100 @@
program mrcc
implicit none
if (.not.read_wf) then
print *, 'read_wf has to be true.'
stop 1
endif
double precision, allocatable :: energy(:)
allocate (energy(N_states))
read_wf = .True.
SOFT_TOUCH read_wf
call print_cas_coefs
call run_mrcc
call set_generators_bitmasks_as_holes_and_particles
call run(N_states,energy)
if(do_pt2_end)then
call run_pt2(N_states,energy)
endif
deallocate(energy)
end
subroutine run(N_st,energy)
implicit none
integer, intent(in) :: N_st
double precision, intent(out) :: energy(N_st)
integer :: i
double precision :: E_new, E_old, delta_e
integer :: iteration
double precision :: E_past(4), lambda
E_new = 0.d0
delta_E = 1.d0
iteration = 0
lambda = 1.d0
do while (delta_E > thresh_mrcc)
iteration += 1
print *, '==========================='
print *, 'MRCC Iteration', iteration
print *, '==========================='
print *, ''
E_old = sum(ci_energy_dressed)
call write_double(6,ci_energy_dressed(1),"MRCC energy")
call diagonalize_ci_dressed(lambda)
E_new = sum(ci_energy_dressed)
delta_E = dabs(E_new - E_old)
call save_wavefunction
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
if (iteration > n_it_mrcc_max) then
exit
endif
enddo
call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
energy(:) = ci_energy_dressed(:)
end
subroutine run_pt2(N_st,energy)
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer, intent(in) :: N_st
double precision, intent(in) :: energy(N_st)
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
pt2 = 0.d0
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
N_det_generators = lambda_mrcc_pt2(0)
do i=1,N_det_generators
j = lambda_mrcc_pt2(i)
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
enddo
enddo
SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', energy
print *, 'E+PT2 = ', energy+pt2
print *, '-----'
call ezfio_set_full_ci_energy_pt2(energy+pt2)
deallocate(pt2,norm_pert)
end
subroutine print_cas_coefs
implicit none
@ -18,7 +105,7 @@ subroutine print_cas_coefs
print *, psi_cas_coef(i,:)
call debug_det(psi_cas(1,1,i),N_int)
enddo
call write_double(6,ci_energy(1),"Initial CI energy")
end

View File

@ -0,0 +1,91 @@
program mrcc_noiter
implicit none
double precision, allocatable :: energy(:)
allocate (energy(N_states))
read_wf = .True.
threshold_generators = .9999d0
SOFT_TOUCH read_wf threshold_generators
call print_cas_coefs
call set_generators_bitmasks_as_holes_and_particles
call run(N_states,energy)
if(do_pt2_end)then
call run_pt2(N_states,energy)
endif
deallocate(energy)
end
subroutine run(N_st,energy)
implicit none
integer, intent(in) :: N_st
double precision, intent(out) :: energy(N_st)
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
enddo
enddo
SOFT_TOUCH psi_coef ci_energy_dressed
call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
call save_wavefunction
energy(:) = ci_energy_dressed(:)
end
subroutine run_pt2(N_st,energy)
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer, intent(in) :: N_st
double precision, intent(in) :: energy(N_st)
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
pt2 = 0.d0
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
N_det_generators = lambda_mrcc_pt2(0)
do i=1,N_det_generators
j = lambda_mrcc_pt2(i)
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
enddo
enddo
SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', energy
print *, 'E+PT2 = ', energy+pt2
print *, '-----'
call ezfio_set_full_ci_energy_pt2(energy+pt2)
deallocate(pt2,norm_pert)
end
subroutine print_cas_coefs
implicit none
integer :: i,j
print *, 'CAS'
print *, '==='
do i=1,N_det_cas
print *, psi_cas_coef(i,:)
call debug_det(psi_cas(1,1,i),N_int)
enddo
call write_double(6,ci_energy(1),"Initial CI energy")
end

View File

@ -24,5 +24,12 @@ s.data["size_max"] = "3072"
print s
s = H_apply_zmq("mrcc_PT2")
s.energy = "ci_electronic_energy_dressed"
s.set_perturbation("epstein_nesbet_2x2")
s.unset_openmp()
print s
END_SHELL

View File

@ -189,10 +189,10 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
! Davidson iterations
! ===================
converged = .False.
integer :: iteration
converged = .False.
do while (.not.converged)
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st)
do k=1,N_st
@ -206,6 +206,7 @@ subroutine davidson_diag_hjj_mrcc(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nin
do iter=1,davidson_sze_max-1
! Compute W_k = H |u_k>
! ----------------------

View File

@ -24,7 +24,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
double precision, intent(inout) :: delta_ii_(Nstates,Ndet_ref)
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,l
integer :: i,j,k,l,m
integer :: degree_alpha(psi_det_size)
integer :: idx_alpha(0:psi_det_size)
logical :: good, fullMatch
@ -48,9 +48,14 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
integer :: N_miniList, ni, leng
double precision, allocatable :: hij_cache(:)
integer(bit_kind), allocatable :: microlist(:,:,:), microlist_zero(:,:,:)
integer, allocatable :: idx_microlist(:), N_microlist(:), ptr_microlist(:), idx_microlist_zero(:)
integer :: mobiles(2), smallerlist
leng = max(N_det_generators, N_det_non_ref)
allocate(miniList(Nint, 2, leng), idx_miniList(leng), hij_cache(N_det_non_ref))
allocate(miniList(Nint, 2, leng), idx_minilist(leng), hij_cache(N_det_non_ref))
!create_minilist_find_previous(key_mask, fullList, miniList, N_fullList, N_miniList, fullMatch, Nint)
call create_minilist_find_previous(key_mask, psi_det_generators, miniList, i_generator-1, N_miniList, fullMatch, Nint)
@ -59,8 +64,21 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
return
end if
allocate(ptr_microlist(0:mo_tot_num*2+1), &
N_microlist(0:mo_tot_num*2) )
allocate( microlist(Nint,2,N_minilist*4), &
idx_microlist(N_minilist*4))
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
if(key_mask(1,1) /= 0) then
call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
call find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
else
call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_minilist)
end if
deallocate(microlist, idx_microlist)
allocate (dIa_hla(Nstates,Ndet_non_ref))
@ -69,17 +87,101 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
! |alpha>
if(N_tq > 0) then
call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
call create_minilist(key_mask, psi_non_ref, miniList, idx_minilist, N_det_non_ref, N_minilist, Nint)
if(N_minilist == 0) return
if(key_mask(1,1) /= 0) then !!!!!!!!!!! PAS GENERAL !!!!!!!!!
allocate(microlist_zero(Nint,2,N_minilist), idx_microlist_zero(N_minilist))
allocate( microlist(Nint,2,N_minilist*4), &
idx_microlist(N_minilist*4))
call create_microlist(miniList, N_minilist, key_mask, microlist, idx_microlist, N_microlist, ptr_microlist, Nint)
do i=0,mo_tot_num*2
do k=ptr_microlist(i),ptr_microlist(i+1)-1
idx_microlist(k) = idx_minilist(idx_microlist(k))
end do
end do
do l=1,N_microlist(0)
do k=1,Nint
microlist_zero(k,1,l) = microlist(k,1,l)
microlist_zero(k,2,l) = microlist(k,2,l)
enddo
idx_microlist_zero(l) = idx_microlist(l)
enddo
end if
end if
do i_alpha=1,N_tq
! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
if(key_mask(1,1) /= 0) then
call getMobiles(tq(1,1,i_alpha), key_mask, mobiles, Nint)
if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
smallerlist = mobiles(1)
else
smallerlist = mobiles(2)
end if
do j=1,idx_alpha(0)
idx_alpha(j) = idx_miniList(idx_alpha(j))
end do
do l=0,N_microlist(smallerlist)-1
microlist_zero(:,:,ptr_microlist(1) + l) = microlist(:,:,ptr_microlist(smallerlist) + l)
idx_microlist_zero(ptr_microlist(1) + l) = idx_microlist(ptr_microlist(smallerlist) + l)
end do
call get_excitation_degree_vector(microlist_zero,tq(1,1,i_alpha),degree_alpha,Nint,N_microlist(smallerlist)+N_microlist(0),idx_alpha)
do j=1,idx_alpha(0)
idx_alpha(j) = idx_microlist_zero(idx_alpha(j))
end do
! i = 1
! j = 2
! do j = 2, idx_alpha_tmp(0)
! if(idx_alpha_tmp(j) < idx_alpha_tmp(j-1)) exit
! end do
!
! m = j
!
! idx_alpha(0) = idx_alpha_tmp(0)
!
! do l = 1, idx_alpha(0)
! if(j > idx_alpha_tmp(0)) then
! k = i
! i += 1
! else if(i >= m) then
! k = j
! j += 1
! else if(idx_alpha_tmp(i) < idx_alpha_tmp(j)) then
! k = i
! i += 1
! else
! k = j
! j += 1
! end if
! ! k=l
! idx_alpha(l) = idx_alpha_tmp(k)
! degree_alpha(l) = degree_alpha_tmp(k)
! end do
!
else
call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
do j=1,idx_alpha(0)
idx_alpha(j) = idx_miniList(idx_alpha(j))
end do
end if
! call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
! do j=1,idx_alpha(0)
! idx_alpha(j) = idx_miniList(idx_alpha(j))
! end do
!print *, idx_alpha(:idx_alpha(0))
do l_sd=1,idx_alpha(0)
k_sd = idx_alpha(l_sd)
@ -133,32 +235,17 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
tmp_det(k,1) = psi_ref(k,1,i_I)
tmp_det(k,2) = psi_ref(k,2,i_I)
enddo
! Hole (see list_to_bitstring)
iint = ishft(h1-1,-bit_kind_shift) + 1
ipos = h1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibclr(tmp_det(iint,s1),ipos)
! Particle
iint = ishft(p1-1,-bit_kind_shift) + 1
ipos = p1-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s1) = ibset(tmp_det(iint,s1),ipos)
if (degree_alpha(k_sd) == 2) then
! Hole (see list_to_bitstring)
iint = ishft(h2-1,-bit_kind_shift) + 1
ipos = h2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibclr(tmp_det(iint,s2),ipos)
! Particle
iint = ishft(p2-1,-bit_kind_shift) + 1
ipos = p2-ishft((iint-1),bit_kind_shift)-1
tmp_det(iint,s2) = ibset(tmp_det(iint,s2),ipos)
endif
logical :: ok
call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
if(.not. ok) cycle
! <I| \l/ |alpha>
do i_state=1,Nstates
dka(i_state) = 0.d0
enddo
do l_sd=k_sd+1,idx_alpha(0)
call get_excitation_degree(tmp_det,psi_non_ref(1,1,idx_alpha(l_sd)),degree,Nint)
if (degree == 0) then
@ -172,11 +259,12 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
if (.not.loop) then
call get_excitation(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),exc,degree,phase2,Nint)
hIl = hij_mrcc(idx_alpha(l_sd),i_I)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
! call i_h_j(psi_ref(1,1,i_I),psi_non_ref(1,1,idx_alpha(l_sd)),Nint,hIl)
do i_state=1,Nstates
dka(i_state) = hIl * lambda_mrcc(i_state,idx_alpha(l_sd)) * phase * phase2
enddo
endif
exit
endif
enddo
@ -211,21 +299,13 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
call omp_unset_lock( psi_ref_lock(i_I) )
enddo
enddo
deallocate (dIa_hla,hij_cache)
deallocate(miniList, idx_miniList)
!deallocate (dIa_hla,hij_cache)
!deallocate(miniList, idx_miniList)
end
BEGIN_PROVIDER [ integer(bit_kind), gen_det_sorted, (N_int,2,N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_shortcut, (0:N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_version, (N_int, N_det_generators,2) ]
&BEGIN_PROVIDER [ integer, gen_det_idx, (N_det_generators,2) ]
gen_det_sorted(:,:,:,1) = psi_det_generators(:,:,:N_det_generators)
gen_det_sorted(:,:,:,2) = psi_det_generators(:,:,:N_det_generators)
call sort_dets_ab_v(gen_det_sorted(:,:,:,1), gen_det_idx(:,1), gen_det_shortcut(0:,1), gen_det_version(:,:,1), N_det_generators, N_int)
call sort_dets_ba_v(gen_det_sorted(:,:,:,2), gen_det_idx(:,2), gen_det_shortcut(0:,2), gen_det_version(:,:,2), N_det_generators, N_int)
END_PROVIDER
subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq,miniList,N_miniList)
@ -258,6 +338,7 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
N_tq = 0
i_loop : do i=1,N_selected
if(is_connected_to(det_buffer(1,1,i), miniList, Nint, N_miniList)) then
cycle
@ -287,8 +368,84 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
end
subroutine find_triples_and_quadruples_micro(i_generator,n_selected,det_buffer,Nint,tq,N_tq,microlist,ptr_microlist,N_microlist,key_mask)
use bitmasks
implicit none
integer, intent(in) :: i_generator,n_selected, Nint
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k,m
logical :: is_in_wavefunction
integer :: degree(psi_det_size)
integer :: idx(0:psi_det_size)
logical :: good
integer(bit_kind), intent(out) :: tq(Nint,2,n_selected)
integer, intent(out) :: N_tq
integer :: nt,ni
logical, external :: is_connected_to
integer(bit_kind),intent(in) :: microlist(Nint,2,*)
integer,intent(in) :: ptr_microlist(0:*)
integer,intent(in) :: N_microlist(0:*)
integer(bit_kind),intent(in) :: key_mask(Nint, 2)
integer :: mobiles(2), smallerlist
N_tq = 0
i_loop : do i=1,N_selected
call getMobiles(det_buffer(1,1,i), key_mask, mobiles, Nint)
if(N_microlist(mobiles(1)) < N_microlist(mobiles(2))) then
smallerlist = mobiles(1)
else
smallerlist = mobiles(2)
end if
if(N_microlist(smallerlist) > 0) then
if(is_connected_to(det_buffer(1,1,i), microlist(1,1,ptr_microlist(smallerlist)), Nint, N_microlist(smallerlist))) then
cycle
end if
end if
if(N_microlist(0) > 0) then
if(is_connected_to(det_buffer(1,1,i), microlist, Nint, N_microlist(0))) then
cycle
end if
end if
! Select determinants that are triple or quadruple excitations
! from the ref
good = .True.
call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
!good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
do k=1,idx(0)
if (degree(k) < 3) then
good = .False.
exit
endif
enddo
if (good) then
if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
N_tq += 1
do k=1,N_int
tq(k,1,N_tq) = det_buffer(k,1,i)
tq(k,2,N_tq) = det_buffer(k,2,i)
enddo
endif
endif
enddo i_loop
end

View File

@ -11,12 +11,13 @@ subroutine mrcc_iterations
double precision :: E_new, E_old, delta_e
integer :: iteration,i_oscillations
double precision :: E_past(4)
double precision :: E_past(4), lambda
E_new = 0.d0
delta_E = 1.d0
iteration = 0
j = 1
i_oscillations = 0
lambda = 1.d0
do while (delta_E > 1.d-7)
iteration += 1
print *, '==========================='
@ -25,10 +26,15 @@ subroutine mrcc_iterations
print *, ''
E_old = sum(ci_energy_dressed)
call write_double(6,ci_energy_dressed(1),"MRCC energy")
call diagonalize_ci_dressed
call diagonalize_ci_dressed(lambda)
E_new = sum(ci_energy_dressed)
delta_E = dabs(E_new - E_old)
! if (E_new > E_old) then
! lambda = lambda * 0.7d0
! else
! lambda = min(1.d0, lambda * 1.1d0)
! endif
! print *, 'energy lambda ', lambda
E_past(j) = E_new
j +=1
call save_wavefunction

View File

@ -1,45 +1,54 @@
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
implicit none
BEGIN_DOC
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
END_DOC
integer :: i,k
double precision :: ihpsi(N_states),ihpsi_current(N_states)
integer :: i_pert_count
i_pert_count = 0
lambda_mrcc = 0.d0
do i=1,N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef, N_int, N_det_ref,size(psi_ref_coef,1), n_states, ihpsi_current)
do k=1,N_states
if (ihpsi_current(k) == 0.d0) then
ihpsi_current(k) = 1.d-32
endif
if(dabs(ihpsi_current(k) * psi_non_ref_coef(i,k)) < 1d-5) then
i_pert_count +=1
else
lambda_mrcc(k,i) = psi_non_ref_coef(i,k)/ihpsi_current(k)
endif
enddo
enddo
BEGIN_PROVIDER [ double precision, lambda_mrcc, (N_states,psi_det_size) ]
&BEGIN_PROVIDER [ integer, lambda_mrcc_pt2, (0:psi_det_size) ]
implicit none
BEGIN_DOC
! cm/<Psi_0|H|D_m> or perturbative 1/Delta_E(m)
END_DOC
integer :: i,k
double precision :: ihpsi_current(N_states)
integer :: i_pert_count
double precision :: hii, lambda_pert
integer :: N_lambda_mrcc_pt2
i_pert_count = 0
lambda_mrcc = 0.d0
N_lambda_mrcc_pt2 = 0
lambda_mrcc_pt2(0) = 0
do i=1,N_det_non_ref
call i_h_psi(psi_non_ref(1,1,i), psi_ref, psi_ref_coef_normalized, N_int, N_det_ref,&
size(psi_ref_coef,1), N_states,ihpsi_current)
call i_H_j(psi_non_ref(1,1,i),psi_non_ref(1,1,i),N_int,hii)
do k=1,N_states
if (ihpsi_current(k) == 0.d0) then
ihpsi_current(k) = 1.d-32
endif
lambda_mrcc(k,i) = min(0.d0,psi_non_ref_coef(i,k)/ihpsi_current(k) )
lambda_pert = 1.d0 / (psi_ref_energy_diagonalized(k)-hii)
if (lambda_pert / lambda_mrcc(k,i) < 0.5d0) then
i_pert_count += 1
lambda_mrcc(k,i) = 0.d0
if (lambda_mrcc_pt2(N_lambda_mrcc_pt2) /= i) then
N_lambda_mrcc_pt2 += 1
lambda_mrcc_pt2(N_lambda_mrcc_pt2) = i
endif
endif
enddo
enddo
lambda_mrcc_pt2(0) = N_lambda_mrcc_pt2
print*,'N_det_non_ref = ',N_det_non_ref
print*,'Number of ignored determinants = ',i_pert_count
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
print*,'lambda max = ',maxval(dabs(lambda_mrcc))
print*,'N_det_non_ref = ',N_det_non_ref
print*,'Number of ignored determinants = ',i_pert_count
print*,'psi_coef_ref_ratio = ',psi_ref_coef(2,1)/psi_ref_coef(1,1)
END_PROVIDER
!BEGIN_PROVIDER [ double precision, delta_ij_non_ref, (N_det_non_ref, N_det_non_ref,N_states) ]
!implicit none
!BEGIN_DOC
!! Dressing matrix in SD basis
!END_DOC
!delta_ij_non_ref = 0.d0
!call H_apply_mrcc_simple(delta_ij_non_ref,N_det_non_ref)
!END_PROVIDER
BEGIN_PROVIDER [ double precision, hij_mrcc, (N_det_non_ref,N_det_ref) ]
implicit none
@ -174,17 +183,19 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
END_PROVIDER
subroutine diagonalize_CI_dressed
subroutine diagonalize_CI_dressed(lambda)
implicit none
BEGIN_DOC
! Replace the coefficients of the CI states by the coefficients of the
! eigenstates of the CI matrix
END_DOC
double precision, intent(in) :: lambda
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
psi_coef(i,j) = lambda * CI_eigenvectors_dressed(i,j) + (1.d0 - lambda) * psi_coef(i,j)
enddo
call normalize(psi_coef(1,j), N_det)
enddo
SOFT_TOUCH psi_coef

View File

@ -3,7 +3,7 @@ import perturbation
END_SHELL
subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy)
implicit none
BEGIN_DOC
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@ -14,6 +14,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
integer(bit_kind),intent(in) :: key_mask(Nint,2)
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
double precision, intent(in) :: electronic_energy(N_st)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
@ -151,7 +152,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
idx_microlist_zero(ptr_microlist(1)+l) = idx_microlist(ptr_microlist(smallerlist)+l)
enddo
end if
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
c_pert,e_2_pert,H_pert_diag,Nint,N_microlist(smallerlist)+N_microlist(0), &
n_st,microlist_zero,idx_microlist_zero,N_microlist(smallerlist)+N_microlist(0))
else
@ -160,11 +161,11 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
cycle
end if
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
end if
! call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
! call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
! c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
do k = 1,N_st
@ -182,7 +183,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
end
subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp)
subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert,sum_norm_pert,sum_H_pert_diag,N_st,Nint,key_mask,fock_diag_tmp,electronic_energy)
implicit none
BEGIN_DOC
! Applly pertubration ``$PERT`` to the buffer of determinants generated in the H_apply
@ -193,6 +194,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
integer(bit_kind), intent(in) :: buffer(Nint,2,buffer_size)
integer(bit_kind),intent(in) :: key_mask(Nint,2)
double precision, intent(in) :: fock_diag_tmp(2,0:mo_tot_num)
double precision, intent(in) :: electronic_energy(N_st)
double precision, intent(inout) :: sum_norm_pert(N_st),sum_e_2_pert(N_st)
double precision, intent(inout) :: coef_pert_buffer(N_st,buffer_size),e_2_pert_buffer(N_st,buffer_size),sum_H_pert_diag(N_st)
double precision :: c_pert(N_st), e_2_pert(N_st), H_pert_diag(N_st)
@ -241,7 +243,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
cycle
endif
call pt2_$PERT(psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
call pt2_$PERT(electronic_energy,psi_det_generators(1,1,i_generator),buffer(1,1,i), fock_diag_tmp, &
c_pert,e_2_pert,H_pert_diag,Nint,N_minilist,n_st,minilist,idx_minilist,N_minilist)
do k = 1,N_st

View File

@ -29,11 +29,11 @@ subroutine pt2_epstein_nesbet ($arguments)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if(CI_electronic_energy(i)>h.and.CI_electronic_energy(i).ne.0.d0)then
if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_electronic_energy(i) - h)
else if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
@ -66,7 +66,6 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
double precision :: i_H_psi_array(N_st)
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
PROVIDE CI_electronic_energy
!call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
@ -74,7 +73,7 @@ subroutine pt2_epstein_nesbet_2x2 ($arguments)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if (i_H_psi_array(i) /= 0.d0) then
delta_e = h - CI_electronic_energy(i)
delta_e = h - electronic_energy(i)
if (delta_e > 0.d0) then
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
else
@ -165,7 +164,7 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
! : you repeat all the correlation energy already taken into account in electronic_energy(1)
!
! that could be repeated to this determinant.
!
@ -195,16 +194,16 @@ subroutine pt2_epstein_nesbet_SC2_projected ($arguments)
enddo
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
h = h + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
delta_E = 1.d0/(electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
@ -248,7 +247,7 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
!
! that can be repeated by repeating all the double excitations
!
! : you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
! : you repeat all the correlation energy already taken into account in electronic_energy(1)
!
! that could be repeated to this determinant.
!
@ -278,16 +277,16 @@ subroutine pt2_epstein_nesbet_SC2_no_projected ($arguments)
enddo
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
h = h + accu_e_corr
delta_E = 1.d0/(CI_SC2_electronic_energy(1) - h)
delta_E = 1.d0/(electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(CI_SC2_electronic_energy(1) - h)
c_pert(1) = i_H_psi_array(1) /(electronic_energy(1) - h)
e_2_pert(1) = i_H_psi_array(1) * c_pert(1)
do i =2,N_st
H_pert_diag(i) = h
if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(CI_SC2_electronic_energy(i) - h))
if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (-dabs(electronic_energy(i) - h))
e_2_pert(i) = (c_pert(i) * i_H_psi_array(i))
else
c_pert(i) = i_H_psi_array(i)
@ -328,11 +327,11 @@ subroutine pt2_epstein_nesbet_sc2 ($arguments)
h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
do i =1,N_st
if(CI_SC2_electronic_energy(i)>h.and.CI_SC2_electronic_energy(i).ne.0.d0)then
if(electronic_energy(i)>h.and.electronic_energy(i).ne.0.d0)then
c_pert(i) = -1.d0
e_2_pert(i) = selection_criterion*selection_criterion_factor*2.d0
else if (dabs(CI_SC2_electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (CI_SC2_electronic_energy(i) - h)
else if (dabs(electronic_energy(i) - h) > 1.d-6) then
c_pert(i) = i_H_psi_array(i) / (electronic_energy(i) - h)
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
else
@ -348,7 +347,7 @@ end
SUBST [ arguments, declarations ]
det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
electronic_energy,det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minilist,idx_minilist,N_minilist ;
integer, intent(in) :: Nint
integer, intent(in) :: ndet
@ -357,6 +356,7 @@ det_ref,det_pert,fock_diag_tmp,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st,minili
integer(bit_kind), intent(in) :: det_ref (Nint,2)
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(in) :: fock_diag_tmp(2,mo_tot_num+1)
double precision , intent(in) :: electronic_energy(N_st)
double precision , intent(out) :: c_pert(N_st)
double precision , intent(out) :: e_2_pert(N_st)
double precision, intent(out) :: H_pert_diag(N_st)

View File

@ -24,6 +24,21 @@ use bitmasks
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, psi_ref_coef_normalized, (psi_det_size,n_states) ]
implicit none
BEGIN_DOC
! Normalized coefficients of the reference
END_DOC
integer :: i,j,k
do k=1,N_states
do j=1,N_det_ref
psi_ref_coef_normalized(j,k) = psi_ref_coef(j,k)
enddo
call normalize(psi_ref_coef_normalized(1,k), N_det_ref)
enddo
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_ref_restart, (N_int,2,psi_det_size) ]

View File

@ -61,6 +61,7 @@ class H_apply(object):
s["params_post"] = ""
self.selection_pt2 = None
self.energy = "CI_electronic_energy"
self.perturbation = None
self.do_double_exc = do_double_exc
#s["omp_parallel"] = """!$OMP PARALLEL DEFAULT(NONE) &
@ -264,13 +265,13 @@ class H_apply(object):
self.data["keys_work"] = """
! if(check_double_excitation)then
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
"""%(pert)
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s)
"""%(pert,self.energy)
else:
self.data["keys_work"] = """
call perturb_buffer_by_mono_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp)
"""%(pert)
sum_norm_pert,sum_H_pert_diag,N_st,N_int,key_mask,fock_diag_tmp,%s)
"""%(pert,self.energy)
self.data["finalization"] = """

View File

@ -1728,3 +1728,55 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
deallocate (shortcut, sort_idx, sorted, version)
end
subroutine apply_excitation(det, exc, res, ok, Nint)
use bitmasks
implicit none
integer, intent(in) :: Nint
integer, intent(in) :: exc(0:2,2,2)
integer(bit_kind),intent(in) :: det(Nint, 2)
integer(bit_kind),intent(out) :: res(Nint, 2)
logical, intent(out) :: ok
integer :: h1,p1,h2,p2,s1,s2,degree
integer :: ii, pos
ok = .false.
degree = exc(0,1,1) + exc(0,1,2)
if(.not. (degree > 0 .and. degree <= 2)) then
print *, degree
print *, "apply ex"
STOP
endif
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
res = det
ii = (h1-1)/bit_kind_size + 1
pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s1) = ibclr(res(ii, s1), pos)
ii = (p1-1)/bit_kind_size + 1
pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s1) = ibset(res(ii, s1), pos)
if(degree == 2) then
ii = (h2-1)/bit_kind_size + 1
pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s2) = ibclr(res(ii, s2), pos)
ii = (p2-1)/bit_kind_size + 1
pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s2) = ibset(res(ii, s2), pos)
endif
ok = .true.
end subroutine

View File

@ -53,10 +53,10 @@ function test_exe() {
}
function run_HF() {
thresh=1.e-8
thresh=1.e-7
test_exe SCF || skip
ezfio set_file $1
ezfio set hartree_fock thresh_scf 1.e-10
ezfio set hartree_fock thresh_scf 1.e-11
qp_run SCF $1
energy="$(ezfio get hartree_fock energy)"
eq $energy $2 $thresh
@ -155,7 +155,7 @@ function run_all_1h_1p() {
ezfio set determinants read_wf True
qp_run mrcc_cassd $INPUT
energy="$(ezfio get mrcc_cassd energy)"
eq $energy -76.2289109271715 1.E-3
eq $energy -76.2284994316618 1.e-4
}
@ -166,7 +166,7 @@ function run_all_1h_1p() {
}
@test "SCF H2O VDZ pseudo" {
run_HF h2o_pseudo.ezfio -0.169483703904991E+02
run_HF h2o_pseudo.ezfio -16.9483708495521
}
@test "FCI H2O VDZ pseudo" {