Merge branch 'toto' of git://github.com/eginer/quantum_package

TODO

@@ -21,6 +21,8 @@
 
 # User doc:
 
+  * Videos: 
+    +) RHF
   * Renvoyer a la doc des modules : c'est pour les programmeurs au depart!
   * Mettre le mp2 comme exercice
 

docs/source/Makefile

@@ -3,3 +3,4 @@ default:
 	make -C ../ html
 clean:
 	make -C ../ clean
+	rm modules/*.rst

docs/source/_static/links.rst

@@ -39,6 +39,7 @@
 .. |OPAM|  replace:: `OPAM`_
 .. |Python|  replace:: `Python`_
 .. |qp| replace:: *Quantum Package*
+.. |resultsFile| replace:: `resultsFile`_
 .. |SLURM| replace:: `SLURM`_
 .. |ZeroMQ| replace:: `ZeroMQ`_
 

docs/source/index.rst

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

docs/source/intro/intro.rst

@@ -11,8 +11,18 @@ The |qp|
 What it is
 ==========
 
-The |qp| is an open-source programming environment for quantum chemistry,
+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|.  

docs/source/modules/perturbation.rst

@@ -245,6 +245,20 @@ Subroutines / functions
 
 
 
+.. c:function:: perturb_buffer_by_mono_epstein_nesbet_2x2_no_ci_diag
+
+    .. code:: text
+
+        subroutine perturb_buffer_by_mono_epstein_nesbet_2x2_no_ci_diag(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 ``epstein_nesbet_2x2_no_ci_diag`` to the buffer of determinants generated in the H_apply routine.
+
+
+
+
+
 .. c:function:: perturb_buffer_by_mono_h_core
 
     .. code:: text
@@ -343,6 +357,20 @@ Subroutines / functions
 
 
 
+.. c:function:: perturb_buffer_epstein_nesbet_2x2_no_ci_diag
+
+    .. code:: text
+
+        subroutine perturb_buffer_epstein_nesbet_2x2_no_ci_diag(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 ``epstein_nesbet_2x2_no_ci_diag`` to the buffer of determinants generated in the H_apply routine.
+
+
+
+
+
 .. c:function:: perturb_buffer_h_core
 
     .. code:: text
@@ -405,7 +433,7 @@ Subroutines / functions
 
         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_309`
+    File: :file:`pt2_equations.irp.f_template_360`
 
     Dummy perturbation to add all connected determinants.
 
@@ -419,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_309`
+    File: :file:`pt2_equations.irp.f_template_360`
 
     Compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution for the various N_st states. 
 
@@ -439,7 +467,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_309`
+    File: :file:`pt2_equations.irp.f_template_360`
 
     Computes the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution for the various N_st states. 
 
@@ -453,6 +481,28 @@ Subroutines / functions
 
 
 
+.. c:function:: pt2_epstein_nesbet_2x2_no_ci_diag
+
+    .. code:: text
+
+        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_360`
+
+    compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution 
+
+    for the various N_st states. 
+
+    e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  ) 
+
+    c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert> 
+
+
+
+
+
+
+
 .. c:function:: pt2_h_core
 
     .. code:: text
@@ -481,7 +531,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_309`
+    File: :file:`pt2_equations.irp.f_template_360`
 
     Computes the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution for the various N_st states. 
 
@@ -501,7 +551,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_309`
+    File: :file:`pt2_equations.irp.f_template_360`
 
     Computes the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution for the various N_st states. 
 
@@ -521,7 +571,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_309`
+    File: :file:`pt2_equations.irp.f_template_360`
 
     Computes the QDPT first order coefficient and second order energetic contribution for the various N_st states. 
 

docs/source/modules/tools.rst

@@ -73,6 +73,20 @@ Subroutines / functions
 
 
 
+.. c:function:: print_e_conv
+
+    .. code:: text
+
+        subroutine print_e_conv
+
+    File: :file:`print_e_conv.irp.f`
+
+    program that prints in a human readable format the convergence of the CIPSI algorithm
+
+
+
+
+
 .. c:function:: print_wf
 
     .. code:: text

docs/source/programmers_guide/index_providers.rst

@@ -1156,6 +1156,7 @@ Index of Subroutines/Functions
 * :c:func:`perturb_buffer_by_mono_dummy`
 * :c:func:`perturb_buffer_by_mono_epstein_nesbet`
 * :c:func:`perturb_buffer_by_mono_epstein_nesbet_2x2`
+* :c:func:`perturb_buffer_by_mono_epstein_nesbet_2x2_no_ci_diag`
 * :c:func:`perturb_buffer_by_mono_h_core`
 * :c:func:`perturb_buffer_by_mono_moller_plesset`
 * :c:func:`perturb_buffer_by_mono_moller_plesset_general`
@@ -1163,12 +1164,14 @@ Index of Subroutines/Functions
 * :c:func:`perturb_buffer_dummy`
 * :c:func:`perturb_buffer_epstein_nesbet`
 * :c:func:`perturb_buffer_epstein_nesbet_2x2`
+* :c:func:`perturb_buffer_epstein_nesbet_2x2_no_ci_diag`
 * :c:func:`perturb_buffer_h_core`
 * :c:func:`perturb_buffer_moller_plesset`
 * :c:func:`perturb_buffer_moller_plesset_general`
 * :c:func:`perturb_buffer_qdpt`
 * :c:func:`primitive_value`
 * :c:func:`print_det`
+* :c:func:`print_e_conv`
 * :c:func:`print_extrapolated_energy`
 * :c:func:`print_generators_bitmasks_holes`
 * :c:func:`print_generators_bitmasks_holes_for_one_generator`
@@ -1184,6 +1187,7 @@ Index of Subroutines/Functions
 * :c:func:`pt2_dummy`
 * :c:func:`pt2_epstein_nesbet`
 * :c:func:`pt2_epstein_nesbet_2x2`
+* :c:func:`pt2_epstein_nesbet_2x2_no_ci_diag`
 * :c:func:`pt2_find_sample`
 * :c:func:`pt2_find_sample_lr`
 * :c:func:`pt2_h_core`

docs/source/programmers_guide/programming.rst

@@ -11,7 +11,7 @@ in order to write a new program is the name of the required |IRPF90| entities
 which may already exist in other modules.  For example, writing a program which
 prints the Hartree-Fock energy is as simple as:
 
-.. code:: irpf90
+.. code:: fortran
 
     program print_hf_energy
       implicit none

docs/source/users_guide/qp_plugins.rst

@@ -25,15 +25,15 @@ Usage
 
     List all the available plugins.
 
-.. option:: list -i 
+.. option:: -i 
 
     List all the *installed* plugins.
 
-.. option:: list -u 
+.. option:: -u
 
     List all the *uninstalled* plugins.
 
-.. option:: list -q 
+.. option:: -q 
 
     List all the downloaded repositories.
 
@@ -53,12 +53,12 @@ Usage
 
     Uninstall the plugin ``plugin_name``.
 
-.. option:: create -n <plugin_name>
+.. option:: -n <plugin_name>
 
-    Create a new plugin named ``plugin_name`` in local repository.
+    Create a new plugin named ``plugin_name`` (in local repository by default).
 
-.. option:: create -n <plugin_name> -r <repository>
+.. option:: -r <repository>
 
-    Create a new plugin named ``plugin_name`` in the specified repository.
+    Specify in which repository the new plugin will be created.
 
 

docs/source/users_guide/quickstart.rst

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

scripts/qp_e_conv_fci

@@ -1,67 +1,46 @@
 #!/bin/bash
 file=$1
-out=${file}.conv
 
- Ndet=`grep "N_det             =" $file | cut -d "=" -f 2` 
+qp=$QP_ROOT
- Evar=`grep "E               =" $file | cut -d "=" -f 2` 
+source ${qp}/quantum_package.rc
- 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`
 
-Ndetarray=[]
+qp_run print_e_conv $1 
-j=0
+nstates=`cat ${1}/determinants/n_states`
-for i in $Ndet
+echo $nstates
+
+
+for i in `seq 1 $nstates`
  do 
-  Ndetarray[$j]=$i
+ out=${1}.${i}.conv
-  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
-
 cat << EOF > ${out}.plt
-set term eps
+set term pdf
-set output "$out.eps"
+set output "$out.pdf"
 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 pdf
+set output "$out.pdf"
+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

src/iterations/NEED

@@ -1 +0,0 @@
-

src/kohn_sham_rs/61.rsks.bats

@@ -6,33 +6,33 @@ source $QP_ROOT/tests/bats/common.bats.sh
 function run() {
   thresh=1.e-8
   qp_edit -c $1
+  functional=$2
   ezfio set_file $1
   ezfio set scf_utils thresh_scf 1.e-10
-  ezfio set dft_keywords exchange_functional "short_range_PBE" 
+  ezfio set dft_keywords exchange_functional $functional
-  ezfio set dft_keywords correlation_functional "short_range_PBE" 
+  ezfio set dft_keywords correlation_functional $functional
   ezfio set ao_two_e_erf_ints mu_erf 0.5 
   ezfio set becke_numerical_grid grid_type_sgn 1 
   qp_run rs_ks_scf $1
   energy="$(ezfio get kohn_sham_rs energy)"
-  eq $energy $2 $thresh
+  eq $energy $3 $thresh
 }
 
 
-@test "H3COH" { # 11.4566s
+@test "H3COH" {
-  run h3coh.ezfio -115.50238225208
+  run h3coh.ezfio short_range_PBE -115.50238225208
 }
 
-@test "N2" { # 18.2364s
+@test "HCN" {
-  run n2.ezfio -109.404692225719
+  run hcn.ezfio short_range_PBE -93.26674673761752
 }
 
-@test "HCN" { # 28.801s
+@test "N2" {
-  run hcn.ezfio -93.26674673761752
+  run n2.ezfio  short_range_PBE -109.404692225719
 }
 
-@test "SiH2_3B1" { # 82.3904s
+@test "SiH2_3B1" {
-  [[ -n $TRAVIS ]] && skip
+  run sih2_3b1.ezfio short_range_LDA -289.4398733527755
-  run sih2_3b1.ezfio -290.372258160809
 }
 
 

src/mo_basis/track_orb.irp.f

@@ -47,11 +47,9 @@ subroutine reorder_active_orb
 
  double precision :: x
  integer :: i1,i2
- print*, 'swapping the active MOs'
  do j = 1, n_act_orb
   i1 = list_act(j)
   i2 = index_active_orb(j)
-  print*, i1,i2
   do i=1,ao_num
     x = mo_coef(i,i1)
     mo_coef(i,i1) = mo_coef(i,i2)

src/perturbation/pt2_equations.irp.f

@@ -146,6 +146,57 @@ end
 
 
 
+subroutine pt2_epstein_nesbet_2x2_no_ci_diag($arguments)
+  use bitmasks
+  implicit none
+  $declarations
+  
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
+  !
+  ! for the various N_st states.
+  !
+  ! e_2_pert(i) = 0.5 * (( <det_pert|H|det_pert> -  E(i) )  - sqrt( ( <det_pert|H|det_pert> -  E(i)) ^2 + 4 <psi(i)|H|det_pert>^2  )
+  !
+  ! c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem_fock,delta_e, h
+  double precision               :: i_H_psi_array(N_st)
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  PROVIDE psi_energy
+
+   call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
+  
+  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  do i =1,N_st
+    if (i_H_psi_array(i) /= 0.d0) then
+      delta_e = h - psi_energy(i)
+      if (delta_e > 0.d0) then
+        e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      else
+        e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
+      endif
+      if (dabs(i_H_psi_array(i)) > 1.d-6) then
+        c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
+      else
+        c_pert(i) = 0.d0
+      endif
+      H_pert_diag(i) = h*c_pert(i)*c_pert(i)
+    else
+      e_2_pert(i) = 0.d0
+      c_pert(i) = 0.d0
+      H_pert_diag(i) = 0.d0
+    endif
+  enddo
+
+end
+
+
+
 subroutine pt2_moller_plesset ($arguments)
   use bitmasks
   implicit none

src/tools/print_e_conv.irp.f

@@ -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

src/tools/print_wf.irp.f

@@ -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