modify prints in print_summary and added expectation value of s2 independant from eigenvectors

This commit is contained in:
Emmanuel Giner 2019-01-17 10:57:27 +01:00
parent c461a7840d
commit fb46e55f18
11 changed files with 86 additions and 78 deletions

View File

@ -99,6 +99,20 @@ Subroutines / functions
.. c:function:: run
.. code:: text
subroutine run
File: :file:`pt2.irp.f`
.. c:function:: save_energy
.. code:: text

View File

@ -202,17 +202,3 @@ Subroutines / functions
Create a MO guess if no MOs are present in the EZFIO directory
.. c:function:: run
.. code:: text
subroutine run
File: :file:`scf.irp.f`
Run SCF calculation

View File

@ -26,7 +26,7 @@ plot "$out" w lp title "E_{var} state $i", "$out" u 1:3 w lp title "E_{var} + PT
EOF
gnuplot ${out}.plt
rm ${out}.plt
#rm ${out}.plt
done
@ -43,5 +43,5 @@ plot "$out" w lp title "Delta E_{var} state $i", "$out" u 1:3 w lp title "Delta
EOF
gnuplot ${out}.plt
rm ${out}.plt
# rm ${out}.plt
done

View File

@ -89,7 +89,7 @@ subroutine run_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2)
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,psi_s2)
do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k))
@ -143,7 +143,7 @@ subroutine run_cipsi
rpt2(:) = pt2(:)/(1.d0 + norm(k))
enddo
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,psi_s2)
call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)
call print_extrapolated_energy()

View File

@ -118,7 +118,7 @@ subroutine run_slave_main
IRP_ENDIF
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
psi_energy(1:N_states) = energy(1:N_states)
TOUCH psi_energy state_average_weight threshold_generators
TOUCH psi_energy psi_s2 state_average_weight threshold_generators
if (mpi_master) then
print *, 'N_det', N_det
@ -222,7 +222,7 @@ subroutine run_slave_main
IRP_ENDIF
if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
psi_energy(1:N_states) = energy(1:N_states)
TOUCH psi_energy state_average_weight pt2_stoch_istate threshold_generators
TOUCH psi_energy psi_s2 state_average_weight pt2_stoch_istate threshold_generators
if (mpi_master) then
print *, 'N_det', N_det
print *, 'N_det_generators', N_det_generators

View File

@ -83,7 +83,7 @@ subroutine run_stochastic_cipsi
call save_energy(psi_energy_with_nucl_rep, pt2)
call write_double(6,correlation_energy_ratio, 'Correlation ratio')
call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states)
call print_summary(psi_energy_with_nucl_rep,pt2,error,variance,norm,N_det,N_occ_pattern,N_states,psi_s2)
do k=1,N_states
rpt2(:) = pt2(:)/(1.d0 + norm(k))
@ -125,7 +125,7 @@ subroutine run_stochastic_cipsi
rpt2(:) = pt2(:)/(1.d0 + norm(k))
enddo
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,N_states)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,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()

View File

