Merge branch 'master' of github.com:LCPQ/quantum_package

Conflicts:
	src/CISD/README.rst
	src/CISD/SC2.irp.f
	src/CISD/cisd.irp.f
	src/Dets/README.rst
	src/Dets/determinants.irp.f
	src/Dets/slater_rules.irp.f
	src/Perturbation/README.rst

scripts/run_tests.py

@@ -23,10 +23,13 @@ if __name__ == '__main__':
     env = os.environ
     
     verbosity = 1
+    if os.environ["OMP_NUM_THREADS"]=="1":
       try:
         nproc = int(subprocess.check_output("cat /proc/cpuinfo | grep processor | wc -l", shell=True))
       except:
         nproc=4
+    else:
+      nproc=1
 
     testfiles = []
     for f in os.listdir(os.getcwd()):

src/CISD/README.rst

@@ -47,7 +47,10 @@ Documentation
   .br
   Initial guess vectors are not necessarily orthonormal
 
-`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/SC2.irp.f#L137>`_
+`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/SC2.irp.f#L143>`_
+  Undocumented
+
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd_sc2.irp.f#L1>`_
   Undocumented
 
 

src/CISD/SC2.irp.f

@@ -36,6 +36,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
   double precision               :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0
   double precision               :: e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),hij_double(sze)
   double precision               :: e_corr_double_before,accu,cpu_2,cpu_1
+  integer                        :: degree_exc(sze)
   integer                        :: i_ok
   
   call write_time(output_CISD)
@@ -59,6 +60,7 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
   e_corr_double = 0.d0
   do i = 1, sze
     call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint)
+    degree_exc(i) = degree
     if(degree==0)then
       index_hf=i
     else if (degree == 2)then
@@ -81,23 +83,28 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
   enddo
   e_corr = e_corr * inv_c0
   e_corr_double = e_corr_double * inv_c0
-  
   converged = .False.
   e_corr_double_before = e_corr_double
   iter = 0
   do while (.not.converged)
-    
     iter +=1
-    !  call cpu_time(cpu_1)
     do i=1,sze
       H_jj_dressed(i) = H_jj_ref(i)
       if (i==index_hf)cycle
       accu = 0.d0
+      if(degree_exc(i)==1)then
         do j=1,N_double
-        call repeat_excitation(dets_in(1,1,i),ref_bitmask,dets_in(1,1,index_double(j)),i_ok,Nint)
+          call get_excitation_degree(dets_in(1,1,i),dets_in(1,1,index_double(j)),degree,N_int)
-        if (i_ok==1)cycle! you check if the excitation is possible
+          if (degree<=2)cycle
           accu += e_corr_array(j)
         enddo
+      else
+        do j=1,N_double
+          call get_excitation_degree(dets_in(1,1,i),dets_in(1,1,index_double(j)),degree,N_int)
+          if (degree<=3)cycle
+          accu += e_corr_array(j)
+        enddo
+      endif
       H_jj_dressed(i) += accu
     enddo
     
@@ -131,7 +138,6 @@ subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint)
   
   call write_time(output_CISD)
   
-  
 end
 
 subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)
@@ -180,4 +186,3 @@ subroutine repeat_excitation(key_in,key_1,key_2,i_ok,Nint)
     return
   endif
 end
-

src/CISD/cisd.irp.f

@@ -5,6 +5,9 @@ program cisd
   PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
 
   double precision :: pt2(10), norm_pert(10), H_pert_diag
+  double precision,allocatable :: H_jj(:)
+  double precision :: diag_h_mat_elem
+  integer,allocatable :: iorder(:)
 
 ! N_states = 3
 ! touch N_states
@@ -13,10 +16,20 @@ program cisd
   print *,  'N_det = ', N_det
   print *,  'N_states = ', N_states
   psi_coef = - 1.d-4
+  allocate(H_jj(n_det),iorder(n_det))
+  do i = 1, N_det
+   H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
+   iorder(i) = i
+  enddo
+  call dsort(H_jj,iorder,n_det)
+
   do k=1,N_states
-    psi_coef(k,k) = 1.d0
+    psi_coef(iorder(k),k) = 1.d0
   enddo
   call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
+  do i = 1, N_states
+   print*,'eigvalues(i) = ',eigvalues(i)
+  enddo
 
   print *,  '---'
   print *,  'HF:', HF_energy
@@ -25,5 +38,9 @@ program cisd
    print *,  'energy(i)     = ',eigvalues(i) + nuclear_repulsion
    print *,  'E_corr        = ',eigvalues(i) - ref_bitmask_energy
   enddo
+! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
+! do i = 1, N_states
+!  print*,'eigvalues(i) = ',eigvalues(i)
+! enddo
   deallocate(eigvalues,eigvectors)
 end

src/CISD/cisd_lapack.irp.f

@@ -1,6 +1,7 @@
 program cisd
   implicit none
   integer :: i
+  PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
   call H_apply_cisd
   double precision, allocatable :: eigvalues(:),eigvectors(:,:)
   allocate(eigvalues(n_det),eigvectors(n_det,n_det))

src/CISD/cisd_sc2.irp.f

