Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

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_convert_output_to_ezfio

@@ -227,8 +227,8 @@ def write_ezfio(res, filename):
           shell_index += [nshell_tot] * len(b.prim)
 
     shell_num = len(ang_mom)
-    assert(shell_index[0] = 1)
+    assert(shell_index[0] == 1)
-    assert(shell_index[-1] = shell_num)
+    assert(shell_index[-1] == shell_num)
 
     # ~#~#~#~#~ #
     # W r i t e #

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)
 
 

configure

@@ -9,7 +9,7 @@ echo "QP_ROOT="$QP_ROOT
 unset CC
 unset CCXX
 
-TREXIO_VERSION=2.3.2
+TREXIO_VERSION=2.4.2
 
 # Force GCC instead of ICC for dependencies
 export CC=gcc
@@ -219,7 +219,7 @@ EOF
               tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
               cd trexio-${VERSION}
               ./configure --prefix=\${QP_ROOT} --without-hdf5 CFLAGS='-g'
-              make -j 8 && make -j 8 check && make -j 8 install
+              (make -j 8 || make) && make check && make -j 8 install
               tar -zxvf "\${QP_ROOT}"/external/qp2-dependencies/${ARCHITECTURE}/ninja.tar.gz 
               mv ninja "\${QP_ROOT}"/bin/
 EOF
@@ -233,7 +233,7 @@ EOF
               tar -zxf trexio-${VERSION}.tar.gz && rm trexio-${VERSION}.tar.gz
               cd trexio-${VERSION}
               ./configure --prefix=\${QP_ROOT} CFLAGS="-g"
-              make -j 8 && make -j 8 check && make -j 8 install
+              (make -j 8 || make) && make check && make -j 8 install
 EOF
     elif [[ ${PACKAGE} = qmckl ]] ; then
 
@@ -245,7 +245,7 @@ EOF
               tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
               cd qmckl-${VERSION}
               ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc CFLAGS='-g'
-              make && make -j 4 check && make install
+              (make -j 8 || make) && make check && make install
 EOF
     elif [[ ${PACKAGE} = qmckl-intel ]] ; then
 
@@ -257,7 +257,7 @@ EOF
               tar -zxf qmckl-${VERSION}.tar.gz && rm qmckl-${VERSION}.tar.gz
               cd qmckl-${VERSION}
               ./configure --prefix=\${QP_ROOT} --enable-hpc --disable-doc --with-icc --with-ifort CFLAGS='-g'
-              make && make -j 4 check && make install
+              (make -j 8 || make) && make check && make install
 EOF
  
 

external/irpf90

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

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
+        if (imin > bounds(2,isample)) then
-          energy_det = energy_det + sum(memo(bounds(1,isample):bounds(2,isample)))
+          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))
+        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/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
 
@@ -553,7 +562,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
           do l = 1, N_st
             do k = 1, N_st_diag
               do i = 1, sze
-                overlp(k+j-1,l) += U(i,l) * U(i,shift2+k)
+                overlp(k+j-1,l) += u_in(i,l) * U(i,shift2+k)
               enddo
             enddo
           enddo
@@ -576,7 +585,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
             enddo
 
             ! Maximum overlap
-            if (dabs(overlp(k,l)) > omax .and. .not. used .and. state_ok(k)) then
+            if ((dabs(overlp(k,l)) > omax) .and. (.not. used) .and. state_ok(k)) then
               omax = dabs(overlp(k,l))
               idx = k
             endif
@@ -615,6 +624,9 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
         endif
       enddo
 
+      ! Swapped eigenvectors
+      prev_y = y
+
 !      if (state_following) then
 !
 !        overlap = -1.d0
@@ -677,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
@@ -792,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/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
+    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/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

@@ -61,6 +61,58 @@ subroutine export_trexio(update,full_path)
 
 
 
+! ------------------------------------------------------------------------------
+
+! 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
+
 ! ------------------------------------------------------------------------------
 
 ! Electrons

src/utils/memory.irp.f

@@ -107,7 +107,7 @@ subroutine check_mem(rss_in,routine)
   double precision, intent(in) :: rss_in
   character*(*) :: routine
   double precision :: mem
-  call total_memory(mem)
+  call resident_memory(mem)
   mem += rss_in
   if (mem > qp_max_mem) then
     call print_memory_usage()