Fix dos files in qp_create

bugfix davidson recontraction + update

Makefile

@@ -2,4 +2,4 @@ default: build.ninja
 	bash -c "source quantum_package.rc ; ninja"
 
 build.ninja:
-	@bash -c ' echo '' ; echo xxxxxxxxxxxxxxxxxx ; echo "The QP is not configured yet. Please run the ./configure command" ; echo xxxxxxxxxxxxxxxxxx ; echo ''  ; ./configure --help' | more
+	@bash -c ' echo '' ; echo xxxxxxxxxxxxxxxxxx ; echo "QP is not configured yet. Please run the ./configure command" ; echo xxxxxxxxxxxxxxxxxx ; echo ''  ; ./configure --help' | more

bin/qp_plugins

@@ -97,7 +97,7 @@ end
 
 def get_repositories():
     l_result = [f for f in os.listdir(QP_PLUGINS) \
-                  if f not in [".gitignore", "local"] ]
+                  if f not in [".gitignore", "local", "README.rst"] ]
     return sorted(l_result)
 
 

bin/qp_set_frozen_core

@@ -83,6 +83,7 @@ def main(arguments):
                 elif charge <= 118: n_frozen += 43
 
         elif arguments["--small"]:
+            for charge in ezfio.nuclei_nucl_charge:
                 if   charge <=   4: pass
                 elif charge <=  18: n_frozen += 1
                 elif charge <=  36: n_frozen += 5

external/irpf90

@@ -1 +1 @@
-Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
+Subproject commit beac615343f421bd6c0571a408ba389a6d5a32ac

ocaml/Atom.ml

@@ -22,10 +22,15 @@ let of_string ~units s =
     }
   | [ name; x; y; z ] ->
     let e = Element.of_string name in
-    { element = e ;
-      charge  = Element.to_charge e;
-      coord   = Point3d.of_string ~units (String.concat " " [x; y; z])
-    }
+    begin
+      try
+      { element = e ;
+        charge  = Element.to_charge e;
+        coord   = Point3d.of_string ~units (String.concat " " [x; y; z])
+      }
+      with
+      | err -> (Printf.eprintf "name = \"%s\"\nxyz = (%s,%s,%s)\n%!" name x y z ; raise err)
+    end
   | _ -> raise (AtomError s)
 
 

ocaml/Molecule.ml

@@ -142,13 +142,21 @@ let of_xyz_string
   result
 
 
+let regexp_r = Str.regexp {|
|}
+let regexp_t = Str.regexp {|	|}
+
 
 let of_xyz_file
     ?(charge=(Charge.of_int 0)) ?(multiplicity=(Multiplicity.of_int 1))
     ?(units=Units.Angstrom)
     filename =
   let lines =
-    match Io_ext.input_lines filename with
+    Io_ext.input_lines filename
+    |> List.map (fun s -> Str.global_replace regexp_r "" s)
+    |> List.map (fun s -> Str.global_replace regexp_t " " s)
+  in
+  let lines =
+    match lines with
     | natoms :: title :: rest ->
           let natoms = 
             try
@@ -173,6 +181,8 @@ let of_zmt_file
     ?(units=Units.Angstrom)
     filename =
   Io_ext.read_all filename
+  |> Str.global_replace regexp_r ""
+  |> Str.global_replace regexp_t " "
   |> Zmatrix.of_string
   |> Zmatrix.to_xyz_string
   |> of_xyz_string ~charge ~multiplicity ~units

ocaml/Point3d.ml

@@ -24,7 +24,9 @@ let of_string ~units s =
   let l = s
           |> String_ext.split ~on:' '
           |> List.filter (fun x -> x <> "")
-          |> list_map float_of_string
+          |> list_map (fun x ->
+                 try float_of_string x with
+                 | Failure msg -> (Printf.eprintf "Bad string: \"%s\"\n%!" x ; failwith msg) )
           |> Array.of_list
   in
   { x = l.(0) *. f ;

scripts/compilation/qp_create_ninja

@@ -802,8 +802,12 @@ if __name__ == "__main__":
     pickle_path = os.path.join(QP_ROOT, "config", "qp_create_ninja.pickle")
 
     if arguments["update"]:
+         try:
             with open(pickle_path, 'rb') as handle:
                 arguments = pickle.load(handle)
+         except FileNotFoundError:
+            print("\n-----\nError: Please run 'configure -c config/<config_file>'\n-----\n")
+            raise
 
     elif arguments["create"]:
 

src/ao_one_e_ints/ao_ortho_canonical.irp.f

@@ -138,6 +138,8 @@ END_PROVIDER
     deallocate(S)
   endif
 
+  FREE ao_overlap
+
 END_PROVIDER
 
 BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]