@@ -2,9 +2,12 @@ program cisd
   implicit none
   integer :: i,k
   double precision, allocatable  :: eigvalues(:),eigvectors(:,:)
-  PROVIDE ref_bitmask_energy
+  PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
 
   double precision :: pt2(10), norm_pert(10), H_pert_diag
+  double precision,allocatable :: H_jj(:)
+  double precision :: diag_h_mat_elem
+  integer,allocatable :: iorder(:)
 
 ! N_states = 3
 ! touch N_states
@@ -13,10 +16,20 @@ program cisd
   print *,  'N_det = ', N_det
   print *,  'N_states = ', N_states
   psi_coef = - 1.d-4
+  allocate(H_jj(n_det),iorder(n_det))
+  do i = 1, N_det
+   H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
+   iorder(i) = i
+  enddo
+  call dsort(H_jj,iorder,n_det)
+
   do k=1,N_states
-    psi_coef(k,k) = 1.d0
+    psi_coef(iorder(k),k) = 1.d0
   enddo
   call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
+  do i = 1, N_states
+   print*,'eigvalues(i) = ',eigvalues(i)
+  enddo
 
   print *,  '---'
   print *,  'HF:', HF_energy
@@ -26,5 +39,8 @@ program cisd
    print *,  'E_corr        = ',eigvalues(i) - ref_bitmask_energy
   enddo
   call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
+  do i = 1, N_states
+   print*,'eigvalues(i) = ',eigvalues(i)
+  enddo
   deallocate(eigvalues,eigvectors)
 end

src/CIS_dressed/Makefile

@@ -0,0 +1,8 @@
+default: all
+
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/CIS_dressed/NEEDED_MODULES

@@ -0,0 +1 @@
+AOs BiInts Bitmask Dets Electrons Ezfio_files Hartree_Fock MonoInts MOs Nuclei Output Utils DensityMatrix

src/CIS_dressed/README.rst

