10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-09-27 03:51:01 +02:00

added qp_e_conv_fci and print_e_conv.irp.f

This commit is contained in:
Emmanuel Giner 2019-01-04 02:17:55 +01:00
parent a45a5a9690
commit d924386eed
15 changed files with 135 additions and 216 deletions

View File

@ -12,8 +12,8 @@
:caption: Introduction
:hidden:
intro/selected_ci
intro/install
intro/selected_ci
.. toctree::
:maxdepth: 1

View File

@ -11,8 +11,18 @@ The |qp|
What it is
==========
The |qp| is an open-source programming environment for quantum chemistry,
especially for `wave function theory <https://en.wikipedia.org/wiki/Ab_initio_quantum_chemistry_methods>`_ (|WFT|).
The |qp| is an open-source **programming environment** for quantum chemistry.
It has been built from the **developper** point of view in order to help
the design of new quantum chemistry methods,
especially for `wave function theory <https://en.wikipedia.org/wiki/Ab_initio_quantum_chemistry_methods>`_ (|WFT|).
From the **user** point of view, the |qp| proposes a stand-alone path
to use optimized selected configuration interaction |sCI| based on the
|CIPSI| algorithm that can efficiently reach near-full configuration interaction
|FCI| quality for relatively large systems (see for instance :cite:`Caffarel_2016,Caffarel_2016.2,Loos_2018,Scemama_2018,Dash_2018,Garniron_2017.2,Loos_2018,Garniron_2018,Giner2018Oct`).
To have a simple example of how to use the |CIPSI| program, go to the `users_guide/quickstart.
The main goal is the development of selected configuration interaction |sCI|
methods and multi-reference perturbation theory |MRPT| in the
determinant-driven paradigm. It also contains the very basics of Kohn-Sham `density functional theory <https://en.wikipedia.org/wiki/Density_functional_theory>`_ |KS-DFT| and `range-separated hybrids <https://aip.scitation.org/doi/10.1063/1.1383587>`_ |RSH|.

View File

@ -786,7 +786,7 @@ Subroutines / functions
File: :file:`pot_ao_erf_ints.irp.f`
subroutine that returs all integrals over r of type erf(mu_in * | r-C_center |)/| r-C_center |
subroutine that returs all integrals over r of type erf(mu_in * | r-C_center | )/| r-C_center |
@ -890,7 +890,7 @@ Subroutines / functions
File: :file:`pot_ao_erf_ints.irp.f`
computes the following integral : int[-infty;+infty] dr AO_i_ao (r) AO_j_ao(r) erf(mu_in * | r-C_center |)/| r-C_center |
computes the following integral : int[-infty;+infty] dr AO_i_ao (r) AO_j_ao(r) erf(mu_in * | r-C_center | )/| r-C_center |

View File