src/ccsd/ccsd_t_space_orb_stoch.irp.f

@@ -110,6 +110,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
   double precision :: eocc
   double precision :: norm
   integer :: isample
+  PROVIDE nthreads_pt2
 
 
   ! Prepare table of triplets (a,b,c)
@@ -124,7 +125,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
     do b = a+1, nV
       do c = b+1, nV
         Nabc = Nabc + 1_8
-        Pabc(Nabc) = -1.d0/(f_v(a) + f_v(b) + f_v(c))
+        Pabc(Nabc) =  f_v(a) + f_v(b) + f_v(c)
         abc(1,Nabc) = int(a,2)
         abc(2,Nabc) = int(b,2)
         abc(3,Nabc) = int(c,2)
@@ -134,13 +135,13 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
       abc(1,Nabc) = int(a,2)
       abc(2,Nabc) = int(b,2)
       abc(3,Nabc) = int(a,2)
-      Pabc(Nabc) = -1.d0/(2.d0*f_v(a) + f_v(b))
+      Pabc(Nabc) =  2.d0*f_v(a) + f_v(b)
 
       Nabc = Nabc + 1_8
       abc(1,Nabc) = int(b,2)
       abc(2,Nabc) = int(a,2)
       abc(3,Nabc) = int(b,2)
-      Pabc(Nabc) = -1.d0/(f_v(a) + 2.d0*f_v(b))
+      Pabc(Nabc) =  f_v(a) + 2.d0*f_v(b)
     enddo
   enddo
 
@@ -149,6 +150,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
   enddo
 
   ! Sort triplets in decreasing Pabc
+  Pabc(:) = -1.d0/max(0.2d0,Pabc(:))
   call dsort_big(Pabc, iorder, Nabc)
 
   ! Normalize
@@ -163,7 +165,6 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
 
   call i8set_order_big(abc, iorder, Nabc)
 
-
   ! Cumulative distribution for sampling
   waccu(Nabc) = 0.d0
   do i8=Nabc-1,1,-1
@@ -181,8 +182,8 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
   integer :: nbuckets
   nbuckets = 100
 
+  double precision, allocatable :: ED(:)
   double precision, allocatable :: wsum(:)
-  allocate(wsum(nbuckets))
 
   converged = .False.
   Ncomputed = 0_8
@@ -197,7 +198,8 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
   iright = Nabc
   integer*8, allocatable :: bounds(:,:)
 
-  allocate (bounds(2,nbuckets))
+  allocate(wsum(nbuckets), ED(nbuckets), bounds(2,nbuckets))
+  ED(:) = 0.d0
   do isample=1,nbuckets
     eta = 1.d0/dble(nbuckets) * dble(isample)
     ieta = binary_search(waccu,eta,Nabc)
@@ -215,11 +217,12 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
   print '(A)', ' ======================= ============== =========='
 
 
+  call set_multiple_levels_omp(.False.)
   call wall_time(t00)
   imin = 1_8
   !$OMP PARALLEL                                                     &
       !$OMP PRIVATE(ieta,eta,a,b,c,kiter,isample)                    &
-      !$OMP DEFAULT(SHARED)
+      !$OMP DEFAULT(SHARED) NUM_THREADS(nthreads_pt2)
 
   do kiter=1,Nabc
 
@@ -233,7 +236,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
     enddo
 
     ! Deterministic part
-    if (imin < Nabc) then
+    if (imin <= Nabc) then
       ieta=imin
       sampled(ieta) = 0_8
       a = abc(1,ieta)
@@ -254,7 +257,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
     ! Stochastic part
     call random_number(eta)
     do isample=1,nbuckets
-      if (imin >= bounds(2,isample)) then
+      if (imin > bounds(2,isample)) then
         cycle
       endif
       ieta = binary_search(waccu,(eta + dble(isample-1))/dble(nbuckets),Nabc)+1