@@ -24,5 +24,210 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`coefs_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L105>`_
+  the first states of the CIS matrix
+
+`diag_elements <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L317>`_
+  Array of the energy of the CIS determinants ordered in the CIS matrix
+
+`eigenvalues_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L104>`_
+  the first states of the CIS matrix
+
+`eigenvalues_cis_dress_d <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L141>`_
+  The first states of the CIS matrix dressed by the doubles
+
+`eigenvalues_cis_dress_d_dt <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L198>`_
+  The first states of the CIS matrix dressed by the doubles and the disconnected triples
+
+`eigenvectors_cis_dress_d <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L143>`_
+  The first states of the CIS matrix dressed by the doubles
+
+`eigenvectors_cis_dress_d_dt <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L200>`_
+  The first states of the CIS matrix dressed by the doubles and the disconnected triples
+
+`h_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L8>`_
+  key of the CIS-matrix
+
+`overlap_d <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L144>`_
+  The first states of the CIS matrix dressed by the doubles
+
+`overlap_ddt <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L201>`_
+  The first states of the CIS matrix dressed by the doubles and the disconnected triples
+
+`psi_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L3>`_
+  key of the CIS-matrix
+
+`psi_cis_adress <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L7>`_
+  key of the CIS-matrix
+
+`psi_cis_holes <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L4>`_
+  key of the CIS-matrix
+
+`psi_cis_particl <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L5>`_
+  key of the CIS-matrix
+
+`psi_cis_spin <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L6>`_
+  key of the CIS-matrix
+
+`s_2_cis_dress_d <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L142>`_
+  The first states of the CIS matrix dressed by the doubles
+
+`s_2_cis_dress_d_dt <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/CIS_providers.irp.f#L199>`_
+  The first states of the CIS matrix dressed by the doubles and the disconnected triples
+
+`diag_elements_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L420>`_
+  Array of the energy of the CIS determinants sorted by energy and
+  Index in the CIS matrix
+
+`diag_elements_sorted_index <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L421>`_
+  Array of the energy of the CIS determinants sorted by energy and
+  Index in the CIS matrix
+
+`diexcitation <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L1563>`_
+  Undocumented
+
+`dress_by_doubles <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L672>`_
+  Undocumented
+
+`dress_t_con <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L1048>`_
+  Generating all the Triples and, in the loops, the connected Singles
+
+`dress_t_discon <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L621>`_
+  Calculation of the dressing by the disconnected Triples, via the impossible
+
+`dress_t_discon_array_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L622>`_
+  Calculation of the dressing by the disconnected Triples, via the impossible
+
+`eigenvalues_dressed_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L451>`_
+  The first states of the dressed CIS matrix
+
+`en2_corr_energy <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L142>`_
+  Calculation of the EN2 correlation energy (EN2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN)
+
+`h_imp <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L3>`_
+  Calculation of the MP2 correlation energy (MP2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp, h_imp, hp_imp)
+
+`h_imp_en <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L144>`_
+  Calculation of the EN2 correlation energy (EN2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN)
+
+`hp_imp <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L4>`_
+  Calculation of the MP2 correlation energy (MP2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp, h_imp, hp_imp)
+
+`hp_imp_en <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L145>`_
+  Calculation of the EN2 correlation energy (EN2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN)
+
+`mp2_corr_energy <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L1>`_
+  Calculation of the MP2 correlation energy (MP2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp, h_imp, hp_imp)
+
+`p_imp <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L2>`_
+  Calculation of the MP2 correlation energy (MP2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp, h_imp, hp_imp)
+
+`p_imp_en <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L143>`_
+  Calculation of the EN2 correlation energy (EN2_corr_energy)
+  and calculation of the contribution of the disconnected Triples on the
+  Singles, via the impossible (p_imp_EN, h_imp_EN, hp_imp_EN)
+
+`size_psi_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/MP2.irp.f#L407>`_
+  Definition of the size of the CIS vector
+
+`get_dm_from_psi <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/density_matrix_suroutine.irp.f#L2>`_
+  Alpha and beta one-body density matrix
+  .br
+  dets_in   ::  bitsrings corresponding to the determinants in the wave function
+  .br
+  u_in      ::  coefficients of the wave function
+  .br
+  sze       ::  number of determinants in the wave function
+  .br
+  dim_in    ::  physical dimension of the array u_in and dets_in
+  .br
+  Nint      ::  should be equal to N_int
+  .br
+  dm_alpha  ::  alpha one body density matrix
+  .br
+  dm_beta   ::  beta  one body density matrix
+
+`cis_states_properties <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/natural_particl_orbitals.irp.f#L2>`_
+  properties of the natural orbital of the particle of the various n_state_cis eigenvectors of the CIS matrix
+  .br
+  You first build the density matrix of the one eigenvector and you take off the Hartree Fock density matrix
+  .br
+  particl(i,j)(state = k) ==   dm(i,j)(Hartree Fock)  - dm(i,j)(state = k)
+  .br
+  you diagonalize particl(i,j) and the first eigenvector is the natural orbital corresponding to the particl
+  .br
+  that is specific to the excitation in the CIS state
+  .br
+  particle_natural_orb_CIS_properties(i,1) = <phi_i|x|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,2) = <phi_i|y|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,3) = <phi_i|z|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,5) = <phi_i|x^2|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,6) = <phi_i|y^2|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,7) = <phi_i|z^2|phi_i>
+  .br
+  CIS_states_properties(i,1:6) = the same but for the hole state i
+
+`get_properties_from_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/natural_particl_orbitals.irp.f#L71>`_
+  Undocumented
+
+`particle_natural_orb_cis_properties <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/natural_particl_orbitals.irp.f#L1>`_
+  properties of the natural orbital of the particle of the various n_state_cis eigenvectors of the CIS matrix
+  .br
+  You first build the density matrix of the one eigenvector and you take off the Hartree Fock density matrix
+  .br
+  particl(i,j)(state = k) ==   dm(i,j)(Hartree Fock)  - dm(i,j)(state = k)
+  .br
+  you diagonalize particl(i,j) and the first eigenvector is the natural orbital corresponding to the particl
+  .br
+  that is specific to the excitation in the CIS state
+  .br
+  particle_natural_orb_CIS_properties(i,1) = <phi_i|x|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,2) = <phi_i|y|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,3) = <phi_i|z|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,5) = <phi_i|x^2|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,6) = <phi_i|y^2|phi_i>
+  .br
+  particle_natural_orb_CIS_properties(i,7) = <phi_i|z^2|phi_i>
+  .br
+  CIS_states_properties(i,1:6) = the same but for the hole state i
+
+`mp2_dressing <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/options.irp.f#L55>`_
+  Number of states asked for the CIS vector
+
+`n_act_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/options.irp.f#L37>`_
+  Number of states asked for the CIS vector
+
+`n_core_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/options.irp.f#L19>`_
+  Number of states asked for the CIS vector
+
+`n_state_cis <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/options.irp.f#L1>`_
+  Number of states asked for the CIS vector
+
+`standard_doubles <http://github.com/LCPQ/quantum_package/tree/master/src/CIS_dressed/options.irp.f#L73>`_
+  Number of states asked for the CIS vector
+
 
 

src/DensityMatrix/README.rst

@@ -8,37 +8,37 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`iunit_two_body_dm_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/BEGIN_PROVIDER [ integer, iunit_two_body_dm_aa ]/;">`_
+`iunit_two_body_dm_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L2>`_
   Temporary files for 2-body dm calculation
 
-`iunit_two_body_dm_ab <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ integer, iunit_two_body_dm_ab ]/;">`_
+`iunit_two_body_dm_ab <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L3>`_
   Temporary files for 2-body dm calculation
 
-`iunit_two_body_dm_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ integer, iunit_two_body_dm_bb ]/;">`_
+`iunit_two_body_dm_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L4>`_
   Temporary files for 2-body dm calculation
 
-`one_body_dm_a <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, one_body_dm_a, (mo_tot_num_align,mo_tot_num) ]/;">`_
+`one_body_dm_a <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L213>`_
   Alpha and beta one-body density matrix
 
-`one_body_dm_b <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, one_body_dm_b, (mo_tot_num_align,mo_tot_num) ]/;">`_
+`one_body_dm_b <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L214>`_
   Alpha and beta one-body density matrix
 
