From fb46e55f18cd23c156d1ad9d4dcf810ad8d2c3b3 Mon Sep 17 00:00:00 2001 From: Emmanuel Giner Date: Thu, 17 Jan 2019 10:57:27 +0100 Subject: [PATCH] modify prints in print_summary and added expectation value of s2 independant from eigenvectors --- docs/source/modules/fci.rst | 14 ++++ docs/source/modules/hartree_fock.rst | 14 ---- scripts/qp_e_conv_fci | 4 +- src/cipsi/cipsi.irp.f | 4 +- src/cipsi/slave_cipsi.irp.f | 4 +- src/cipsi/stochastic_cipsi.irp.f | 4 +- src/davidson/u0_h_u0.irp.f | 106 +++++++++++++++------------ src/dressing/energy.irp.f | 4 - src/fci/pt2.irp.f | 2 +- src/iterations/print_summary.irp.f | 7 +- src/kohn_sham_rs/61.rsks.bats | 1 + 11 files changed, 86 insertions(+), 78 deletions(-) diff --git a/docs/source/modules/fci.rst b/docs/source/modules/fci.rst index 9f08fb65..fc4a9f20 100644 --- a/docs/source/modules/fci.rst +++ b/docs/source/modules/fci.rst @@ -99,6 +99,20 @@ Subroutines / functions +.. c:function:: run + + .. code:: text + + subroutine run + + File: :file:`pt2.irp.f` + + + + + + + .. c:function:: save_energy .. code:: text diff --git a/docs/source/modules/hartree_fock.rst b/docs/source/modules/hartree_fock.rst index e115ba57..c81b9aca 100644 --- a/docs/source/modules/hartree_fock.rst +++ b/docs/source/modules/hartree_fock.rst @@ -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 - - diff --git a/scripts/qp_e_conv_fci b/scripts/qp_e_conv_fci index 78f05777..beaea930 100755 --- a/scripts/qp_e_conv_fci +++ b/scripts/qp_e_conv_fci @@ -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 diff --git a/src/cipsi/cipsi.irp.f b/src/cipsi/cipsi.irp.f index cd13fedf..5ee8f24f 100644 --- a/src/cipsi/cipsi.irp.f +++ b/src/cipsi/cipsi.irp.f @@ -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() diff --git a/src/cipsi/slave_cipsi.irp.f b/src/cipsi/slave_cipsi.irp.f index 4236e13e..755717b8 100644 --- a/src/cipsi/slave_cipsi.irp.f +++ b/src/cipsi/slave_cipsi.irp.f @@ -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 diff --git a/src/cipsi/stochastic_cipsi.irp.f b/src/cipsi/stochastic_cipsi.irp.f index 5871d0e2..f92c68fd 100644 --- a/src/cipsi/stochastic_cipsi.irp.f +++ b/src/cipsi/stochastic_cipsi.irp.f @@ -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() diff --git a/src/davidson/u0_h_u0.irp.f b/src/davidson/u0_h_u0.irp.f index faea2558..0dd5337f 100644 --- a/src/davidson/u0_h_u0.irp.f +++ b/src/davidson/u0_h_u0.irp.f @@ -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 - diff --git a/src/dressing/energy.irp.f b/src/dressing/energy.irp.f index 7e223a2e..66db2fd5 100644 --- a/src/dressing/energy.irp.f +++ b/src/dressing/energy.irp.f @@ -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) diff --git a/src/fci/pt2.irp.f b/src/fci/pt2.irp.f index b196a860..8fe013aa 100644 --- a/src/fci/pt2.irp.f +++ b/src/fci/pt2.irp.f @@ -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 diff --git a/src/iterations/print_summary.irp.f b/src/iterations/print_summary.irp.f index dd050d18..a8037982 100644 --- a/src/iterations/print_summary.irp.f +++ b/src/iterations/print_summary.irp.f @@ -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 *, '' diff --git a/src/kohn_sham_rs/61.rsks.bats b/src/kohn_sham_rs/61.rsks.bats index 2705181f..7fcdecb7 100644 --- a/src/kohn_sham_rs/61.rsks.bats +++ b/src/kohn_sham_rs/61.rsks.bats @@ -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