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

Conflicts:
	src/MOGuess/NEEDED_MODULES

scripts/generate_h_apply.py

@@ -27,7 +27,7 @@ filter_integrals
 
 class H_apply(object):
 
-  def __init__(self,sub,SingleRef=False):
+  def __init__(self,sub,SingleRef=False,do_mono_exc=True, do_double_exc=True):
     s = {}
     for k in keywords:
       s[k] = ""
@@ -56,6 +56,9 @@ class H_apply(object):
     s["omp_do"]           = "!$OMP DO SCHEDULE (static)"
     s["omp_enddo"]        = "!$OMP ENDDO NOWAIT"
 
+    d = { True : '.True.', False : '.False.'}
+    s["do_mono_excitations"] = d[do_mono_exc]
+    s["do_double_excitations"] = d[do_double_exc]
     s["keys_work"]  += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)"
 
     s["filter_integrals"] = "array_pairs = .True."

src/AOs/aos.irp.f

@@ -93,46 +93,6 @@ END_PROVIDER
  enddo
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, ao_overlap, (ao_num_align,ao_num) ]
- implicit none
- BEGIN_DOC
-! matrix of the overlap for tha AOs 
-! S(i,j) = overlap between the ith and the jth atomic basis function
- END_DOC
- integer :: i,j,k,l,nz,num_i,num_j,powA(3),powB(3)
- double precision :: accu,overlap_x,overlap_y,overlap_z,A_center(3),B_center(3),norm
- nz=100
- do i = 1, ao_num 
-  num_i = ao_nucl(i)
-  powA(1) = ao_power(i,1)
-  powA(2) = ao_power(i,2)
-  powA(3) = ao_power(i,3)
-  A_center(1)=nucl_coord(num_i,1)
-  A_center(2)=nucl_coord(num_i,2)
-  A_center(3)=nucl_coord(num_i,3)
-  do j = i, ao_num
-   num_j = ao_nucl(j)
-   powB(1) = ao_power(j,1)
-   powB(2) = ao_power(j,2)
-   powB(3) = ao_power(j,3)
-   B_center(1)=nucl_coord(num_j,1)
-   B_center(2)=nucl_coord(num_j,2)
-   B_center(3)=nucl_coord(num_j,3)
-   accu = 0.d0
-   do k = 1, ao_prim_num(i)
-    do l = 1, ao_prim_num(j)
-     call overlap_gaussian_xyz(A_center,B_center,ao_expo(i,k),ao_expo(j,l),powA,powB,overlap_x,overlap_y,overlap_z,norm,nz)
-     accu = accu + norm * ao_coef(i,k) * ao_coef(j,l)
-    enddo
-   enddo
-   ao_overlap(i,j) = accu
-   ao_overlap(j,i) = accu
-  enddo
- enddo
-
-END_PROVIDER
-
-
  BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ]
 &BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ]
  implicit none

src/CID/H_apply.irp.f

@@ -3,7 +3,7 @@
 
 BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import H_apply
-H = H_apply("cisd")
+H = H_apply("cisd",do_double_exc=True,do_mono_exc=False)
 print H
 END_SHELL
 

src/CID/README.rst

@@ -36,7 +36,7 @@ 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/CISD/cisd_lapack.irp.f#L1>`_
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CID/cid_lapack.irp.f#L1>`_
   Undocumented
 
 

src/CID/do_mono_double.irp.f