-`two_body_dm_diag_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/BEGIN_PROVIDER [ double precision, two_body_dm_diag_aa, (mo_tot_num_align,mo_tot_num)]/;">`_
+`two_body_dm_diag_aa <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L170>`_
   diagonal part of the two body density matrix
 
-`two_body_dm_diag_ab <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, two_body_dm_diag_ab, (mo_tot_num_align,mo_tot_num)]/;">`_
+`two_body_dm_diag_ab <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L172>`_
   diagonal part of the two body density matrix
 
-`two_body_dm_diag_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L/&BEGIN_PROVIDER [ double precision, two_body_dm_diag_bb, (mo_tot_num_align,mo_tot_num)]/;">`_
+`two_body_dm_diag_bb <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/density_matrix.irp.f#L171>`_
   diagonal part of the two body density matrix
 
-`det_coef_provider <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L/&BEGIN_PROVIDER [ double precision , det_coef_provider, (det_num) ]/;">`_
+`det_coef_provider <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L8>`_
   Undocumented
 
-`det_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L/BEGIN_PROVIDER [integer, det_num]/;">`_
+`det_num <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L3>`_
   Undocumented
 
-`det_provider <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L/BEGIN_PROVIDER [ integer(bit_kind), det_provider, (N_int,2,det_num)]/;">`_
+`det_provider <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix/det_num.irp.f#L7>`_
   Undocumented
 
 
@@ -61,5 +61,4 @@ Needed Modules
 * `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
 * `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
 * `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
-* `DensityMatrix <http://github.com/LCPQ/quantum_package/tree/master/src/DensityMatrix>`_
 

src/Dets/README.rst

@@ -161,9 +161,34 @@ Documentation
 
 `filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L2>`_
   Filters out the determinants that are not connected by H
+  .br
+  returns the array idx which contains the index of the
+  .br
+  determinants in the array key1 that interact
+  .br
+  via the H operator with key2.
+  .br
+  idx(0) is the number of determinants that interact with key1
 
-`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L86>`_
+`filter_connected_i_h_psi0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L94>`_
-  Undocumented
+  returns the array idx which contains the index of the
+  .br
+  determinants in the array key1 that interact
+  .br
+  via the H operator with key2.
+  .br
+  idx(0) is the number of determinants that interact with key1
+
+`filter_connected_i_h_psi0_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/filter_connected.irp.f#L193>`_
+  standard filter_connected_i_H_psi but returns in addition
+  .br
+  the array of the index of the non connected determinants to key1
+  .br
+  in order to know what double excitation can be repeated on key1
+  .br
+  idx_repeat(0) is the number of determinants that can be used
+  .br
+  to repeat the excitations
 
 `get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L1>`_
   Returns <S^2>
@@ -177,10 +202,10 @@ Documentation
 `s_z2_sz <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/s2.irp.f#L37>`_
   Undocumented
 
-`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L671>`_
+`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L721>`_
   Needed for diag_H_mat_elem
 
-`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L716>`_
+`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L766>`_
   Needed for diag_H_mat_elem
 
 `decode_exc <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L76>`_
@@ -190,7 +215,7 @@ Documentation
   s1,s2 : Spins (1:alpha, 2:beta)
   degree : Degree of excitation
 
-`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L609>`_
+`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L659>`_
   Computes <i|H|i>
 
 `get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L141>`_
@@ -202,16 +227,16 @@ Documentation
 `get_excitation_degree <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L1>`_
   Returns the excitation degree between two determinants
 
-`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L525>`_
+`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L575>`_
   Applies get_excitation_degree to an array of determinants
 
 `get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L274>`_
   Returns the excitation operator between two singly excited determinants and the phase
 
-`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L764>`_
+`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L814>`_
   Returns a list of occupation numbers from a bitstring
 
-`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L780>`_
+`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L830>`_
   Computes v_0 = H|u_0>
   .br
   n : number of determinants
@@ -222,7 +247,20 @@ Documentation
   Returns <i|H|j> where i and j are determinants
 
 `i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L491>`_
-  Undocumented
+  <key|H|psi> for the various Nstate
+
+`i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/slater_rules.irp.f#L527>`_
+  <key|H|psi> for the various Nstate
+  .br
+  returns in addition
+  .br
+  the array of the index of the non connected determinants to key1
+  .br
+  in order to know what double excitation can be repeated on key1
+  .br
+  idx_repeat(0) is the number of determinants that can be used
+  .br
+  to repeat the excitations
 
 `h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/utils.irp.f#L1>`_
   H matrix on the basis of the slater deter;inants defined by psi_det

src/Dets/filter_connected.irp.f

@@ -4,6 +4,14 @@ subroutine filter_connected(key1,key2,Nint,sze,idx)
   implicit none
   BEGIN_DOC
   ! Filters out the determinants that are not connected by H
+  !
+  ! returns the array idx which contains the index of the 
+  !
+  ! determinants in the array key1 that interact 
+  !
+  ! via the H operator with key2.
+  !
+  ! idx(0) is the number of determinants that interact with key1
   END_DOC
   integer, intent(in)            :: Nint, sze
   integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
@@ -85,6 +93,15 @@ end
 
 subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
   use bitmasks
