9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-01 10:15:18 +02:00
qp2/src/cipsi/stochastic_cipsi.irp.f
2020-09-01 11:43:27 -05:00

167 lines
4.6 KiB
Fortran

subroutine run_stochastic_cipsi
implicit none
BEGIN_DOC
! Selected Full Configuration Interaction with Stochastic selection and PT2.
END_DOC
integer :: i,j,k
double precision, allocatable :: pt2(:), variance(:), norm2(:), rpt2(:), zeros(:)
integer :: to_select
logical, external :: qp_stop
double precision :: rss
double precision, external :: memory_of_double
PROVIDE H_apply_buffer_allocated
N_iter = 1
threshold_generators = 1.d0
SOFT_TOUCH threshold_generators
rss = memory_of_double(N_states)*4.d0
call check_mem(rss,irp_here)
allocate (pt2(N_states), zeros(N_states), rpt2(N_states), norm2(N_states), variance(N_states))
double precision :: hf_energy_ref
logical :: has
double precision :: relative_error
relative_error=PT2_relative_error
zeros = 0.d0
pt2 = -huge(1.e0)
rpt2 = -huge(1.e0)
norm2 = 0.d0
variance = huge(1.e0)
if (s2_eig) then
call make_s2_eigenfunction
endif
if (is_complex) then
call diagonalize_ci_complex
else
call diagonalize_ci
endif
call save_wavefunction
call ezfio_has_hartree_fock_energy(has)
if (has) then
call ezfio_get_hartree_fock_energy(hf_energy_ref)
else
hf_energy_ref = ref_bitmask_energy_with_nucl_rep
endif
if (N_det > N_det_max) then
psi_det = psi_det_sorted
if (is_complex) then
psi_coef_complex = psi_coef_sorted_complex
N_det = N_det_max
soft_touch N_det psi_det psi_coef_complex
else
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
endif
if (s2_eig) then
call make_s2_eigenfunction
endif
if (is_complex) then
call diagonalize_ci_complex
else
call diagonalize_CI
endif
call save_wavefunction
endif
double precision :: correlation_energy_ratio
double precision :: error(N_states)
correlation_energy_ratio = 0.d0
do while ( &
(N_det < N_det_max) .and. &
(maxval(abs(rpt2(1:N_states))) > pt2_max) .and. &
(maxval(abs(variance(1:N_states))) > variance_max) .and. &
(correlation_energy_ratio <= correlation_energy_ratio_max) &
)
write(*,'(A)') '--------------------------------------------------------------------------------'
to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
to_select = max(N_states_diag, to_select)
pt2 = 0.d0
variance = 0.d0
norm2 = 0.d0
call ZMQ_pt2(psi_energy_with_nucl_rep,pt2,relative_error,error, variance, &
norm2, to_select) ! Stochastic PT2 and selection
do k=1,N_states
rpt2(k) = pt2(k)/(1.d0 + norm2(k))
enddo
correlation_energy_ratio = (psi_energy_with_nucl_rep(1) - hf_energy_ref) / &
(psi_energy_with_nucl_rep(1) + rpt2(1) - hf_energy_ref)
correlation_energy_ratio = min(1.d0,correlation_energy_ratio)
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2)
call save_energy(psi_energy_with_nucl_rep, rpt2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
call print_extrapolated_energy()
N_iter += 1
if (qp_stop()) exit
! Add selected determinants
call copy_h_apply_buffer_to_wf()
! call save_wavefunction
if (is_complex) then
PROVIDE psi_coef_complex
else
PROVIDE psi_coef
endif
PROVIDE psi_det
PROVIDE psi_det_sorted
if (is_complex) then
call diagonalize_ci_complex
else
call diagonalize_CI
endif
call save_wavefunction
call save_energy(psi_energy_with_nucl_rep, zeros)
if (qp_stop()) exit
enddo
if (.not.qp_stop()) then
if (N_det < N_det_max) then
if (is_complex) then
call diagonalize_ci_complex
else
call diagonalize_CI
endif
call save_wavefunction
call save_energy(psi_energy_with_nucl_rep, zeros)
endif
pt2(:) = 0.d0
variance(:) = 0.d0
norm2(:) = 0.d0
call ZMQ_pt2(psi_energy_with_nucl_rep, pt2,relative_error,error,variance, &
norm2,0) ! Stochastic PT2
do k=1,N_states
rpt2(k) = pt2(k)/(1.d0 + norm2(k))
enddo
call save_energy(psi_energy_with_nucl_rep, rpt2)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm2,N_det,N_occ_pattern,N_states,psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
call print_extrapolated_energy()
endif
end