@@ -280,7 +283,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
     enddo
 
     call wall_time(t01)
-    if ((t01-t00 > 1.0d0).or.(imin >= Nabc)) then
+    if ((t01-t00 > 1.0d0).or.(imin > Nabc)) then
 
       !$OMP TASKWAIT
       call wall_time(t01)
@@ -300,8 +303,11 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
 
 
       do isample=1,nbuckets
-        if (imin >= bounds(2,isample)) then
-          energy_det = energy_det + sum(memo(bounds(1,isample):bounds(2,isample)))
+        if (imin > bounds(2,isample)) then
+          if (ED(isample) == 0.d0) then
+             ED(isample) = sum(memo(bounds(1,isample):bounds(2,isample)))
+          endif
+          energy_det = energy_det + ED(isample)
           scale = scale - wsum(isample)
         else
           exit
@@ -310,12 +316,14 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
 
       isample = min(isample,nbuckets)
       do ieta=bounds(1,isample), Nabc
-          w = dble(max(sampled(ieta),0_8))
-          tmp = w * memo(ieta) * Pabc(ieta)
-          ET = ET + tmp
-          ET2 = ET2 + tmp * memo(ieta) * Pabc(ieta)
-          norm = norm + w
+        if (sampled(ieta) < 0_8) cycle
+        w = dble(sampled(ieta))
+        tmp = w * memo(ieta) * Pabc(ieta)
+        ET = ET + tmp
+        ET2 = ET2 + tmp * memo(ieta) * Pabc(ieta)
+        norm = norm + w
       enddo
+
       norm = norm/scale
       if (norm > 0.d0) then
         energy_stoch = ET / norm
@@ -327,7 +335,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
       print '(''   '',F20.8, ''   '', ES12.4,''   '', F8.2,''  '')', eccsd+energy, dsqrt(variance/(norm-1.d0)), 100.*real(Ncomputed)/real(Nabc)
     endif
     !$OMP END MASTER
-    if (imin >= Nabc) exit
+    if (imin > Nabc) exit
   enddo
 
   !$OMP END PARALLEL

src/cipsi_utils/pt2_stoch_routines.irp.f