+  BEGIN_DOC
+  ! returns the array idx which contains the index of the 
+  !
+  ! determinants in the array key1 that interact 
+  !
+  ! via the H operator with key2.
+  !
+  ! idx(0) is the number of determinants that interact with key1
+  END_DOC
   implicit none
   integer, intent(in)            :: Nint, sze
   integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
@@ -173,3 +190,215 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
   idx(0) = l-1
 end
 
+subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat)
+  use bitmasks
+  BEGIN_DOC
+  ! standard filter_connected_i_H_psi but returns in addition
+  !
+  ! the array of the index of the non connected determinants to key1
+  !
+  ! in order to know what double excitation can be repeated on key1
+  !
+  ! idx_repeat(0) is the number of determinants that can be used 
+  ! 
+  ! to repeat the excitations 
+  END_DOC
+  implicit none
+  integer, intent(in)            :: Nint, sze
+  integer(bit_kind), intent(in)  :: key1(Nint,2,sze)
+  integer(bit_kind), intent(in)  :: key2(Nint,2)
+  integer, intent(out)           :: idx(0:sze)
+  integer, intent(out)           :: idx_repeat(0:sze)
+  
+  integer                        :: i,l,l_repeat
+  integer                        :: degree_x2
+
+  ASSERT (Nint > 0)
+  ASSERT (Nint == N_int)
+  ASSERT (sze > 0)
+  
+  l=1
+  l_repeat=1
+  call get_excitation_degree(ref_bitmask,key2,degree,Nint)
+  integer :: degree
+  
+  if(degree == 2)then
+   if (Nint==1) then
+     
+ 
+      !DIR$ LOOP COUNT (1000)
+      do i=1,sze
+        degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) +             &
+            popcnt(xor( key1(1,2,i), key2(1,2)))
+        if (degree_x2 < 5) then
+          if(degree_x2 .ne. 0)then
+            idx(l) = i
+            l = l+1
+          endif
+        elseif(degree>6)then
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+        endif
+      enddo
+     
+   else if (Nint==2) then
+     
+     !DIR$ LOOP COUNT (1000)
+     do i=1,sze
+       degree_x2 =  popcnt(xor( key1(1,1,i), key2(1,1))) +            &
+           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
+           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
+           popcnt(xor( key1(2,2,i), key2(2,2)))
+       if (degree_x2 < 5) then
+         if(degree_x2 .ne. 0)then
+           idx(l) = i
+           l = l+1
+         endif
+       elseif(degree>6)then
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       endif
+     enddo
+     
+   else if (Nint==3) then
+     
+     !DIR$ LOOP COUNT (1000)
+     do i=1,sze
+       degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) +             &
+           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
+           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
+           popcnt(xor( key1(2,2,i), key2(2,2))) +                     &
+           popcnt(xor( key1(3,1,i), key2(3,1))) +                     &
+           popcnt(xor( key1(3,2,i), key2(3,2)))
+       if (degree_x2 < 5) then
+         if(degree_x2 .ne. 0)then
+           idx(l) = i
+           l = l+1
+         endif
+       elseif(degree>6)then
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       endif
+     enddo
+     
+   else
+     
+     !DIR$ LOOP COUNT (1000)
+     do i=1,sze
+       degree_x2 = 0
+       !DEC$ LOOP COUNT MIN(4)
+       do l=1,Nint
+         degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +&
+             popcnt(xor( key1(l,2,i), key2(l,2)))
+         if (degree_x2 > 4) then
+           exit
+         endif
+       enddo
+       if (degree_x2 <= 5) then
+         if(degree_x2 .ne. 0)then
+           idx(l) = i
+           l = l+1
+         endif
+       elseif(degree>6)then
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       endif
+     enddo
+     
+   endif
+  elseif(degree==1)then
+   if (Nint==1) then
+     
+ 
+      !DIR$ LOOP COUNT (1000)
+      do i=1,sze
+        degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) +             &
+            popcnt(xor( key1(1,2,i), key2(1,2)))
+        if (degree_x2 < 5) then
+          if(degree_x2 .ne. 0)then
+            idx(l) = i
+            l = l+1
+          endif
+        else
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+        endif
+      enddo
+     
+   else if (Nint==2) then
+     
+     !DIR$ LOOP COUNT (1000)
+     do i=1,sze
+       degree_x2 =  popcnt(xor( key1(1,1,i), key2(1,1))) +            &
+           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
+           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
+           popcnt(xor( key1(2,2,i), key2(2,2)))
+       if (degree_x2 < 5) then
+         if(degree_x2 .ne. 0)then
+           idx(l) = i
+           l = l+1
+         endif
+       else
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       endif
+     enddo
+     
+   else if (Nint==3) then
+     
+     !DIR$ LOOP COUNT (1000)
+     do i=1,sze
+       degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) +             &
+           popcnt(xor( key1(1,2,i), key2(1,2))) +                     &
+           popcnt(xor( key1(2,1,i), key2(2,1))) +                     &
+           popcnt(xor( key1(2,2,i), key2(2,2))) +                     &
+           popcnt(xor( key1(3,1,i), key2(3,1))) +                     &
+           popcnt(xor( key1(3,2,i), key2(3,2)))
+       if (degree_x2 < 5) then
+         if(degree_x2 .ne. 0)then
+           idx(l) = i
+           l = l+1
+         endif
+        else
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       endif
+     enddo
+     
+   else
+     
+     !DIR$ LOOP COUNT (1000)
+     do i=1,sze
+       degree_x2 = 0
+       !DEC$ LOOP COUNT MIN(4)
+       do l=1,Nint
+         degree_x2 = degree_x2+ popcnt(xor( key1(l,1,i), key2(l,1))) +&
+             popcnt(xor( key1(l,2,i), key2(l,2)))
+         if (degree_x2 > 4) then
+           exit
+         endif
+       enddo
+       if (degree_x2 <= 5) then
+         if(degree_x2 .ne. 0)then
+           idx(l) = i
+           l = l+1
+         endif
+        else
+         idx_repeat(l_repeat) = i
+         l_repeat = l_repeat + 1
+       endif
+     enddo
+     
+   endif
+
+  else 
+   print*,'more than a double excitation, can not apply the '
+   print*,'SC2 dressing of the diagonal element .....'
+   print*,'stop !!'
+   print*,'degree = ',degree
+   stop
+  endif
+  idx(0) = l-1
+  idx_repeat(0) = l_repeat-1
+end
+