@ -146,9 +146,6 @@ Providers
File: :file:`grid_becke_vector.irp.f`
final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
<<<<<<< HEAD
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
=======
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
@ -187,7 +184,6 @@ Providers
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
@ -195,48 +191,6 @@ Providers
.. c:var:: final_weight_at_r
.. code:: text
double precision, allocatable :: final_weight_at_r (n_points_integration_angular,n_points_radial_grid,nucl_num)
File: :file:`grid_becke.irp.f`
<<<<<<< HEAD
Total weight on each grid point which takes into account all Lebedev, Voronoi and radial weights.
=======
final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
final_weight_functions_at_final_grid_points(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
.. c:var:: final_weight_at_r_vector
.. code:: text
double precision, allocatable :: final_grid_points (3,n_points_final_grid)
double precision, allocatable :: final_weight_at_r_vector (n_points_final_grid)
integer, allocatable :: index_final_points (3,n_points_final_grid)
integer, allocatable :: index_final_points_reverse (n_points_integration_angular,n_points_radial_grid,nucl_num)
File: :file:`grid_becke_vector.irp.f`
final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
.. c:var:: grid_points_per_atom
.. code:: text
@ -276,13 +230,9 @@ Providers
File: :file:`grid_becke_vector.irp.f`
final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
<<<<<<< HEAD
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
=======
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
@ -302,13 +252,9 @@ Providers
File: :file:`grid_becke_vector.irp.f`
final_grid_points(1:3,j) = (/ x, y, z /) of the jth grid point
<<<<<<< HEAD
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
=======
final_weight_at_r_vector(i) = Total weight function of the ith grid point which contains the Lebedev, Voronoi and radial weights contributions
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
index_final_points(1:3,i) = gives the angular, radial and atomic indices associated to the ith grid point
index_final_points_reverse(i,j,k) = index of the grid point having i as angular, j as radial and l as atomic indices
@ -365,13 +311,9 @@ Providers
File: :file:`grid_becke.irp.f`
n_points_radial_grid = number of radial grid points per atom
<<<<<<< HEAD
n_points_integration_angular = number of angular grid points per atom
=======
n_points_integration_angular = number of angular grid points per atom
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
These numbers are automatically set by setting the grid_type_sgn parameter
@ -387,37 +329,15 @@ Providers
File: :file:`grid_becke.irp.f`
n_points_radial_grid = number of radial grid points per atom
<<<<<<< HEAD
n_points_integration_angular = number of angular grid points per atom
=======
n_points_integration_angular = number of angular grid points per atom
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
These numbers are automatically set by setting the grid_type_sgn parameter
.. c:var:: weight_at_r
<<<<<<< HEAD
=======
.. code:: text
double precision, allocatable :: weight_at_r (n_points_integration_angular,n_points_radial_grid,nucl_num)
File: :file:`grid_becke.irp.f`
Weight function at grid points : w_n(r) according to the equation (22) of Becke original paper (JCP, 88, 1988)
The "n" discrete variable represents the nucleis which in this array is represented by the last dimension and the points are labelled by the other dimensions.
.. c:var:: weight_functions_at_grid_points
>>>>>>> 36234f0822384b0bce01530d210522e7975b759b
.. code:: text

View File

@ -100,7 +100,9 @@ Providers
File: :file:`hf_energy.irp.f`
Extra contribution to the SCF energy coming from the density.
For a Hartree-Fock calculation: extra_e_contrib_density = 0
For a Kohn-Sham or Range-separated Kohn-Sham: the exchange/correlation - trace of the V_xc potential

View File

@ -106,23 +106,6 @@ Providers
.. c:var:: extra_e_contrib_density
.. code:: text
double precision :: extra_e_contrib_density
File: :file:`ks_enery.irp.f`
Extra contribution to the SCF energy coming from the density.
For a Hartree-Fock calculation: extra_e_contrib_density = 0
For a Kohn-Sham or Range-separated Kohn-Sham: the exchange/correlation - 1/2 trace of the V_xc potential
.. c:var:: fock_matrix_alpha_no_xc_ao
.. code:: text

View File

@ -52,7 +52,7 @@ EZFIO parameters
.. option:: energy
SCF energy
Energy range separated hybrid

View File

@ -203,20 +203,6 @@ Subroutines / functions
.. c:function:: perturb_buffer_by_mono_decontracted
.. code:: text
subroutine perturb_buffer_by_mono_decontracted(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)
File: :file:`perturbation.irp.f_shell_13`
Applly pertubration ``decontracted`` to the buffer of determinants generated in the H_apply routine.
.. c:function:: perturb_buffer_by_mono_dummy
.. code:: text
@ -329,20 +315,6 @@ Subroutines / functions
.. c:function:: perturb_buffer_decontracted
.. code:: text
subroutine perturb_buffer_decontracted(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)
File: :file:`perturbation.irp.f_shell_13`
Applly pertubration ``decontracted`` to the buffer of determinants generated in the H_apply routine.
.. c:function:: perturb_buffer_dummy
.. code:: text
@ -455,27 +427,13 @@ Subroutines / functions
.. c:function:: pt2_decontracted
.. code:: text
subroutine pt2_decontracted (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)
File: :file:`pt2_equations.irp.f_template_412`
.. c:function:: pt2_dummy
.. code:: text
subroutine pt2_dummy (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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
Dummy perturbation to add all connected determinants.
@ -489,7 +447,7 @@ Subroutines / functions
subroutine pt2_epstein_nesbet (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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
@ -511,7 +469,7 @@ Subroutines / functions
subroutine pt2_epstein_nesbet_2x2 (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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
@ -533,7 +491,7 @@ Subroutines / functions
subroutine pt2_epstein_nesbet_2x2_no_ci_diag(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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
@ -577,7 +535,7 @@ Subroutines / functions
subroutine pt2_moller_plesset (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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
@ -599,7 +557,7 @@ Subroutines / functions
subroutine pt2_moller_plesset_general (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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
compute the general Moller-Plesset perturbative first order coefficient and second order energetic contribution
@ -621,7 +579,7 @@ Subroutines / functions
subroutine pt2_qdpt (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)
File: :file:`pt2_equations.irp.f_template_412`
File: :file:`pt2_equations.irp.f_template_365`
compute the QDPT first order coefficient and second order energetic contribution

View File

@ -404,9 +404,6 @@ Index of Providers
* :c:data:`n_det_max_full`
* :c:data:`n_det_non_cas`
* :c:data:`n_det_print_wf`
* :c:data:`n_det_non_ref`
* :c:data:`n_det_print_wf`
* :c:data:`n_det_ref`
* :c:data:`n_det_selectors`
* :c:data:`n_double_exc_bitmasks`
* :c:data:`n_double_selectors`
@ -1134,7 +1131,6 @@ Index of Subroutines/Functions
* :c:func:`overlap_x_abs`
* :c:func:`past_d1`
* :c:func:`past_d2`
* :c:func:`perturb_buffer_by_mono_decontracted`
* :c:func:`perturb_buffer_by_mono_dummy`
* :c:func:`perturb_buffer_by_mono_epstein_nesbet`
* :c:func:`perturb_buffer_by_mono_epstein_nesbet_2x2`
@ -1143,7 +1139,6 @@ Index of Subroutines/Functions
* :c:func:`perturb_buffer_by_mono_moller_plesset`
* :c:func:`perturb_buffer_by_mono_moller_plesset_general`
* :c:func:`perturb_buffer_by_mono_qdpt`
* :c:func:`perturb_buffer_decontracted`
* :c:func:`perturb_buffer_dummy`
* :c:func:`perturb_buffer_epstein_nesbet`
* :c:func:`perturb_buffer_epstein_nesbet_2x2`
@ -1168,7 +1163,6 @@ Index of Subroutines/Functions
* :c:func:`provide_everything`
* :c:func:`pt2`
* :c:func:`pt2_collector`
* :c:func:`pt2_decontracted`
* :c:func:`pt2_dummy`
* :c:func:`pt2_epstein_nesbet`
* :c:func:`pt2_epstein_nesbet_2x2`
@ -1368,4 +1362,4 @@ Index of Subroutines/Functions
* :c:func:`zmq_put_psi_det_beta_unique`
* :c:func:`zmq_put_psi_det_size`
* :c:func:`zmq_selection`
* :c:func:`zmq_set_running`
* :c:func:`zmq_set_running`

View File

@ -2,8 +2,9 @@
Quick-start guide
=================
This tutorial should talk you through everything you need to get started with
the |qp|. As an example, we will run a |CIPSI| calculation on the HCN molecule.
This tutorial should teach you everything you need to get started with
the the basics of the |qp|.
As an example, we will run a |CIPSI| calculation on the HCN molecule in the 631-G basis set.
Demo video

View File

@ -1,58 +1,17 @@
#!/bin/bash
file=$1
out=${file}.conv
Ndet=`grep "N_det =" $file | cut -d "=" -f 2`
Evar=`grep "E =" $file | cut -d "=" -f 2`
EPT2=`grep "E+rPT2 =" $file | cut -d "=" -f 2 | cut -d "+" -f 1`
err=`grep "E+rPT2 =" $file | cut -d "=" -f 2 | cut -d "+" -f 2 | cut -d "/" -f 2 | cut -d " " -f 5`
qp=$QP_ROOT
source ${qp}/quantum_package.rc
Ndetarray=[]
j=0
for i in $Ndet
qp_run print_e_conv $1
nstates=`cat ${1}/determinants/n_states`
echo $nstates
for i in `seq 1 $nstates`
do
Ndetarray[$j]=$i
let "j=j+1"
done
Nmax=${#Ndetarray[*]}
let "Nmax-=1"
Evararray=[]
j=0
for i in $Evar
do
Evararray[$j]=$i
let "j=j+1"
done
EPT2array=[]
j=0
for i in $EPT2
do
EPT2array[$j]=$i
let "j=j+1"
done
errarray=[]
j=0
for i in $err
do
errarray[$j]=$i
let "j=j+1"
done
echo "#Ndet E_var E+PT2 statistical error " > $out
for i in `seq 0 $Nmax`
do
echo ${Ndetarray[$i]} ${Evararray[$i]} ${EPT2array[$i]} ${errarray[$i]} >> $out
done
out=${1}.${i}.conv
cat << EOF > ${out}.plt
set term eps
set output "$out.eps"
@ -60,8 +19,28 @@ set log x
set xlabel "Number of determinants"
set ylabel "Total Energy (a.u.)"
plot "$out" w lp title "E_{var}", "$out" u 1:3:4 w errorlines title "E_{var} + PT2"
plot "$out" w lp title "E_{var} state $i", "$out" u 1:3 w lp title "E_{var} + PT2 state $i"
EOF
gnuplot ${out}.plt
rm ${out}.plt
done
for i in `seq 2 $nstates`
do
out=${1}.${i}.delta_e.conv
cat << EOF > ${out}.plt
set term eps
set output "$out.eps"
set log x
set xlabel "Number of determinants"
set ylabel "Energy difference (a.u.)"
plot "$out" w lp title "Delta E_{var} state $i", "$out" u 1:3 w lp title "Delta E_{var} + PT2 state $i"
EOF
gnuplot ${out}.plt
rm ${out}.plt
done

View File

@ -1,7 +1,7 @@
subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center)
implicit none
BEGIN_DOC
! subroutine that returs all integrals over r of type erf(mu_in * | r-C_center |)/| r-C_center |
! subroutine that returs all integrals over r of type erf(mu_in * | r-C_center | )/| r-C_center |
END_DOC
double precision, intent(in) :: mu_in,C_center(3)
double precision, intent(out) :: integrals_ao(ao_num,ao_num)
@ -19,7 +19,7 @@ double precision function NAI_pol_mult_erf_ao(i_ao,j_ao,mu_in,C_center)
implicit none
BEGIN_DOC
! computes the following integral :
! int[-infty;+infty] dr AO_i_ao (r) AO_j_ao(r) erf(mu_in * | r-C_center |)/| r-C_center |
! int[-infty;+infty] dr AO_i_ao (r) AO_j_ao(r) erf(mu_in * | r-C_center | )/| r-C_center |
END_DOC
integer, intent(in) :: i_ao,j_ao
double precision, intent(in) :: mu_in, C_center(3)

View File

@ -1 +0,0 @@

View File

@ -0,0 +1,71 @@
program print_e_conv
implicit none
BEGIN_DOC
! program that prints in a human readable format the convergence of the CIPSI algorithm
END_DOC
provide ezfio_filename
call routine
end
subroutine routine
implicit none
integer :: N_iter_tmp
integer :: i,istate
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
character*(128) :: filename
integer, allocatable :: n_det_tmp(:)
call ezfio_get_iterations_N_iter(N_iter_tmp)
print*,'N_iter_tmp = ',N_iter_tmp
double precision, allocatable :: e(:,:),pt2(:,:)
allocate(e(N_states, 100),pt2(N_states, 100),n_det_tmp(100))
call ezfio_get_iterations_energy_iterations(e)
call ezfio_get_iterations_pt2_iterations(pt2)
call ezfio_get_iterations_n_det_iterations(n_det_tmp)
do istate = 1, N_states
if (istate.lt.10)then
write (filename, "(I1)")istate
else
write (filename, "(I2)")istate
endif
print*,filename
output=trim(ezfio_filename)//'.'//trim(filename)//'.conv'
output=trim(output)
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,*)'# N_det E_var E_var + PT2'
do i = 1, N_iter_tmp
write(i_unit_output,'(I9,X,3(F16.10,X))')n_det_tmp(i),e(istate,i),e(istate,i) + pt2(istate,i)
enddo
enddo
if(N_states.gt.1)then
double precision, allocatable :: deltae(:,:),deltae_pt2(:,:)
allocate(deltae(N_states,100),deltae_pt2(N_states,100))
do i = 1, N_iter_tmp
do istate = 1, N_states
deltae(istate,i) = e(istate,i) - e(1,i)
deltae_pt2(istate,i) = e(istate,i) + pt2(istate,i) - (e(1,i) + pt2(1,i))
enddo
enddo
do istate = 2, N_states
if (istate.lt.10)then
write (filename, "(I1)")istate
else
write (filename, "(I2)")istate
endif
output=trim(ezfio_filename)//'.'//trim(filename)//'.delta_e.conv'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,*)'# N_det Delta E_var Delta (E_var + PT2)'
do i = 1, N_iter_tmp
write(i_unit_output,'(I9,X,100(F16.10,X))')n_det_tmp(i),deltae(istate,i),deltae_pt2(istate,i)
enddo
enddo
endif
end

View File

@ -59,6 +59,8 @@ subroutine routine
else
coef_2_2 = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * hij * hij )) /hij
endif
else
coef_2_2 = 0.d0
endif
call get_excitation(psi_det(1,1,1),psi_det(1,1,i),exc,degree,phase,N_int)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
@ -83,7 +85,7 @@ subroutine routine
print*,'hdouble = ',hdouble
print*,'hmono+hdouble = ',hmono+hdouble
print*,'hij = ',hij
else
else if(degree ==2)then
print*,'s1',s1
print*,'h1,p1 = ',h1,p1
print*,'s2',s2