@@ -543,27 +543,59 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
         ! 1/(N-1.5) : see  Brugger, The American Statistician (23) 4 p. 32 (1969)
         if(c > 2) then
           eqt = dabs((pt2_data_S2(t) % pt2(pt2_stoch_istate) / c) - (pt2_data_S(t) % pt2(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
-          eqt = sqrt(eqt / (dble(c) - 1.5d0))
+          eqt = dsqrt(eqt / (dble(c) - 1.5d0))
           pt2_data_err % pt2(pt2_stoch_istate) = eqt
 
           eqt = dabs((pt2_data_S2(t) % variance(pt2_stoch_istate) / c) - (pt2_data_S(t) % variance(pt2_stoch_istate)/c)**2) ! dabs for numerical stability
-          eqt = sqrt(eqt / (dble(c) - 1.5d0))
+          eqt = dsqrt(eqt / (dble(c) - 1.5d0))
           pt2_data_err % variance(pt2_stoch_istate) = eqt
 
           eqta(:) = dabs((pt2_data_S2(t) % overlap(:,pt2_stoch_istate) / c) - (pt2_data_S(t) % overlap(:,pt2_stoch_istate)/c)**2) ! dabs for numerical stability
-          eqta(:) = sqrt(eqta(:) / (dble(c) - 1.5d0))
+          eqta(:) = dsqrt(eqta(:) / (dble(c) - 1.5d0))
           pt2_data_err % overlap(:,pt2_stoch_istate) = eqta(:)
 
 
           if ((time - time1 > 1.d0) .or. (n==N_det_generators)) then
             time1 = time
-            print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.4)', c, &
-              pt2_data     % pt2(pt2_stoch_istate) +E, &
-              pt2_data_err % pt2(pt2_stoch_istate), &
-              pt2_data     % variance(pt2_stoch_istate), &
-              pt2_data_err % variance(pt2_stoch_istate), &
-              pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate), &
-              pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate), &
+
+            value1 =  pt2_data     % pt2(pt2_stoch_istate) + E
+            error1 =  pt2_data_err % pt2(pt2_stoch_istate)
+            value2 =  pt2_data     % pt2(pt2_stoch_istate)
+            error2 =  pt2_data_err % pt2(pt2_stoch_istate)
+            value3 =  pt2_data     % variance(pt2_stoch_istate)
+            error3 =  pt2_data_err % variance(pt2_stoch_istate)
+            value4 =  pt2_data     % overlap(pt2_stoch_istate,pt2_stoch_istate)
+            error4 =  pt2_data_err % overlap(pt2_stoch_istate,pt2_stoch_istate)
+
+            ! Max size of the values (FX.Y) with X=size
+            size1 = 15
+            size2 = 12
+            size3 = 12
+            size4 = 12
+
+            ! To generate the format: number(error)
+            call format_w_error(value1,error1,size1,8,format_value1,str_error1)
+            call format_w_error(value2,error2,size2,8,format_value2,str_error2)
+            call format_w_error(value3,error3,size3,8,format_value3,str_error3)
+            call format_w_error(value4,error4,size4,8,format_value4,str_error4)
+
+            ! value > string with the right format
+            write(str_value1,'('//format_value1//')') value1
+            write(str_value2,'('//format_value2//')') value2
+            write(str_value3,'('//format_value3//')') value3
+            write(str_value4,'('//format_value4//')') value4
+
+            ! Convergence criterion
+            conv_crit = dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &
+                  (1.d-20 + dabs(pt2_data % pt2(pt2_stoch_istate)) )
+            write(str_conv,'(G10.3)') conv_crit
+
+            write(*,'(I10,X,X,A20,X,A16,X,A16,X,A16,X,A12,X,F10.1)') c,&
+            adjustl(adjustr(str_value1)//'('//str_error1(1:1)//')'),&
+            adjustl(adjustr(str_value2)//'('//str_error2(1:1)//')'),&
+            adjustl(adjustr(str_value3)//'('//str_error3(1:1)//')'),&
+            adjustl(adjustr(str_value4)//'('//str_error4(1:1)//')'),&
+            adjustl(str_conv),&
             time-time0
             if (stop_now .or. (                                      &
                   (do_exit .and. (dabs(pt2_data_err % pt2(pt2_stoch_istate)) /    &

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -139,7 +139,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   integer                        :: iter2, itertot
   double precision, allocatable  :: y(:,:), h(:,:), h_p(:,:), lambda(:), s2(:)
   real, allocatable              :: y_s(:,:)
-  double precision, allocatable  :: s_(:,:), s_tmp(:,:)
+  double precision, allocatable  :: s_(:,:), s_tmp(:,:), prev_y(:,:)
   double precision               :: diag_h_mat_elem
   double precision, allocatable  :: residual_norm(:)
   character*(16384)              :: write_buffer
@@ -288,6 +288,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       h(N_st_diag*itermax,N_st_diag*itermax),                        &
 !      h_p(N_st_diag*itermax,N_st_diag*itermax),                      &
       y(N_st_diag*itermax,N_st_diag*itermax),                        &
+      prev_y(N_st_diag*itermax,N_st_diag*itermax),                        &
       s_(N_st_diag*itermax,N_st_diag*itermax),                       &
       s_tmp(N_st_diag*itermax,N_st_diag*itermax),                    &
       residual_norm(N_st_diag),                                      &
@@ -301,6 +302,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
   s_ = 0.d0
   s_tmp = 0.d0
 
+  prev_y = 0.d0
+  do i = 1, N_st_diag*itermax
+    prev_y(i,i) = 1d0
+  enddo
 
   ASSERT (N_st > 0)
   ASSERT (N_st_diag >= N_st)
@@ -479,6 +484,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        if (info > 0) then
          ! Numerical errors propagate. We need to reduce the number of iterations
          itermax = iter-1
+
+         ! eigenvectors of the previous iteration
+         y = prev_y
+         shift2 = shift2 - N_st_diag
          exit
        endif
 
@@ -522,6 +531,84 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         enddo
       endif
 
+      if (state_following) then
+         if (.not. only_expected_s2) then
+           print*,''
+           print*,'!!! State following only available with only_expected_s2 = .True. !!!'
+           STOP
+         endif
+       endif
+
+      if (state_following) then
+
+        integer ::  state(N_st), idx
+        double precision :: omax
+        logical :: used
+        logical, allocatable :: ok(:)
+        double precision, allocatable :: overlp(:,:)
+
+        allocate(overlp(shift2,N_st),ok(shift2))
+
+        overlp = 0d0
+        do j = 1, shift2-1, N_st_diag
+
+          ! Computes some states from the guess vectors
+          ! Psi(:,j:j+N_st_diag) = U y(:,j:j+N_st_diag) and put them
+          ! in U(1,shift2+1:shift2+1+N_st_diag) as temporary array
+          call dgemm('N','N', sze, N_st_diag, shift2,                    &
+          1.d0, U, size(U,1), y(1,j), size(y,1), 0.d0, U(1,shift2+1), size(U,1))
+
+          ! Overlap
+          do l = 1, N_st
+            do k = 1, N_st_diag
+              do i = 1, sze
+                overlp(k+j-1,l) += u_in(i,l) * U(i,shift2+k)
+              enddo
+            enddo
+          enddo
+
+        enddo
+
+        state = 0
+        do l = 1, N_st
+
+          omax = 0d0
+          idx = 0
+          do k = 1, shift2
+
+            ! Already used ?
+            used = .False.
+            do i = 1, N_st
+              if (state(i) == k) then
+                used = .True.
+              endif
+            enddo
+
+            ! Maximum overlap
+            if ((dabs(overlp(k,l)) > omax) .and. (.not. used) .and. state_ok(k)) then
+              omax = dabs(overlp(k,l))
+              idx = k
+            endif
+          enddo
+
+          state(l) = idx
+        enddo
+
+        ! tmp array before setting state_ok
+        ok = .False.
+        do l = 1, N_st
+          ok(state(l)) = .True.
+        enddo
+
+        do k = 1, shift2
+          if (.not. ok(k)) then
+            state_ok(k) = .False.
+          endif
+        enddo
+
+        deallocate(overlp,ok)
+      endif
+
       do k=1,shift2
         if (.not. state_ok(k)) then
           do l=k+1,shift2
@@ -537,46 +624,49 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         endif
       enddo
 
-      if (state_following) then
+      ! Swapped eigenvectors
+      prev_y = y
 
-        overlap = -1.d0
-        do k=1,shift2
-          do i=1,shift2
-            overlap(k,i) = dabs(y(k,i))
-          enddo
-        enddo
-        do k=1,N_st
-          cmax = -1.d0
-          do i=1,N_st
-            if (overlap(i,k) > cmax) then
-              cmax = overlap(i,k)
-              order(k) = i
-            endif
-          enddo
-          do i=1,N_st_diag
-            overlap(order(k),i) = -1.d0
-          enddo
-        enddo
-        overlap = y
-        do k=1,N_st
-          l = order(k)
-          if (k /= l) then
-            y(1:shift2,k) = overlap(1:shift2,l)
-          endif
-        enddo
-        do k=1,N_st
-          overlap(k,1) = lambda(k)
-          overlap(k,2) = s2(k)
-        enddo
-        do k=1,N_st
-          l = order(k)
-          if (k /= l) then
-            lambda(k) = overlap(l,1)
-            s2(k) = overlap(l,2)
-          endif
-        enddo
-
-      endif
+!      if (state_following) then
+!
+!        overlap = -1.d0
+!        do k=1,shift2
+!          do i=1,shift2
+!            overlap(k,i) = dabs(y(k,i))
+!          enddo
+!        enddo
+!        do k=1,N_st
+!          cmax = -1.d0
+!          do i=1,N_st
+!            if (overlap(i,k) > cmax) then
+!              cmax = overlap(i,k)
+!              order(k) = i
+!            endif
+!          enddo
+!          do i=1,N_st_diag
+!            overlap(order(k),i) = -1.d0
+!          enddo
+!        enddo
+!        overlap = y
+!        do k=1,N_st
+!          l = order(k)
+!          if (k /= l) then
+!            y(1:shift2,k) = overlap(1:shift2,l)
+!          endif
+!        enddo
+!        do k=1,N_st
+!          overlap(k,1) = lambda(k)
+!          overlap(k,2) = s2(k)
+!        enddo
+!        do k=1,N_st
+!          l = order(k)
+!          if (k /= l) then
+!            lambda(k) = overlap(l,1)
+!            s2(k) = overlap(l,2)
+!          endif
+!        enddo
+!
+!      endif
 
 
       ! Express eigenvectors of h in the determinant basis
@@ -599,7 +689,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         do i=1,sze
           U(i,shift2+k) =  &
             (lambda(k) * U(i,shift2+k) - W(i,shift2+k) )      &
-              /max(H_jj(i) - lambda (k),1.d-2)
+              /max(dabs(H_jj(i) - lambda (k)),1.d-2) * dsign(1d0,H_jj(i) - lambda (k))
         enddo
 
         if (k <= N_st) then
@@ -714,7 +804,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
       residual_norm,                                                 &
       U, overlap,                                                    &
       h, y_s, S_d,                                                   &
-      y, s_, s_tmp,                                                  &
+      y, s_, s_tmp, prev_y,                                          &
       lambda                                                         &
       )
   FREE nthreads_davidson

src/davidson/diagonalize_ci.irp.f

@@ -123,6 +123,7 @@ END_PROVIDER
         endif
 
      enddo
+
      if (N_states_diag > N_states_diag_save) then
        N_states_diag = N_states_diag_save
        TOUCH N_states_diag
@@ -133,24 +134,101 @@ END_PROVIDER
      print *,  'Diagonalization of H using Lapack'
      allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
      allocate (eigenvalues(N_det))
+
      if (s2_eig) then
+
        double precision, parameter :: alpha = 0.1d0
        allocate (H_prime(N_det,N_det) )
+
        H_prime(1:N_det,1:N_det) = H_matrix_all_dets(1:N_det,1:N_det) +  &
          alpha * S2_matrix_all_dets(1:N_det,1:N_det)
+
        do j=1,N_det
          H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
        enddo
+
        call lapack_diag(eigenvalues,eigenvectors,H_prime,size(H_prime,1),N_det)
        call nullify_small_elements(N_det,N_det,eigenvectors,size(eigenvectors,1),1.d-12)
+
        CI_electronic_energy(:) = 0.d0
        i_state = 0
+
        allocate (s2_eigvalues(N_det))
        allocate(index_good_state_array(N_det),good_state_array(N_det))
+
        good_state_array = .False.
        call u_0_S2_u_0(s2_eigvalues,eigenvectors,N_det,psi_det,N_int,&
          N_det,size(eigenvectors,1))
-       if (only_expected_s2) then
+
+       if (state_following) then
+         if (.not. only_expected_s2) then
+           print*,''
+           print*,'!!! State following only available with only_expected_s2 = .True. !!!'
+           STOP
+         endif
+         if (N_det < N_states) then
+           print*,''
+           print*,'!!! State following requires at least N_states determinants to be activated !!!'
+           STOP
+         endif
+       endif
+
+       if (state_following .and. only_expected_s2) then
+
+         integer :: state(N_states), idx,l
+         double precision :: omax
+         double precision, allocatable :: overlp(:)
+         logical :: used
+         logical, allocatable :: ok(:)
+
+         allocate(overlp(N_det), ok(N_det))
+
+         i_state = 0
+         state = 0
+         do l = 1, N_states
+
+           ! Overlap wrt each state
+           overlp = 0d0
+           do k = 1, N_det
+             do i = 1, N_det
+               overlp(k) = overlp(k) + psi_coef(i,l) * eigenvectors(i,k)
+             enddo
+           enddo
+
+           ! Idx of the state with the maximum overlap not already "used"
+           omax = 0d0
+           idx = 0
+           do k = 1, N_det
+
+             ! Already used ?
+             used = .False.
+             do i = 1, N_states
+               if (state(i) == k) then
+                 used = .True.
+               endif
+             enddo
+
+             ! Maximum overlap
+             if (dabs(overlp(k)) > omax .and. .not. used) then
+               if (dabs(s2_eigvalues(k)-expected_s2) > 0.5d0) cycle
+               omax = dabs(overlp(k))
+               idx = k
+             endif
+           enddo
+
+           state(l) = idx
+           i_state +=1
+         enddo
+
+         deallocate(overlp, ok)
+
+         do i = 1, i_state
+           index_good_state_array(i) = state(i)
+           good_state_array(i) = .True.
+         enddo
+
+       else if (only_expected_s2) then
+
          do j=1,N_det
            ! Select at least n_states states with S^2 values closed to "expected_s2"
            if(dabs(s2_eigvalues(j)-expected_s2).le.0.5d0)then
@@ -158,17 +236,23 @@ END_PROVIDER
              index_good_state_array(i_state) = j
              good_state_array(j) = .True.
            endif
+
            if(i_state.eq.N_states) then
              exit
            endif
          enddo
+
        else
+
          do j=1,N_det
            index_good_state_array(j) = j
            good_state_array(j) = .True.
          enddo
+
        endif
+
        if(i_state .ne.0)then
+
          ! Fill the first "i_state" states that have a correct S^2 value
          do j = 1, i_state
            do i=1,N_det
@@ -177,6 +261,7 @@ END_PROVIDER
            CI_electronic_energy(j) = eigenvalues(index_good_state_array(j))
            CI_s2(j) = s2_eigvalues(index_good_state_array(j))
          enddo
+
          i_other_state = 0
          do j = 1, N_det
            if(good_state_array(j))cycle
@@ -201,6 +286,7 @@ END_PROVIDER
          print*,'  as the CI_eigenvectors'
          print*,'  You should consider more states and maybe ask for s2_eig to be .True. or just enlarge the CI space'
          print*,''
+
          do j=1,min(N_states_diag,N_det)
            do i=1,N_det
              CI_eigenvectors(i,j) = eigenvectors(i,j)
@@ -209,14 +295,18 @@ END_PROVIDER
            CI_s2(j) = s2_eigvalues(j)
          enddo
        endif
+
        deallocate(index_good_state_array,good_state_array)
        deallocate(s2_eigvalues)
+
      else
+
        call lapack_diag(eigenvalues,eigenvectors,                    &
            H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
        CI_electronic_energy(:) = 0.d0
        call u_0_S2_u_0(CI_s2,eigenvectors,N_det,psi_det,N_int,       &
            min(N_det,N_states_diag),size(eigenvectors,1))
+
        ! Select the "N_states_diag" states of lowest energy
        do j=1,min(N_det,N_states_diag)
          do i=1,N_det
@@ -224,7 +314,9 @@ END_PROVIDER
          enddo
          CI_electronic_energy(j) = eigenvalues(j)
        enddo
+
      endif
+
      do k=1,N_states_diag
        CI_electronic_energy(k) = 0.d0
        do j=1,N_det
@@ -235,6 +327,7 @@ END_PROVIDER
          enddo
        enddo
      enddo
+
      deallocate(eigenvectors,eigenvalues)
    endif
 

src/mo_basis/utils.irp.f

@@ -228,7 +228,11 @@ subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label)
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
 
   do i=1,m
-    eig(i) = D(i)
+    if (eig(i) > 1.d-20) then
+      eig(i) = D(i)
+    else
+      eig(i) = 0.d0
+    endif
   enddo
 
   deallocate(A,mo_coef_new,U,Vt,D)

src/mp2/H_apply.irp.f

@@ -0,0 +1,15 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python3 ]
+from generate_h_apply import *
+from perturbation import perturbations
+
+s = H_apply("mp2")
+s.set_perturbation("Moller_plesset")
+#s.set_perturbation("epstein_nesbet")
+print(s)
+
+s = H_apply("mp2_selection")
+s.set_selection_pt2("Moller_Plesset")
+print(s)
+END_SHELL
+

src/mp2/NEED

@@ -0,0 +1,6 @@
+generators_full
+selectors_full
+determinants
+davidson
+davidson_undressed
+perturbation

src/mp2/README.rst

@@ -0,0 +1,4 @@
+===
+mp2
+===
+

src/mp2/mp2.irp.f

@@ -0,0 +1,21 @@
+program mp2
+  call run
+end
+
+subroutine run
+  implicit none
+  double precision, allocatable  :: pt2(:), norm_pert(:)
+  double precision               :: H_pert_diag, E_old
+  integer                        :: N_st, iter
+  PROVIDE Fock_matrix_diag_mo H_apply_buffer_allocated
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st))
+  E_old = HF_energy
+  call H_apply_mp2(pt2, norm_pert, H_pert_diag,  N_st)
+  print *,  'N_det    = ', N_det
+  print *,  'N_states = ', N_states
+  print *,  'MP2      = ', pt2
+  print *,  'E        = ', E_old
+  print *,  'E+MP2    = ', E_old+pt2
+  deallocate(pt2,norm_pert)
+end

src/scf_utils/diagonalize_fock.irp.f

@@ -47,7 +47,7 @@ BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_mo, (ao_num,mo_num)
        do j = 1, n_core_orb
          jorb = list_core(j)
          F(iorb,jorb) = 0.d0
-         F(jorb,iorb) = 0.d0                                                                                                                 
+         F(jorb,iorb) = 0.d0
        enddo
      enddo
    endif

src/scf_utils/roothaan_hall_scf.irp.f

@@ -13,9 +13,9 @@ END_DOC
   integer                        :: iteration_SCF,dim_DIIS,index_dim_DIIS
 
   logical                        :: converged
-  integer                        :: i,j
+  integer                        :: i,j,m
   logical, external              :: qp_stop
-  double precision, allocatable :: mo_coef_save(:,:)
+  double precision, allocatable :: mo_coef_save(:,:), S(:,:)
 
   PROVIDE ao_md5 mo_occ level_shift
 
@@ -208,9 +208,29 @@ END_DOC
       size(Fock_matrix_mo,2),mo_label,1,.true.)
    call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
    call orthonormalize_mos
-   call save_mos
   endif
 
+
+  ! Identify degenerate MOs and force them on the axes
+  allocate(S(ao_num,ao_num))
+  i=1
+  do while (i<mo_num)
+    j=i
+    m=1
+    do while ( (j<mo_num).and.(fock_matrix_diag_mo(j+1)-fock_matrix_diag_mo(i) < 1.d-8) )
+      j += 1
+      m += 1
+    enddo
+    if (m>1) then
+      call dgemm('N','T',ao_num,ao_num,m,1.d0,mo_coef(1,i),size(mo_coef,1),mo_coef(1,i),size(mo_coef,1),0.d0,S,size(S,1))
+      call pivoted_cholesky( S, m, -1.d0, ao_num, mo_coef(1,i))
+    endif
+    i = j+1
+  enddo
+
+
+  call save_mos
+
   call write_double(6, Energy_SCF, 'SCF energy')
 
   call write_time(6)

src/trexio/export_trexio_routines.irp.f

@@ -59,7 +59,59 @@ subroutine export_trexio(update,full_path)
   enddo
   call ezfio_set_trexio_trexio_file(trexio_filename)
 
-    
+
+
+! ------------------------------------------------------------------------------
+
+! Metadata
+! --------
+
+  integer :: code_num, author_num
+  character*(64) :: code(100), author(100), user
+  character*(64), parameter :: qp2_code = "QuantumPackage"
+
+  call getenv("USER",user)
+  do k=1,N_states
+    rc = trexio_read_metadata_code_num(f(k), code_num)
+    if (rc == TREXIO_ATTR_MISSING) then
+      i = 1
+      code(:) = ""
+    else
+      rc = trexio_read_metadata_code(f(k), code, 64)
+      do i=1, code_num
+        if (trim(code(i)) == trim(qp2_code)) then
+          exit
+        endif
+      enddo
+    endif
+    if (i == code_num+1) then
+      code(i) = qp2_code
+      rc = trexio_write_metadata_code_num(f(k), i)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+      rc = trexio_write_metadata_code(f(k), code, 64)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
+
+    rc = trexio_read_metadata_author_num(f(k), author_num)
+    if (rc == TREXIO_ATTR_MISSING) then
+      i = 1
+      author(:) = ""
+    else
+      rc = trexio_read_metadata_author(f(k), author, 64)
+      do i=1, author_num
+        if (trim(author(i)) == trim(user)) then
+          exit
+        endif
+      enddo
+    endif
+    if (i == author_num+1) then
+      author(i) = user
+      rc = trexio_write_metadata_author_num(f(k), i)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+      rc = trexio_write_metadata_author(f(k), author, 64)
+      call trexio_assert(rc, TREXIO_SUCCESS)
+    endif
+  enddo
 
 ! ------------------------------------------------------------------------------