src/Dets/slater_rules.irp.f

@@ -502,6 +502,9 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
   integer                        :: exc(0:2,2,2)
   double precision               :: hij
   integer                        :: idx(0:Ndet)
+  BEGIN_DOC
+  ! <key|H|psi> for the various Nstate
+  END_DOC
   
   ASSERT (Nint > 0)
   ASSERT (N_int == Nint)
@@ -521,6 +524,53 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
 end
 
 
+subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat)
+  use bitmasks
+  BEGIN_DOC
+  ! <key|H|psi> for the various Nstate
+  ! 
+  ! returns in addition 
+  !
+  ! the array of the index of the non connected determinants to key1
+  !
+  ! in order to know what double excitation can be repeated on key1
+  !
+  ! idx_repeat(0) is the number of determinants that can be used 
+  ! 
+  ! to repeat the excitations 
+  END_DOC
+  implicit none
+  integer, intent(in)            :: Nint, Ndet,Ndet_max,Nstate
+  integer(bit_kind), intent(in)  :: keys(Nint,2,Ndet)
+  integer(bit_kind), intent(in)  :: key(Nint,2)
+  double precision, intent(in)   :: coef(Ndet_max,Nstate)
+  double precision, intent(out)  :: i_H_psi_array(Nstate)
+  integer         , intent(out)  :: idx_repeat(0:Ndet)
+  
+  integer                        :: i, ii,j
+  double precision               :: phase
+  integer                        :: exc(0:2,2,2)
+  double precision               :: hij
+  integer                        :: idx(0:Ndet)
+  
+  ASSERT (Nint > 0)
+  ASSERT (N_int == Nint)
+  ASSERT (Nstate > 0)
+  ASSERT (Ndet > 0)
+  ASSERT (Ndet_max >= Ndet)
+  i_H_psi_array = 0.d0
+  call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat)
+  do ii=1,idx(0)
+    i = idx(ii)
+    !DEC$ FORCEINLINE
+    call i_H_j(keys(1,1,i),key,Nint,hij)
+    do j = 1, Nstate
+      i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij
+    enddo
+  enddo
+end
+
+
 
 subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx)
   use bitmasks

src/Makefile.config.example

@@ -0,0 +1,29 @@
+OPENMP  =1
+PROFILE =0
+DEBUG  = 0
+
+IRPF90_FLAGS+= --align=32
+FC     = ifort -g 
+FCFLAGS= 
+FCFLAGS+= -xHost
+#FCFLAGS+= -xAVX
+FCFLAGS+= -O2 
+FCFLAGS+= -ip 
+FCFLAGS+= -opt-prefetch
+FCFLAGS+= -ftz
+MKL=-mkl=parallel 
+
+ifeq ($(PROFILE),1)
+FC    += -p -g
+CXX   += -pg
+endif
+
+ifeq ($(OPENMP),1)
+FC    += -openmp
+CXX   += -fopenmp
+endif
+
+ifeq ($(DEBUG),1)
+FC    += -C -traceback -fpe0 
+#FCFLAGS =-O0
+endif

src/Perturbation/Moller_plesset.irp.f

@@ -0,0 +1,41 @@
+subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,n_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
+  double precision               :: i_H_psi_array(N_st)
+  
+  BEGIN_DOC
+  ! compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various n_st states.
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies) 
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies) 
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem
+  integer                        :: exc(0:2,2,2)
+  integer                        :: degree
+  double precision               :: phase,delta_e
+  integer                        :: h1,h2,p1,p2,s1,s2
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call get_excitation(det_pert,ref_bitmask,exc,degree,phase,Nint)
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+  delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
+            Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2) 
+  delta_e = 1.d0/delta_e
+
+  call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
+  H_pert_diag = diag_H_mat_elem(det_pert,Nint)
+  do i =1,n_st
+    c_pert(i) = i_H_psi_array(i) *delta_e
+    e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+  enddo
+  
+end