@@ -1,19 +0,0 @@
-BEGIN_PROVIDER [logical, do_double_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the double excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_double_excitations = .True.
-
-END_PROVIDER
-
-BEGIN_PROVIDER [logical, do_mono_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the mono excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_mono_excitations = .False.
-
-END_PROVIDER

src/CID_SC2_selected/H_apply.irp.f

@@ -3,7 +3,7 @@ BEGIN_SHELL [ /usr/bin/env python ]
 from generate_h_apply import *
 from perturbation import perturbations
 
-s = H_apply("PT2",SingleRef=True)
+s = H_apply("PT2",SingleRef=True,do_mono_exc=False,do_double_exc=True)
 s.set_perturbation("epstein_nesbet_sc2_projected")
 print s
 END_SHELL

src/CID_SC2_selected/README.rst

@@ -8,7 +8,7 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`cisd_sc2_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2_selected/cisd_sc2_selection.irp.f#L1>`_
+`cisd_sc2_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CID_SC2_selected/cid_sc2_selection.irp.f#L1>`_
   Undocumented
 
 

src/CID_SC2_selected/do_mono_double.irp.f

@@ -1,19 +0,0 @@
-BEGIN_PROVIDER [logical, do_double_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the double excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_double_excitations = .True.
-
-END_PROVIDER
-
-BEGIN_PROVIDER [logical, do_mono_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the mono excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_mono_excitations = .False.
-
-END_PROVIDER

src/CID_selected/H_apply.irp.f

@@ -4,7 +4,7 @@ from generate_h_apply import *
 from perturbation import perturbations
 
 for perturbation in perturbations:
-  s = H_apply("cisd_selection_"+perturbation)
+  s = H_apply("cisd_selection_"+perturbation,do_mono_exc=False)
   s.set_selection_pt2(perturbation)
   print s
 END_SHELL

src/CID_selected/README.rst

@@ -8,10 +8,10 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/H_apply.irp.f#L13>`_
+`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/CID_selected/H_apply.irp.f#L13>`_
   Undocumented
 
-`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected/cisd_selection.irp.f#L1>`_
+`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CID_selected/cid_selection.irp.f#L1>`_
   Undocumented
 
 

src/CID_selected/do_mono_double.irp.f

@@ -1,19 +0,0 @@
-BEGIN_PROVIDER [logical, do_double_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the double excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_double_excitations = .True.
-
-END_PROVIDER
-
-BEGIN_PROVIDER [logical, do_mono_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the mono excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_mono_excitations = .False.
-
-END_PROVIDER

src/CISD/do_mono_double.irp.f

@@ -1,19 +0,0 @@
-BEGIN_PROVIDER [logical, do_double_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the double excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_double_excitations = .True.
-
-END_PROVIDER
-
-BEGIN_PROVIDER [logical, do_mono_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the mono excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_mono_excitations = .True.
-
-END_PROVIDER

src/CISD_SC2_selected/do_mono_double.irp.f

@@ -1,19 +0,0 @@
-BEGIN_PROVIDER [logical, do_double_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the double excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_double_excitations = .True.
-
-END_PROVIDER
-
-BEGIN_PROVIDER [logical, do_mono_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the mono excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_mono_excitations = .True.
-
-END_PROVIDER

src/CISD_selected/do_mono_double.irp.f

@@ -1,19 +0,0 @@
-BEGIN_PROVIDER [logical, do_double_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the double excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_double_excitations = .True.
-
-END_PROVIDER
-
-BEGIN_PROVIDER [logical, do_mono_excitations]
- implicit none
- BEGIN_DOC
- ! if True then the mono excitations are performed in the calculation
- ! always true in the CISD
- END_DOC
- do_mono_excitations = .True.
-
-END_PROVIDER

src/DensityMatrix/density_matrix.irp.f

@@ -209,51 +209,3 @@ END_PROVIDER
      enddo
    enddo
 END_PROVIDER
- 
- BEGIN_PROVIDER [ double precision, one_body_dm_a, (mo_tot_num_align,mo_tot_num) ]
-&BEGIN_PROVIDER [ double precision, one_body_dm_b, (mo_tot_num_align,mo_tot_num) ]
-   implicit none
-   BEGIN_DOC
-   ! Alpha and beta one-body density matrix
-   END_DOC
-   
-   integer                        :: j,k,l
-   integer                        :: occ(N_int*bit_kind_size,2)
-   double precision               :: ck, cl, ckl
-   double precision               :: phase
-   integer                        :: h1,h2,p1,p2,s1,s2, degree
-   integer                        :: exc(0:2,2,2),n_occ_alpha
-   one_body_dm_a = 0.d0
-   one_body_dm_b = 0.d0
-
-   do k=1,det_num
-     call bitstring_to_list(det_provider(1,1,k), occ(1,1), n_occ_alpha, N_int)
-     call bitstring_to_list(det_provider(1,2,k), occ(1,2), n_occ_alpha, N_int)
-     ck = det_coef_provider(k)
-     do l=1,elec_alpha_num
-       j = occ(l,1)
-       one_body_dm_a(j,j) += ck*ck
-     enddo
-     do l=1,elec_beta_num
-       j = occ(l,2)
-       one_body_dm_b(j,j) += ck*ck
-     enddo
-     do l=1,k-1
-       call get_excitation_degree(det_provider(1,1,k),det_provider(1,1,l),degree,N_int)
-       if (degree /= 1) then
-         cycle
-       endif
-       call get_mono_excitation(det_provider(1,1,k),det_provider(1,1,l),exc,phase,N_int)
-       call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-       ckl = ck * det_coef_provider(l) * phase
-       if (s1==1) then
-         one_body_dm_a(h1,p1) += ckl
-         one_body_dm_a(p1,h1) += ckl
-       else
-         one_body_dm_b(h1,p1) += ckl
-         one_body_dm_b(p1,h1) += ckl
-       endif
-     enddo
-   enddo
-END_PROVIDER
- 

src/Dets/H_apply_template.f

@@ -389,7 +389,7 @@ subroutine $subroutine($params_main)
   !$OMP PARALLEL DEFAULT(SHARED) &
   !$OMP PRIVATE(i_generator,wall_2,ispin,k,mask) 
   allocate( mask(N_int,2,6) )
-  !$OMP DO SCHEDULE(guided)
+  !$OMP DO SCHEDULE(dynamic,4)
   do i_generator=1,nmax
     if (abort_here) then
       cycle
@@ -420,13 +420,13 @@ subroutine $subroutine($params_main)
       enddo
     enddo
 
-    if(do_double_excitations)then
+    if($do_double_excitations)then
     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
         mask(1,1,d_hole1), mask(1,1,d_part1),                        &
         mask(1,1,d_hole2), mask(1,1,d_part2),                        &
         i_generator $params_post)
     endif
-    if(do_mono_excitations)then
+    if($do_mono_excitations)then
     call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
         mask(1,1,s_hole ), mask(1,1,s_part ),                        &
         i_generator $params_post)
@@ -475,13 +475,13 @@ subroutine $subroutine($params_main)
             not(psi_generators(k,ispin,i_generator)) )
       enddo
     enddo
-    if(do_double_excitations)then
+    if($do_double_excitations)then
      call $subroutine_diexc(psi_generators(1,1,i_generator),          &
          mask(1,1,d_hole1), mask(1,1,d_part1),                        &
          mask(1,1,d_hole2), mask(1,1,d_part2),                        &
          i_generator $params_post)
     endif
-    if(do_mono_excitations)then
+    if($do_mono_excitations)then
      call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
          mask(1,1,s_hole ), mask(1,1,s_part ),                        &
          i_generator $params_post)

src/Dets/README.rst

@@ -50,11 +50,11 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L113>`_
+`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L112>`_
   Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det
   after calling this function.
 
-`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L199>`_
+`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L198>`_
   Fill the H_apply buffer with determiants for CISD
 
 `h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L14>`_
@@ -64,7 +64,7 @@ Documentation
 `h_apply_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L44>`_
   Theshold on | <Di|H|Dj> |
 
-`resize_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L63>`_
+`resize_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L62>`_
   Undocumented
 
 `cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/SC2.irp.f#L1>`_
@@ -152,55 +152,55 @@ Documentation
 `davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/davidson.irp.f#L374>`_
   Can be : [  energy  | residual | both | wall_time | cpu_time | iterations ]
 
-`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L273>`_
+`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L303>`_
   Return an integer*8 corresponding to a determinant index for searching
 
 `n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L20>`_
   Number of determinants in the wave function
 
-`n_det_max_jacobi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L38>`_
+`n_det_max_jacobi <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L48>`_
   Maximum number of determinants diagonalized my jacobi
 
 `n_states <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L3>`_
   Number of states to consider
 
-`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L186>`_
+`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L216>`_
   Contribution of determinants to the state-averaged density
 
-`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L207>`_
+`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L237>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L150>`_
+`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L171>`_
   The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L206>`_
+`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L236>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L237>`_
+`psi_coef_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L267>`_
   Determinants on which we apply <i|H|psi> for perturbation.
   o They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a determinant
 
-`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L64>`_
+`psi_det <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L84>`_
   The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
   is empty
 
-`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L56>`_
+`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L66>`_
   Size of the psi_det/psi_coef arrays
 
-`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L205>`_
+`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L235>`_
   Wave function sorted by determinants contribution to the norm (state-averaged)
 
-`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L236>`_
+`psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L266>`_
   Determinants on which we apply <i|H|psi> for perturbation.
   o They are sorted by determinants interpreted as integers. Useful
   to accelerate the search of a determinant
 
-`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L107>`_
+`read_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L127>`_
   Reads the determinants from the EZFIO file
 
-`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L288>`_
+`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants.irp.f#L318>`_
   Save the wave function into the EZFIO file
 
 `double_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/determinants_bitmasks.irp.f#L40>`_
@@ -237,19 +237,22 @@ Documentation
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
-`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L19>`_
+`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L27>`_
   Eigenvectors/values of the CI matrix
 
-`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L18>`_
+`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L26>`_
   Eigenvectors/values of the CI matrix
 
 `ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L1>`_
   N_states lowest eigenvalues of the CI matrix
 
-`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L38>`_
+`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L46>`_
   Replace the coefficients of the CI states by the coefficients of the
   eigenstates of the CI matrix
 
+`threshold_convergence_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/diagonalize_CI_SC2.irp.f#L18>`_
+  convergence of the correlation energy of SC2 iterations
+
 `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

src/Dets/density_matrix.irp.f

@@ -0,0 +1,106 @@
+ BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ]
+&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Alpha and beta one-body density matrix for each state
+   END_DOC
+
+   integer                        :: j,k,l,m
+   integer                        :: occ(N_int*bit_kind_size,2)
+   double precision               :: ck, cl, ckl
+   double precision               :: phase
+   integer                        :: h1,h2,p1,p2,s1,s2, degree
+   integer                        :: exc(0:2,2,2),n_occ_alpha
+   double precision, allocatable  :: tmp_a(:,:), tmp_b(:,:)
+
+   one_body_dm_mo_alpha = 0.d0
+   one_body_dm_mo_beta  = 0.d0
+   !$OMP PARALLEL DEFAULT(NONE)                                         &
+      !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, &
+      !$OMP  tmp_a, tmp_b, n_occ_alpha)&
+      !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,&
+      !$OMP  elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,&
+      !$OMP  mo_tot_num)
+   allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) )
+   tmp_a = 0.d0
+   tmp_b = 0.d0
+   !$OMP DO SCHEDULE(dynamic)
+   do k=1,N_det
+     call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int)
+     call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int)
+     do m=1,N_states
+       ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m)
+       do l=1,elec_alpha_num
+         j = occ(l,1)
+         tmp_a(j,j) += ck
+       enddo
+       do l=1,elec_beta_num
+         j = occ(l,2)
+         tmp_b(j,j) += ck
+       enddo
+     enddo
+     do l=1,k-1
+       call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int)
+       if (degree /= 1) then
+         cycle
+       endif
+       call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int)
+       call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+       do m=1,N_states
+         ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m)
+         if (s1==1) then
+           tmp_a(h1,p1) += ckl
+           tmp_a(p1,h1) += ckl
+         else
+           tmp_b(h1,p1) += ckl
+           tmp_b(p1,h1) += ckl
+         endif
+       enddo
+     enddo
+   enddo
+   !$OMP END DO NOWAIT
+   !$OMP CRITICAL
+   one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a
+   !$OMP END CRITICAL
+   !$OMP CRITICAL
+   one_body_dm_mo_beta  = one_body_dm_mo_beta  + tmp_b
+   !$OMP END CRITICAL
+   deallocate(tmp_a,tmp_b)
+   !$OMP BARRIER
+   !$OMP END PARALLEL
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! One-body density matrix
+ END_DOC
+ one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta
+END_PROVIDER
+
+subroutine save_natural_mos
+ implicit none
+ BEGIN_DOC
+ ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
+ END_DOC
+ character*(64) :: label
+ double precision, allocatable :: tmp(:,:)
+ allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2)))
+
+ ! Negation to have the occupied MOs first after the diagonalization
+ tmp = -one_body_dm_mo
+ label = "Natural"
+ call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label)
+ deallocate(tmp)
+ call save_mos
+
+end
+
+BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Weights in the state-average calculation of the density matrix
+ END_DOC
+ state_average_weight = 1.d0/dble(N_states)
+END_PROVIDER
+

src/Dets/save_natorb.irp.f

@@ -0,0 +1,4 @@
+program save_natorb
+  call save_natural_mos
+end
+

src/Full_CI/full_ci.irp.f

@@ -25,14 +25,14 @@ program cisd
       exit
     endif
   enddo
-  print *,  'Final step'
+!  print *,  'Final step'
-  call remove_small_contributions
+!  call remove_small_contributions
-  call diagonalize_CI
+!  call diagonalize_CI
-  print *,  'N_det    = ', N_det
+!  print *,  'N_det    = ', N_det
-  print *,  'N_states = ', N_states
+!  print *,  'N_states = ', N_states
-  print *,  'PT2      = ', pt2
+!  print *,  'PT2      = ', pt2
-  print *,  'E        = ', CI_energy
+!  print *,  'E        = ', CI_energy
-  print *,  'E+PT2    = ', CI_energy+pt2
+!  print *,  'E+PT2    = ', CI_energy+pt2
-  print *,  '-----'
+!  print *,  '-----'
   deallocate(pt2,norm_pert)
 end

src/Generators_full/generators.irp.f

@@ -47,7 +47,13 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
  ! For Single reference wave functions, the generator is the
  ! Hartree-Fock determinant
  END_DOC
- psi_generators = psi_det_sorted
+ integer                        :: i, k
+ do i=1,N_det
+   do k=1,N_int
+     psi_generators(k,1,i) = psi_det_sorted(k,1,i)
+     psi_generators(k,2,i) = psi_det_sorted(k,2,i)
+   enddo
+ enddo
 
 END_PROVIDER
 

src/Hartree_Fock/Fock_matrix.irp.f

@@ -217,3 +217,46 @@ BEGIN_PROVIDER [ double precision, HF_energy ]
 END_PROVIDER
 
 
+BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Fock matrix in AO basis set
+ END_DOC
+ 
+ if (elec_alpha_num == elec_beta_num) then
+ integer :: i,j
+   do j=1,ao_num
+    !DIR$ VECTOR ALIGNED
+    do i=1,ao_num_align
+     Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j)
+    enddo
+   enddo
+ else
+   double precision, allocatable :: T(:,:), M(:,:)
+   ! F_ao = S C F_mo C^t S
+   allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num))
+   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
+       ao_overlap, size(ao_overlap,1),                               &
+       mo_coef, size(mo_coef,1),                                     &
+       0.d0,                                                         &
+       M, size(M,1))
+   call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,           &
+       M, size(M,1),                                                 &
+       Fock_matrix_mo, size(Fock_matrix_mo,1),                       &
+       0.d0,                                                         &
+       T, size(T,1))
+   call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,           &
+       T, size(T,1),                                                 &
+       mo_coef, size(mo_coef,1),                                     &
+       0.d0,                                                         &
+       M, size(M,1))
+   call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                   &
+       M, size(M,1),                                                 &
+       ao_overlap, size(ao_overlap,1),                               &
+       0.d0,                                                         &
+       Fock_matrix_ao, size(Fock_matrix_ao,1))
+
+   deallocate(T)
+ endif
+END_PROVIDER
+

src/Hartree_Fock/HF_density_matrix_ao.irp.f

@@ -1,46 +1,27 @@
- BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
+BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_alpha, (ao_num_align,ao_num) ]
-&BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]
    implicit none
    BEGIN_DOC
-   ! Alpha and Beta density matrix in the AO basis
+   ! Alpha density matrix in the AO basis
    END_DOC
-   integer                        :: i,j,k,l1,l2
-   integer, allocatable           :: mo_occ(:,:)
    
-   allocate ( mo_occ(elec_alpha_num,2) )
+   call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
-   call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
+        mo_coef, size(mo_coef,1), &
-   ASSERT ( j==elec_alpha_num )
+        mo_coef, size(mo_coef,1), 0.d0, &
+        HF_density_matrix_ao_alpha, size(HF_density_matrix_ao_alpha,1))
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta,  (ao_num_align,ao_num) ]
+   implicit none
+   BEGIN_DOC
+   ! Beta density matrix in the AO basis
+   END_DOC
    
-   call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
+   call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
-   ASSERT ( j==elec_beta_num )
+        mo_coef, size(mo_coef,1), &
-   
+        mo_coef, size(mo_coef,1), 0.d0, &
-   do j=1,ao_num
+        HF_density_matrix_ao_beta, size(HF_density_matrix_ao_beta,1))
-     !DIR$ VECTOR ALIGNED
+
-     do i=1,ao_num_align
-       HF_density_matrix_ao_alpha(i,j) = 0.d0
-       HF_density_matrix_ao_beta (i,j) = 0.d0
-     enddo
-     do k=1,elec_beta_num
-       l1 = mo_occ(k,1)
-       l2 = mo_occ(k,2)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
-             mo_coef(i,l1) * mo_coef(j,l1)
-         HF_density_matrix_ao_beta (i,j) = HF_density_matrix_ao_beta (i,j) +&
-             mo_coef(i,l2) * mo_coef(j,l2)
-       enddo
-     enddo
-     do k=elec_beta_num+1,elec_alpha_num
-       l1 = mo_occ(k,1)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao_alpha(i,j) = HF_density_matrix_ao_alpha(i,j) +&
-             mo_coef(i,l1) * mo_coef(j,l1)
-       enddo
-     enddo
-   enddo
-   deallocate(mo_occ)
 END_PROVIDER
  
 BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
@@ -48,40 +29,13 @@ BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num_align,ao_num) ]
    BEGIN_DOC
    ! Density matrix in the AO basis
    END_DOC
-   integer                        :: i,j,k,l1,l2
+   ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
-   integer, allocatable           :: mo_occ(:,:)
+   if (elec_alpha_num== elec_beta_num) then
+     HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_alpha
+   else
+     ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_beta ,1))
+     HF_density_matrix_ao = HF_density_matrix_ao_alpha + HF_density_matrix_ao_beta
+   endif
    
-   allocate ( mo_occ(elec_alpha_num,2) )
-   call bitstring_to_list( HF_bitmask(1,1), mo_occ(1,1), j, N_int)
-   ASSERT ( j==elec_alpha_num )
-   
-   call bitstring_to_list( HF_bitmask(1,2), mo_occ(1,2), j, N_int)
-   ASSERT ( j==elec_beta_num )
-   
-   do j=1,ao_num
-     !DIR$ VECTOR ALIGNED
-     do i=1,ao_num_align
-       HF_density_matrix_ao(i,j) = 0.d0
-     enddo
-     do k=1,elec_beta_num
-       l1 = mo_occ(k,1)
-       l2 = mo_occ(k,2)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) +           &
-             mo_coef(i,l1) * mo_coef(j,l1) +                         &
-             mo_coef(i,l2) * mo_coef(j,l2)
-       enddo
-     enddo
-     do k=elec_beta_num+1,elec_alpha_num
-       l1 = mo_occ(k,1)
-       !DIR$ VECTOR ALIGNED
-       do i=1,ao_num
-         HF_density_matrix_ao(i,j) = HF_density_matrix_ao(i,j) +           &
-             mo_coef(i,l1) * mo_coef(j,l1)
-       enddo
-     enddo
-   enddo
-   deallocate(mo_occ)
 END_PROVIDER
  

src/Hartree_Fock/README.rst

@@ -32,6 +32,9 @@ Documentation
 `fock_matrix_alpha_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L172>`_
   Fock matrix on the MO basis
 
+`fock_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L220>`_
+  Fock matrix in AO basis set
+
 `fock_matrix_beta_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L84>`_
   Alpha Fock matrix in AO basis set
 
@@ -71,14 +74,35 @@ Documentation
 `hf_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/Fock_matrix.irp.f#L211>`_
   Hartree-Fock energy
 
-`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L46>`_
+`hf_density_matrix_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L27>`_
   Density matrix in the AO basis
 
 `hf_density_matrix_ao_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L1>`_
-  Alpha and Beta density matrix in the AO basis
+  Alpha density matrix in the AO basis
 
-`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L2>`_
+`hf_density_matrix_ao_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/HF_density_matrix_ao.irp.f#L14>`_
-  Alpha and Beta density matrix in the AO basis
+  Beta density matrix in the AO basis
+
+`fock_mo_to_ao <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L31>`_
+  Undocumented
+
+`insert_new_scf_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L19>`_
+  Undocumented
+
+`it_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L1>`_
+  Number of the current SCF iteration
+
+`scf_density_matrices <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L9>`_
+  Density matrices at every SCF iteration
+
+`scf_energies <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L10>`_
+  Density matrices at every SCF iteration
+
+`scf_interpolation_step <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L64>`_
+  Undocumented
+
+`scf_iterations <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/SCF.irp.f#L89>`_
+  Undocumented
 
 `diagonal_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L1>`_
   Diagonal Fock matrix in the MO basis
@@ -86,13 +110,16 @@ Documentation
 `eigenvectors_fock_matrix_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/diagonalize_fock.irp.f#L2>`_
   Diagonal Fock matrix in the MO basis
 
-`scf_iteration <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L1>`_
+`run <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L7>`_
   Undocumented
 
-`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L41>`_
+`scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/mo_SCF_iterations.irp.f#L2>`_
+  Undocumented
+
+`do_diis <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L39>`_
   If True, compute integrals on the fly
 
-`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L22>`_
+`n_it_scf_max <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L21>`_
   Maximum number of SCF iterations
 
 `thresh_scf <http://github.com/LCPQ/quantum_package/tree/master/src/Hartree_Fock/options.irp.f#L1>`_

src/Hartree_Fock/SCF.irp.f

@@ -0,0 +1,108 @@
+BEGIN_PROVIDER [ integer, it_scf ]
+  implicit none
+  BEGIN_DOC
+  ! Number of the current SCF iteration
+  END_DOC
+  it_scf = 0
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, SCF_density_matrices, (ao_num_align,ao_num,2,n_it_scf_max) ]
+&BEGIN_PROVIDER [ double precision, SCF_energies, (n_it_scf_max) ]
+  implicit none
+  BEGIN_DOC
+  ! Density matrices at every SCF iteration
+  END_DOC
+  SCF_density_matrices = 0.d0
+  SCF_energies = 0.d0
+END_PROVIDER
+
+subroutine insert_new_SCF_density_matrix
+  implicit none
+  integer                        :: i,j
+  do j=1,ao_num
+    do i=1,ao_num
+      SCF_density_matrices(i,j,1,it_scf) = HF_density_matrix_ao_alpha(i,j)
+      SCF_density_matrices(i,j,2,it_scf) = HF_density_matrix_ao_beta(i,j)
+    enddo
+  enddo
+  SCF_energies(it_scf) = HF_energy
+end
+
+subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
+  implicit none
+  integer, intent(in)            :: LDFMO ! size(FMO,1)
+  integer, intent(in)            :: LDFAO ! size(FAO,1)
+  double precision, intent(in)   :: FMO(LDFMO,*)
+  double precision, intent(out)  :: FAO(LDFAO,*)
+  
+  double precision, allocatable  :: T(:,:), M(:,:)
+  ! F_ao = S C F_mo C^t S
+  allocate (T(mo_tot_num_align,mo_tot_num),M(ao_num_align,mo_tot_num))
+  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+      ao_overlap, size(ao_overlap,1),                                &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0,                                                          &
+      M, size(M,1))
+  call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0,            &
+      M, size(M,1),                                                  &
+      FMO, size(FMO,1),                                              &
+      0.d0,                                                          &
+      T, size(T,1))
+  call dgemm('N','T', mo_tot_num,ao_num,mo_tot_num, 1.d0,            &
+      T, size(T,1),                                                  &
+      mo_coef, size(mo_coef,1),                                      &
+      0.d0,                                                          &
+      M, size(M,1))
+  call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0,                    &
+      M, size(M,1),                                                  &
+      ao_overlap, size(ao_overlap,1),                                &
+      0.d0,                                                          &
+      FAO, size(FAO,1))
+  deallocate(T,M)
+end
+
+subroutine SCF_interpolation_step
+  implicit none
+  integer                        :: i,j
+  double precision               :: c
+  
+  if (it_scf == 1) then
+    return
+  endif
+  call random_number(c)
+  c = c*0.5d0
+  do j=1,ao_num
+    do i=1,ao_num
+      HF_density_matrix_ao_alpha(i,j) = c*SCF_density_matrices(i,j,1,it_scf-1)+SCF_density_matrices(i,j,1,it_scf) * (1.d0 - c)
+      HF_density_matrix_ao_beta (i,j) = c*SCF_density_matrices(i,j,2,it_scf-1)+SCF_density_matrices(i,j,2,it_scf) * (1.d0 - c)
+    enddo
+  enddo
+  TOUCH HF_density_matrix_ao_alpha HF_density_matrix_ao_beta
+  
+  ! call Fock_mo_to_ao(Fock_matrix_mo_alpha, size(Fock_matrix_mo_alpha,1),&
+      !      Fock_matrix_alpha_ao, size(Fock_matrix_alpha_ao,1) )
+  ! call Fock_mo_to_ao(Fock_matrix_mo_beta, size(Fock_matrix_mo_beta,1),&
+      !      Fock_matrix_beta_ao, size(Fock_matrix_beta_ao,1) )
+  ! SOFT_TOUCH Fock_matrix_alpha_ao Fock_matrix_beta_ao Fock_matrix_mo_alpha Fock_matrix_mo_beta
+end
+
+subroutine scf_iterations
+  implicit none
+  integer                        :: i,j
+  do i=1,n_it_scf_max
+    it_scf += 1
+    SOFT_TOUCH it_scf
+    mo_coef = eigenvectors_Fock_matrix_mo
+    TOUCH mo_coef
+    call insert_new_SCF_density_matrix
+    print *,  HF_energy
+    if (SCF_energies(it_scf)>SCF_energies(it_scf-1)) then
+      call SCF_interpolation_step
+    endif
+    if (it_scf>1 ) then
+      if (dabs(SCF_energies(it_scf)-SCF_energies(it_scf-1)) < thresh_SCF) then
+        exit
+      endif
+    endif
+  enddo
+end

src/Hartree_Fock/diagonalize_fock.irp.f

@@ -5,16 +5,52 @@
    ! Diagonal Fock matrix in the MO basis
    END_DOC
    
-   double precision, allocatable  :: R(:,:)
+   integer                        :: i,j
-   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: R
+   integer                        :: liwork, lwork, n, info
+   integer, allocatable           :: iwork(:)
+   double precision, allocatable  :: work(:), F(:,:), S(:,:)
    
-   allocate(R(mo_tot_num_align,mo_tot_num))
+   allocate(F(ao_num_align,ao_num), S(ao_num_align,ao_num) )
+   do j=1,ao_num
+     do i=1,ao_num
+       S(i,j) = ao_overlap(i,j)
+       F(i,j) = Fock_matrix_ao(i,j)
+     enddo
+   enddo
+
+   n = ao_num
+   lwork = 1+6*n + 2*n*n
+   liwork = 3 + 5*n
    
-   call lapack_diag(diagonal_Fock_matrix_mo,R,Fock_matrix_mo,size(Fock_matrix_mo,1),&
+   allocate(work(lwork), iwork(liwork) )
-       mo_tot_num)
+
-   
+   lwork = -1
-   call dgemm('N','N',ao_num,mo_tot_num,mo_tot_num,1.d0,mo_coef,size(mo_coef,1),&
+   liwork = -1
-       R,size(R,1),0.d0,eigenvectors_Fock_matrix_mo,size(eigenvectors_Fock_matrix_mo,1))
+   call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),&
-   deallocate(R)
+     diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+     print *,  irp_here//' failed'
+     stop 1
+   endif
+   lwork = int(work(1))
+   liwork = iwork(1)
+   deallocate(work,iwork)
+   allocate(work(lwork), iwork(liwork) )
+
+   call dsygvd(1,'v','u',mo_tot_num,F,size(F,1),S,size(S,1),&
+     diagonal_Fock_matrix_mo, work, lwork, iwork, liwork, info)
+
+   if (info /= 0) then
+     print *,  irp_here//' failed'
+     stop 1
+   endif
+   do j=1,ao_num
+     do i=1,ao_num
+       eigenvectors_Fock_matrix_mo(i,j) = F(i,j)
+     enddo
+   enddo
+
+   deallocate(work, iwork, F, S)
 END_PROVIDER
  

src/Hartree_Fock/mo_SCF_iterations.irp.f

@@ -1,35 +1,21 @@
-program scf_iteration
+
+program scf
+  call orthonormalize_mos
+  call run
+end
+
+subroutine run
+
   use bitmasks
   implicit none
   double precision               :: SCF_energy_before,SCF_energy_after,diag_H_mat_elem,get_mo_bielec_integral
   double precision               :: E0
-  integer                        :: i_it
+  integer                        :: i_it, i, j, k
-  
+   
   E0 = HF_energy 
-  i_it = 0
-  n_it_scf_max = 10
-  SCF_energy_before = huge(1.d0)
-  SCF_energy_after = E0
-  print *,  E0
-  mo_label = "Canonical"
-  thresh_SCF = 1.d-10
-  do while (i_it < 40 .and. dabs(SCF_energy_before - SCF_energy_after) > thresh_SCF)
-    SCF_energy_before = SCF_energy_after
-    mo_coef = eigenvectors_Fock_matrix_mo
-    TOUCH mo_coef mo_label
-    call clear_mo_map
-    SCF_energy_after = HF_energy
-    print*,SCF_energy_after, dabs(SCF_energy_before - SCF_energy_after) 
-    i_it +=1
-    if(i_it > n_it_scf_max)exit
-  enddo
   
-  if (i_it >= n_it_scf_max) then
+  thresh_SCF = 1.d-10
-    stop 'Failed'
+  call damping_SCF
-  endif
-  if (SCF_energy_after - E0 > thresh_SCF) then
-    stop 'Failed'
-  endif
   mo_label = "Canonical"
   TOUCH mo_label mo_coef
   call save_mos

src/MOGuess/README.rst

@@ -22,19 +22,19 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L1>`_
+`h_core_guess <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/H_CORE_guess.irp.f#L/subroutine h_core_guess/;">`_
   Undocumented
 
-`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L2>`_
+`ao_ortho_lowdin_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]/;">`_
   matrix of the coefficients of the mos generated by the
   orthonormalization by the S^{-1/2} canonical transformation of the aos
   ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital
 
-`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L26>`_
+`ao_ortho_lowdin_overlap <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/mo_ortho_lowdin.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)]/;">`_
   overlap matrix of the ao_ortho_lowdin
   supposed to be the Identity
 
-`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L1>`_
+`ao_ortho_lowdin_nucl_elec_integral <http://github.com/LCPQ/quantum_package/tree/master/src/MOGuess/pot_mo_ortho_lowdin_ints.irp.f#L/BEGIN_PROVIDER [double precision, ao_ortho_lowdin_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]/;">`_
   Undocumented
 
 

src/MOs/NEEDED_MODULES

@@ -1 +1 @@
-AOs Ezfio_files Nuclei Output Utils
+AOs Ezfio_files Nuclei Output Utils Electrons

src/MOs/README.rst

@@ -27,6 +27,7 @@ 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>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
 
 Documentation
 =============
@@ -34,12 +35,19 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
+`cholesky_mo <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L1>`_
+  Cholesky decomposition of MO Density matrix to
+  generate MOs
+
+`mo_density_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/cholesky_mo.irp.f#L44>`_
+  Density matrix in MO basis
+
 `mo_coef <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L22>`_
   Molecular orbital coefficients on AO basis set
   mo_coef(i,j) = coefficient of the ith ao on the jth mo
   mo_label : Label characterizing the MOS (local, canonical, natural, etc)
 
-`mo_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L61>`_
+`mo_coef_transp <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L60>`_
   Molecular orbital coefficients on AO basis set
 
 `mo_label <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L23>`_
@@ -47,13 +55,16 @@ Documentation
   mo_coef(i,j) = coefficient of the ith ao on the jth mo
   mo_label : Label characterizing the MOS (local, canonical, natural, etc)
 
+`mo_occ <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L78>`_
+  MO occupation numbers
+
 `mo_tot_num <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L1>`_
   Total number of molecular orbitals and the size of the keys corresponding
 
 `mo_tot_num_align <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/mos.irp.f#L12>`_
   Aligned variable for dimensioning of arrays
 
-`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L21>`_
+`mo_as_eigvectors_of_mo_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L22>`_
   Undocumented
 
 `save_mos <http://github.com/LCPQ/quantum_package/tree/master/src/MOs/utils.irp.f#L1>`_

src/MOs/cholesky_mo.irp.f

@@ -0,0 +1,80 @@
+subroutine cholesky_mo(n,m,P,LDP,C,LDC,tol_in,rank)
+ implicit none
+ BEGIN_DOC
+! Cholesky decomposition of AO Density matrix to
+! generate MOs
+ END_DOC
+ integer, intent(in) :: n,m, LDC, LDP
+ double precision, intent(in) :: P(LDP,n)
+ double precision, intent(out) :: C(LDC,m)
+ double precision, intent(in) :: tol_in
+ integer, intent(out) :: rank
+
+ integer :: info
+ integer :: i,k
+ integer :: ipiv(n)
+ double precision:: tol
+ double precision, allocatable :: W(:,:), work(:)
+ !DEC$ ATTRIBUTES ALIGN: 32 :: W
+ !DEC$ ATTRIBUTES ALIGN: 32 :: work
+ !DEC$ ATTRIBUTES ALIGN: 32 :: ipiv
+
+ allocate(W(LDC,n),work(2*n))
+ tol=tol_in
+
+ info = 0
+ do i=1,n
+  do k=1,i
+   W(i,k) = P(i,k)
+  enddo
+  do k=i+1,n
+   W(i,k) = 0.
+  enddo
+ enddo
+ call DPSTRF('L', n, W, LDC, ipiv, rank, tol, work, info )
+ do i=1,n
+  do k=1,min(m,rank)
+   C(ipiv(i),k) = W(i,k)
+  enddo
+ enddo
+
+ deallocate(W,work)
+end
+
+BEGIN_PROVIDER [ double precision, mo_density_matrix, (mo_tot_num_align, mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Density matrix in MO basis
+ END_DOC
+ integer :: i,j,k
+ mo_density_matrix = 0.d0
+ do k=1,mo_tot_num
+   if (mo_occ(k) == 0.d0) then
+     cycle
+   endif
+   do j=1,ao_num
+     do i=1,ao_num
+       mo_density_matrix(i,j) = mo_density_matrix(i,j) + &
+         mo_occ(k) * mo_coef(i,k) * mo_coef(j,k)
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, mo_density_matrix_virtual, (mo_tot_num_align, mo_tot_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Density matrix in MO basis (virtual MOs)
+ END_DOC
+ integer :: i,j,k
+ mo_density_matrix_virtual = 0.d0
+ do k=1,mo_tot_num
+   do j=1,ao_num
+     do i=1,ao_num
+       mo_density_matrix_virtual(i,j) = mo_density_matrix_virtual(i,j) + &
+         (2.d0-mo_occ(k)) * mo_coef(i,k) * mo_coef(j,k)
+     enddo
+   enddo
+ enddo
+END_PROVIDER
+

src/MOs/mos.ezfio_config

@@ -1,5 +1,7 @@
 mo_basis
+  ao_num             integer
   mo_tot_num         integer
   mo_coef            double precision  (ao_basis_ao_num,mo_basis_mo_tot_num)
   mo_label           character*(64)
+  mo_occ             double precision  (mo_basis_mo_tot_num)
 

src/MOs/mos.irp.f

@@ -32,7 +32,7 @@ END_PROVIDER
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer
   logical                        :: exists
   PROVIDE ezfio_filename
-
+  
   !Label
   call ezfio_has_mo_basis_mo_label(exists)
   if (exists) then
@@ -40,7 +40,7 @@ END_PROVIDER
   else
     mo_label = 'no_label'
   endif
-
+  
   ! Coefs
   allocate(buffer(ao_num,mo_tot_num))
   buffer = 0.d0
@@ -75,3 +75,25 @@ BEGIN_PROVIDER [ double precision, mo_coef_transp, (mo_tot_num_align,ao_num) ]
   
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ]
+  implicit none
+  BEGIN_DOC
+  ! MO occupation numbers
+  END_DOC
+  PROVIDE ezfio_filename
+  logical :: exists
+  call ezfio_has_mo_basis_mo_occ(exists)
+  if (exists) then
+    call ezfio_get_mo_basis_mo_occ(mo_occ)
+  else
+    mo_occ = 0.d0
+    integer :: i
+    do i=1,elec_beta_num
+      mo_occ(i) = 2.d0
+    enddo
+    do i=elec_beta_num+1,elec_alpha_num
+      mo_occ(i) = 1.d0
+    enddo
+  endif
+END_PROVIDER
+

src/MOs/utils.irp.f

@@ -14,6 +14,7 @@ subroutine save_mos
     enddo
   enddo
   call ezfio_set_mo_basis_mo_coef(buffer)
+  call ezfio_set_mo_basis_mo_occ(mo_occ)
   deallocate (buffer)
   
 end
@@ -27,21 +28,31 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label)
   double precision, allocatable  :: mo_coef_new(:,:), R(:,:),eigvalues(:)
   !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
   
+  call write_time(output_mos)
   if (m /= mo_tot_num) then
     print *, irp_here, ': Error : m/= mo_tot_num'
+    stop 1
   endif
-  allocate(R(n,m))
+  allocate(R(n,m),mo_coef_new(ao_num_align,m),eigvalues(m))
-  allocate(mo_coef_new(ao_num_align,m),eigvalues(m))
   mo_coef_new = mo_coef
   
   call lapack_diag(eigvalues,R,matrix,size(matrix,1),size(matrix,2))
   integer :: i
+  write (output_mos,'(A)'), 'MOs are now **'//trim(label)//'**'
+  write (output_mos,'(A)'), ''
+  write (output_mos,'(A)'), 'Eigenvalues'
+  write (output_mos,'(A)'), '-----------'
+  write (output_mos,'(A)'), ''
+  write (output_mos,'(A)'), '======== ================'
   do i = 1, m
-   print*,'eigvalues(i) = ',eigvalues(i)
+   write (output_mos,'(I8,X,F16.10)'), i,eigvalues(i)
   enddo
+  write (output_mos,'(A)'), '======== ================'
+  write (output_mos,'(A)'), ''
   
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
   deallocate(mo_coef_new,R,eigvalues)
+  call write_time(output_mos)
   
   mo_label = label
   SOFT_TOUCH mo_coef mo_label

src/MonoInts/NEEDED_MODULES

@@ -1 +1 @@
-AOs Ezfio_files MOs Nuclei Output Utils
+AOs Electrons Ezfio_files MOs Nuclei Output Utils

src/MonoInts/README.rst

@@ -5,6 +5,7 @@ Needed Modules
 .. NEEDED_MODULES file.
 
 * `AOs <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
+* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
 * `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
 * `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
 * `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_

src/Perturbation/README.rst

@@ -82,7 +82,7 @@ Documentation
 .. Do not edit this section. It was auto-generated from the
 .. NEEDED_MODULES file.
 
-`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L1>`_
+`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L/subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_st)/;">`_
   compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
   .br
   for the various n_st states.
@@ -92,7 +92,7 @@ Documentation
   e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies)
   .br
 
-`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
+`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L/subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
   compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
   .br
   for the various N_st states.
@@ -102,7 +102,7 @@ Documentation
   e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
   .br
 
-`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L40>`_
+`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L/subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
   compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
   .br
   for the various N_st states.
@@ -112,20 +112,46 @@ Documentation
   c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
   .br
 
-`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L1>`_
+`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L/subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,N_st)/;">`_
+  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 CI_electronic_energy(1)
+  .br
+  that could be repeated to this determinant.
+  .br
+  In addition, for the perturbative energetic contribution you have the standard second order
+  .br
+  e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
+  .br
+  and also the purely projected contribution
+  .br
+  H_pert_diag = <HF|H|det_pert> c_pert
+
+`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/pert_sc2.irp.f#L/double precision function repeat_all_e_corr(key_in)/;">`_
+  Undocumented
+
+`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/subroutine fill_H_apply_buffer_selection(n_selected,det_buffer,e_2_pert_buffer,coef_pert_buffer, &>`_
   Fill the H_apply buffer with determiants for the selection
 
-`remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L81>`_
+`remove_small_contributions <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/subroutine remove_small_contributions/;">`_
   Remove determinants with small contributions. N_states is assumed to be
   provided.
 
-`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L68>`_
+`selection_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/BEGIN_PROVIDER [ double precision, selection_criterion ]/;">`_
   Threshold to select determinants. Set by selection routines.
 
-`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L70>`_
+`selection_criterion_factor <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/&BEGIN_PROVIDER [ double precision, selection_criterion_factor ]/;">`_
   Threshold to select determinants. Set by selection routines.
 
-`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L69>`_
+`selection_criterion_min <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L/&BEGIN_PROVIDER [ double precision, selection_criterion_min ]/;">`_
   Threshold to select determinants. Set by selection routines.