@ -1,12 +1,16 @@
BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
&BEGIN_PROVIDER [ double precision, psi_s2, (N_states) ]
implicit none
BEGIN_DOC
! Electronic energy of the current wave function
! psi_energy(i) = $\langle \Psi_i | H | \Psi_i \rangle$
!
! psi_s2(i) = $\langle \Psi_i | S^2 | \Psi_i \rangle$
END_DOC
call u_0_H_u_0(psi_energy,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
call u_0_H_u_0(psi_energy,psi_s2,psi_coef,N_det,psi_det,N_int,N_states,psi_det_size)
integer :: i
do i=N_det+1,N_states
psi_energy(i) = 0.d0
psi_s2(i) = 0.d0
enddo
END_PROVIDER
@ -19,6 +23,57 @@ BEGIN_PROVIDER [ double precision, psi_energy_with_nucl_rep, (N_states) ]
END_PROVIDER
subroutine u_0_H_u_0(e_0,s_0,u_0,n,keys_tmp,Nint,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes $E_0 = \frac{\langle u_0|H|u_0 \rangle}{\langle u_0|u_0 \rangle}$
!
! and $S_0 = \frac{\langle u_0|S^2|u_0 \rangle}{\langle u_0|u_0 \rangle}$
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint, N_st, sze
double precision, intent(out) :: e_0(N_st),s_0(N_st)
double precision, intent(inout) :: u_0(sze,N_st)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision, allocatable :: v_0(:,:), s_vec(:,:), u_1(:,:)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
if ((n > 100000).and.distributed_davidson) then
allocate (v_0(n,N_states_diag),s_vec(n,N_states_diag), u_1(n,N_states_diag))
u_1(1:n,1:N_states) = u_0(1:n,1:N_states)
u_1(1:n,N_states+1:N_states_diag) = 0.d0
call H_S2_u_0_nstates_zmq(v_0,s_vec,u_1,N_states_diag,n)
deallocate(u_1)
else
allocate (v_0(n,N_st),s_vec(n,N_st),u_1(n,N_st))
u_1(1:n,:) = u_0(1:n,:)
call H_S2_u_0_nstates_openmp(v_0,s_vec,u_1,N_st,n)
u_0(1:n,:) = u_1(1:n,:)
deallocate(u_1)
endif
double precision :: norm
!$OMP PARALLEL DO PRIVATE(i,norm) DEFAULT(SHARED)
do i=1,N_st
norm = u_dot_u(u_0(1,i),n)
if (norm /= 0.d0) then
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)
s_0(i) = u_dot_v(s_vec(1,i),u_0(1,i),n)
else
e_0(i) = 0.d0
s_0(i) = 0.d0
endif
enddo
!$OMP END PARALLEL DO
deallocate (s_vec, v_0)
end
subroutine H_S2_u_0_nstates_openmp(v_0,s_0,u_0,N_st,sze)
@ -576,48 +631,3 @@ N_int;;
END_TEMPLATE
subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze)
use bitmasks
implicit none
BEGIN_DOC
! Computes $E_0 = \frac{\langle u_0|H|u_0 \rangle}{\langle u_0|u_0 \rangle}$
!
! n : number of determinants
!
END_DOC
integer, intent(in) :: n,Nint, N_st, sze
double precision, intent(out) :: e_0(N_st)
double precision, intent(inout) :: u_0(sze,N_st)
integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n)
double precision, allocatable :: v_0(:,:), s_0(:,:), u_1(:,:)
double precision :: u_dot_u,u_dot_v,diag_H_mat_elem
integer :: i,j
if ((n > 100000).and.distributed_davidson) then
allocate (v_0(n,N_states_diag),s_0(n,N_states_diag), u_1(n,N_states_diag))
u_1(1:n,1:N_states) = u_0(1:n,1:N_states)
u_1(1:n,N_states+1:N_states_diag) = 0.d0
call H_S2_u_0_nstates_zmq(v_0,s_0,u_1,N_states_diag,n)
deallocate(u_1)
else
allocate (v_0(n,N_st),s_0(n,N_st),u_1(n,N_st))
u_1(1:n,:) = u_0(1:n,:)
call H_S2_u_0_nstates_openmp(v_0,s_0,u_1,N_st,n)
u_0(1:n,:) = u_1(1:n,:)
deallocate(u_1)
endif
double precision :: norm
!$OMP PARALLEL DO PRIVATE(i,norm) DEFAULT(SHARED)
do i=1,N_st
norm = u_dot_u(u_0(1,i),n)
if (norm /= 0.d0) then
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)
else
e_0(i) = 0.d0
endif
enddo
!$OMP END PARALLEL DO
deallocate (s_0, v_0)
end

View File

@ -22,10 +22,6 @@ BEGIN_PROVIDER [ double precision, dress_E0_denominator, (N_states) ]
print *, h0_type, ' not implemented'
stop
endif
! call u_0_H_u_0(dress_E0_denominator,psi_coef,N_det,psi_det,N_int,N_states,size(psi_coef,1))
! do i=N_det+1,N_states
! dress_E0_denominator(i) = 0.d0
! enddo
call write_double(6,dress_E0_denominator(1)+nuclear_repulsion, 'dress Energy denominator')
else
dress_E0_denominator = -huge(1.d0)

View File

@ -44,7 +44,7 @@ subroutine run
rpt2(:) = pt2(:)/(1.d0 + norm(k))
enddo
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern)
call print_summary(psi_energy_with_nucl_rep(1:N_states),pt2,error,variance,norm,N_det,N_occ_pattern,psi_s2)
call save_energy(E_CI_before,pt2)
end

View File

@ -1,11 +1,11 @@
subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_st)
subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_st,s2_)
implicit none
BEGIN_DOC
! Print the extrapolated energy in the output
END_DOC
integer, intent(in) :: n_det_, n_occ_pattern_, n_st
double precision, intent(in) :: e_(n_st), pt2_(n_st), variance_(n_st), norm_(n_st), error_(n_st)
double precision, intent(in) :: e_(n_st), pt2_(n_st), variance_(n_st), norm_(n_st), error_(n_st), s2_(n_st)
integer :: i, k
integer :: N_states_p
character*(9) :: pt2_string
@ -65,11 +65,12 @@ subroutine print_summary(e_,pt2_,error_,variance_,norm_,n_det_,n_occ_pattern_,n_
do k=1, N_states_p
print*,'* State ',k
print *, '< S^2 > = ', s2_(k)
print *, 'E = ', e_(k)
print *, 'Variance = ', variance_(k)
print *, 'PT norm = ', dsqrt(norm_(k))
print *, 'PT2 = ', pt2_(k)
print *, 'rPT2 = ', pt2_(k)*f(k)
print *, 'E = ', e_(k)
print *, 'E+PT2 '//pt2_string//' = ', e_(k)+pt2_(k), ' +/- ', error_(k)
print *, 'E+rPT2'//pt2_string//' = ', e_(k)+pt2_(k)*f(k), ' +/- ', error_(k)*f(k)
print *, ''

View File

@ -14,6 +14,7 @@ function run() {
qp set dft_keywords correlation_functional $functional
qp set ao_two_e_erf_ints mu_erf 0.5
qp set becke_numerical_grid grid_type_sgn 1
qp_reset --mos $1
qp run rs_ks_scf
energy="$(ezfio get kohn_sham_rs energy)"
eq $energy $3 $thresh