src/Perturbation/README.rst

@@ -82,6 +82,9 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/cisd_test.irp.f#L1>`_
+  Undocumented
+
 `pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
   compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
@@ -102,6 +105,114 @@ Documentation
   c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
   .br
 
+`pt2_epstein_nesbet_2x2_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L118>`_
+  compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
+  .br
+  for the various n_st states.
+  .br
+  but  with the correction in the denominator
+  .br
+  comming from the interaction of that determinant with all the others determinants
+  .br
+  that can be repeated by repeating all the double excitations
+  .br
+  : you repeat all the correlation energy already taken into account in reference_energy(1)
+  .br
+  that could be repeated to this determinant.
+  .br
+  <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
+  .br
+  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  )
+  .br
+  c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
+  .br
+
+`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L68>`_
+  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  .br
+  for the various n_st states,
+  .br
+  but  with the correction in the denominator
+  .br
+  comming from the interaction of that determinant with all the others determinants
+  .br
+  that can be repeated by repeating all the double excitations
+  .br
+  : you repeat all the correlation energy already taken into account in reference_energy(1)
+  .br
+  that could be repeated to this determinant.
+  .br
+  <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
+  .br
+  c_pert(i) = <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
+  .br
+  e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - (<det_pert|H|det_pert> ) )
+  .br
+
+`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L170>`_
+  compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  .br
+  for the various n_st states,
+  .br
+  but  with the correction in the denominator
+  .br
+  comming from the interaction of that determinant with all the others determinants
+  .br
+  that can be repeated by repeating all the double excitations
+  .br
+  : you repeat all the correlation energy already taken into account in reference_energy(1)
+  .br
+  that could be repeated to this determinant.
+  .br
+  BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection
+  .br
+  <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
+  .br
+  c_pert(1) = 1/c_HF <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
+  .br
+  e_2_pert(1) = <HF|H|det_pert> c_pert(1)
+  .br
+  NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !!
+  .br
+
+`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L1>`_
+  Fill the H_apply buffer with determiants for the selection
+
+`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L67>`_
+  Threshold to select determinants. Set by selection routines.
+
+`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L69>`_
+  Threshold to select determinants. Set by selection routines.
+
+`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L68>`_
+  Threshold to select determinants. Set by selection routines.
+
+`diagonalize <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L18>`_
+  Undocumented
+
+`h_apply_cisd_pt2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L63>`_
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry).
+
+`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L97>`_
+  Calls H_apply on the HF determinant and selects all connected single and double
+  excitations (of the same symmetry).
+
+`n_det_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L36>`_
+  Undocumented
+
+`psi_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L41>`_
+  On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
+
+`psi_ref_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L42>`_
+  On what we apply <i|H|psi> for perturbation. If selection, it may be 0.9 of the norm.
+
+`psi_ref_size <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L32>`_
+  Undocumented
+
+`reference_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L6>`_
+  Reference energy
+
 
 
 Needed Modules

src/Perturbation/epstein_nesbet.irp.f

@@ -64,3 +64,164 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
   enddo
 
 end
+
+subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,n_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: idx_repeat(ndet)
+  
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various n_st states, 
+  ! 
+  ! but  with the correction in the denominator 
+  ! 
+  ! comming from the interaction of that determinant with all the others determinants 
+  ! 
+  ! that can be repeated by repeating all the double excitations
+  !
+  ! : you repeat all the correlation energy already taken into account in reference_energy(1)
+  ! 
+  ! that could be repeated to this determinant.
+  !
+  ! <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
+  !
+  ! c_pert(i) = <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
+  !
+  ! e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - (<det_pert|H|det_pert> ) )
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem,accu_e_corr,hij
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat)
+  accu_e_corr = 0.d0
+  do i = 1, idx_repeat(0)
+   call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij)
+   accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i),1)
+  enddo
+  accu_e_corr = accu_e_corr / psi_ref_coef(1,1)
+  H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+  do i =1,n_st
+    c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag)
+    e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+  enddo
+  
+end
+subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,n_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: idx_repeat(ndet)
+
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
+  !
+  ! for the various n_st states.
+  ! 
+  ! but  with the correction in the denominator 
+  ! 
+  ! comming from the interaction of that determinant with all the others determinants 
+  ! 
+  ! that can be repeated by repeating all the double excitations
+  !
+  ! : you repeat all the correlation energy already taken into account in reference_energy(1)
+  ! 
+  ! that could be repeated to this determinant.
+  !
+  ! <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
+  !
+  ! 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,accu_e_corr,hij,delta_e
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat)
+  accu_e_corr = 0.d0
+  do i = 1, idx_repeat(0)
+   call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij)
+   accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i),1)
+  enddo
+  accu_e_corr = accu_e_corr / psi_ref_coef(1,1)
+  H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+  do i =1,n_st
+    delta_e = H_pert_diag - reference_energy(i)
+    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)))
+    c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
+  enddo
+  
+end
+
+subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: Nint,ndet,n_st
+  integer(bit_kind), intent(in)  :: det_pert(Nint,2)
+  double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
+  double precision               :: i_H_psi_array(N_st)
+  integer                        :: idx_repeat(ndet)
+  
+  BEGIN_DOC
+  ! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
+  !
+  ! for the various n_st states, 
+  ! 
+  ! but  with the correction in the denominator 
+  ! 
+  ! comming from the interaction of that determinant with all the others determinants 
+  ! 
+  ! that can be repeated by repeating all the double excitations
+  !
+  ! : you repeat all the correlation energy already taken into account in reference_energy(1)
+  ! 
+  ! that could be repeated to this determinant.
+  !
+  ! BUT on the contrary with ""pt2_epstein_nesbet_SC2"", you compute the energy by projection 
+  !
+  ! <det_pert|H|det_pert> --->  <det_pert|H|det_pert> + delta_e_corr
+  !
+  ! c_pert(1) = 1/c_HF <psi(i)|H|det_pert>/( E(i) - (<det_pert|H|det_pert> ) )
+  !
+  ! e_2_pert(1) = <HF|H|det_pert> c_pert(1)
+  !
+  ! NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !!
+  !
+  END_DOC
+  
+  integer                        :: i,j
+  double precision               :: diag_H_mat_elem,accu_e_corr,hij,h0j
+  ASSERT (Nint == N_int)
+  ASSERT (Nint > 0)
+  call i_H_psi_SC2(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array,idx_repeat)
+  accu_e_corr = 0.d0
+  call i_H_j(ref_bitmask,det_pert,Nint,h0j)
+  do i = 1, idx_repeat(0)
+   call i_H_j(psi_ref(1,1,idx_repeat(i)),det_pert,Nint,hij)
+   accu_e_corr = accu_e_corr + hij * psi_ref_coef(idx_repeat(i),1)
+  enddo
+  accu_e_corr = accu_e_corr / psi_ref_coef(1,1)
+  H_pert_diag = diag_H_mat_elem(det_pert,Nint) + accu_e_corr
+
+  c_pert(1) = 1.d0/psi_ref_coef(1,1) * i_H_psi_array(1) / (reference_energy(i) - H_pert_diag)
+  e_2_pert(1) = c_pert(i) * h0j
+  do i =2,n_st
+    c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag)
+    e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
+  enddo
+  
+end

src/SingleRefMethod/README.rst

@@ -11,6 +11,14 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`n_det_generators <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod/generators.irp.f#L3>`_
+  For Single reference wave functions, the number of generators is 1 : the
+  Hartree-Fock determinant
+
+`psi_generators <http://github.com/LCPQ/quantum_package/tree/master/src/SingleRefMethod/generators.irp.f#L12>`_
+  For Single reference wave functions, the generator is the
+  Hartree-Fock determinant
+
 
 
 Needed Modules

