10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-04 10:25:58 +02:00
quantum_package/plugins/shiftedbk/shifted_bk_iter.irp.f

160 lines
5.2 KiB
Fortran
Raw Normal View History

2018-05-23 18:16:33 +02:00
program shifted_bk
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:)
integer :: degree
integer :: n_det_before
double precision :: threshold_davidson_in
allocate (pt2(N_states))
double precision :: hf_energy_ref
logical :: has
double precision :: relative_error, absolute_error
integer :: N_states_p
character*(512) :: fmt
PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
PROVIDE psi_bilinear_matrix_rows psi_det_sorted_order psi_bilinear_matrix_order
PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
PROVIDE psi_bilinear_matrix_transp_order
pt2 = -huge(1.e0)
threshold_davidson_in = threshold_davidson
threshold_davidson = threshold_davidson_in * 100.d0
SOFT_TOUCH threshold_davidson
call diagonalize_CI_dressed
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
endif
if (N_det > N_det_max) then
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
call diagonalize_CI_dressed
call save_wavefunction
N_states_p = min(N_det,N_states)
endif
n_det_before = 0
character*(8) :: pt2_string
double precision :: threshold_selectors_save, threshold_generators_save
threshold_selectors_save = threshold_selectors
threshold_generators_save = threshold_generators
double precision :: error(N_states), energy(N_states)
error = 0.d0
threshold_selectors = 1.d0
threshold_generators = 1d0
if (.True.) then
pt2_string = '(sh-Bk) '
do while ( (N_det < N_det_max) )
write(*,'(A)') '--------------------------------------------------------------------------------'
N_det_delta_ij = N_det
do i=1,N_states
energy(i) = psi_energy(i)+nuclear_repulsion
enddo
PROVIDE delta_ij_tmp
call delta_ij_done()
call diagonalize_ci_dressed
do i=1,N_states
pt2(i) = ci_energy_dressed(i) - energy(i)
enddo
N_states_p = min(N_det,N_states)
print *, ''
print '(A,I12)', 'Summary at N_det = ', N_det
print '(A)', '-----------------------------------'
print *, ''
print *, ''
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(12X,', N_states_p, '(6X,A7,1X,I6,10X))'
write(*,fmt) ('State',k, k=1,N_states_p)
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
write(fmt,*) '(A12,', N_states_p, '(1X,F14.8,15X))'
write(*,fmt) '# E ', energy(1:N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', energy(1:N_states_p)-energy(1)
write(*,fmt) '# Excit. (eV)', (energy(1:N_states_p)-energy(1))*27.211396641308d0
endif
write(fmt,*) '(A12,', 2*N_states_p, '(1X,F14.8))'
write(*,fmt) '# PT2'//pt2_string, (pt2(k), error(k), k=1,N_states_p)
write(*,'(A)') '#'
write(*,fmt) '# E+PT2 ', (energy(k)+pt2(k),error(k), k=1,N_states_p)
if (N_states_p > 1) then
write(*,fmt) '# Excit. (au)', ( (energy(k)+pt2(k)-energy(1)-pt2(1)), &
dsqrt(error(k)*error(k)+error(1)*error(1)), k=1,N_states_p)
write(*,fmt) '# Excit. (eV)', ( (energy(k)+pt2(k)-energy(1)-pt2(1))*27.211396641308d0, &
dsqrt(error(k)*error(k)+error(1)*error(1))*27.211396641308d0, k=1,N_states_p)
endif
write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
write(*,fmt)
print *, ''
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
do k=1, N_states_p
print*,'State ',k
print *, 'PT2 = ', pt2(k)
print *, 'E = ', energy(k)
print *, 'E+PT2'//pt2_string//' = ', energy(k)+pt2(k)
enddo
print *, '-----'
if(N_states.gt.1)then
print *, 'Variational Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (energy(i) - energy(1)), &
(energy(i) - energy(1)) * 27.211396641308d0
enddo
print *, '-----'
print*, 'Variational + perturbative Energy difference (au | eV)'
do i=2, N_states_p
print*,'Delta E = ', (energy(i)+ pt2(i) - (energy(1) + pt2(1))), &
(energy(i)+ pt2(i) - (energy(1) + pt2(1))) * 27.211396641308d0
enddo
endif
call ezfio_set_shiftedbk_energy_pt2(energy(1)+pt2(1))
! call dump_fci_iterations_value(N_det,energy,pt2)
n_det_before = N_det
PROVIDE psi_coef
PROVIDE psi_det
PROVIDE psi_det_sorted
if (N_det >= N_det_max) then
threshold_davidson = threshold_davidson_in
end if
call save_wavefunction
call ezfio_set_shiftedbk_energy(energy(1))
call ezfio_set_shiftedbk_energy_pt2(ci_energy_dressed(1))
enddo
endif
end