src/Utils/README.rst

@@ -10,20 +10,41 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`apply_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L146>`_
+`apply_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L169>`_
   Apply the rotation found by find_rotation
 
-`find_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L127>`_
+`find_rotation <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L150>`_
   Find A.C = B
 
-`get_pseudo_inverse <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L73>`_
+`get_pseudo_inverse <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L96>`_
   Find C = A^-1
 
-`lapack_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L158>`_
+`lapack_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L181>`_
   Diagonalize matrix H
+  .br
+  H is untouched between input and ouptut
+  .br
+  eigevalues(i) = ith lowest eigenvalue of the H matrix
+  .br
+  eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
+  .br
 
 `ortho_lowdin <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/LinearAlgebra.irp.f#L1>`_
-  Compute U.S^-1/2 canonical orthogonalization
+  Compute C_new=C_old.S^-1/2 canonical orthogonalization.
+  .br
+  overlap : overlap matrix
+  .br
+  LDA : leftmost dimension of overlap array
+  .br
+  N : Overlap matrix is NxN (array is (LDA,N) )
+  .br
+  C : Coefficients of the vectors to orthogonalize. On exit,
+  orthogonal vectors
+  .br
+  LDC : leftmost dimension of C
+  .br
+  m : Coefficients matrix is MxN, ( array is (LDC,N) )
+  .br
 
 `add_poly <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/integration.irp.f#L243>`_
   Add two polynomials
@@ -133,20 +154,18 @@ Documentation
 `inv_int <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L171>`_
   1/i
 
-`normalize <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L272>`_
+`normalize <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L270>`_
   Normalizes vector u
   u is expected to be aligned in memory.
 
 `nproc <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L197>`_
   Number of current OpenMP threads
 
-`u_dot_u <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L244>`_
+`u_dot_u <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L239>`_
   Compute <u|u>
-  u is expected to be aligned in memory.
 
 `u_dot_v <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L213>`_
   Compute <u|v>
-  u and v are expected to be aligned in memory.
 
 `wall_time <http://github.com/LCPQ/quantum_package/tree/master/src/Utils/util.irp.f#L182>`_
   The equivalent of cpu_time, but for the wall time.

src/Utils/util.irp.f

@@ -174,7 +174,7 @@ BEGIN_PROVIDER [ double precision, inv_int, (128) ]
   ! 1/i
   END_DOC
   integer                        :: i
-  do i=1,size(inv_int)
+  do i=1,128
     inv_int(i) = 1.d0/dble(i)
   enddo
 END_PROVIDER