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

scripts/cipsi_save.sh

@@ -0,0 +1,27 @@
+#!/bin/bash
+#
+# This script runs a CIPSI calculation as a sequence of single CIPSI iterations.
+# After each iteration, the EZFIO directory is saved.
+#
+# Usage: cipsi_save [EZFIO_FILE] [NDET]
+#
+# Example: cipsi_save file.ezfio 10000
+
+EZ=$1
+NDETMAX=$2
+
+qp set_file ${EZ}
+qp reset -d
+qp set determinants read_wf true
+declare -i NDET
+NDET=1
+while [[ ${NDET} -lt ${NDETMAX} ]]
+do
+        NDET=$(($NDET + $NDET))
+        qp set determinants n_det_max $NDET
+        qp run fci > ${EZ}.out
+        NDET=$(qp get determinants n_det)
+	mv ${EZ}.out ${EZ}.${NDET}.out
+        cp -r ${EZ} ${EZ}.${NDET}
+done
+

src/cipsi/selection.irp.f

@@ -195,7 +195,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
   integer                        :: l_a, nmax, idx
   integer, allocatable           :: indices(:), exc_degree(:), iorder(:)
-  double precision, parameter :: norm_thr = 1.d-16
+
+  ! Removed to avoid introducing determinants already presents in the wf
+  !double precision, parameter :: norm_thr = 1.d-16
+
   allocate (indices(N_det),                                          &
       exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
 
@@ -215,10 +218,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
       i = psi_bilinear_matrix_rows(l_a)
       if (nt + exc_degree(i) <= 4) then
         idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a))
-        if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
+        ! Removed to avoid introducing determinants already presents in the wf
+        !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
           indices(k) = idx
           k=k+1
-        endif
+        !endif
       endif
     enddo
   enddo
@@ -242,10 +246,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
         idx = psi_det_sorted_order(                                  &
             psi_bilinear_matrix_order(                               &
             psi_bilinear_matrix_transp_order(l_a)))
-        if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
+        ! Removed to avoid introducing determinants already presents in the wf
+        !if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
           indices(k) = idx
           k=k+1
-        endif
+        !endif
       endif
     enddo
   enddo
@@ -566,6 +571,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
   double precision, external :: diag_H_mat_elem_fock
   double precision :: E_shift
   double precision :: s_weight(N_states,N_states)
+  logical, external :: is_in_wavefunction
   PROVIDE dominant_dets_of_cfgs N_dominant_dets_of_cfgs
   do jstate=1,N_states
     do istate=1,N_states
@@ -830,8 +836,27 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
             endif
 
         end select
+
+        ! To force the inclusion of determinants with a positive pt2 contribution
+        if (e_pert(istate) > 1d-8) then
+          w = -huge(1.0)
+        endif
+
       end do
 
+!!!BEGIN_DEBUG
+!      ! To check if the pt2 is taking determinants already in the wf
+!      if (is_in_wavefunction(det(N_int,1),N_int)) then
+!        print*, 'A determinant contributing to the pt2 is already in'
+!        print*, 'the wave function:'
+!        call  print_det(det(N_int,1),N_int)
+!        print*,'contribution to the pt2 for the states:', e_pert(:)
+!        print*,'error in the filtering in'
+!        print*, 'cipsi/selection.irp.f sub:  selecte_singles_and_doubles'
+!        print*, 'abort'
+!        call abort
+!      endif
+!!!END_DEBUG
 
       integer(bit_kind) :: occ(N_int,2), n
       if (h0_type == 'CFG') then

src/cis/cis.irp.f

@@ -62,6 +62,7 @@ subroutine run
   else
    call H_apply_cis
   endif
+  print*,''
   print *,  'N_det = ', N_det
   print*,'******************************'
   print *,  'Energies  of the states:'
@@ -69,11 +70,13 @@ subroutine run
     print *,  i, CI_energy(i)
   enddo
   if (N_states > 1) then
-    print*,'******************************'
+    print*,''
-    print*,'Excitation energies '
+    print*,'******************************************************'
+    print*,'Excitation energies (au)                     (eV)'
     do i = 2, N_states
-      print*, i ,CI_energy(i) - CI_energy(1)
+      print*, i ,CI_energy(i) - CI_energy(1), (CI_energy(i) - CI_energy(1))/0.0367502d0
     enddo
+    print*,''
   endif
 
   call ezfio_set_cis_energy(CI_energy)

src/csf/sigma_vector.irp.f

@@ -51,11 +51,24 @@
       if(cfg_seniority_index(i+2) > ncfgpersomo) then
           ncfgpersomo = cfg_seniority_index(i+2)
       else
-        k = 0
+        ! l = i+k+2 
-        do while(cfg_seniority_index(i+2+k) < ncfgpersomo)
+        ! Loop over l with a constraint to ensure that l <= size(cfg_seniority_index,1)-1
-          k = k + 2
+        ! Old version commented just below
-          ncfgpersomo = cfg_seniority_index(i+2+k)
+        do l = min(size(cfg_seniority_index,1)-1, i+2), size(cfg_seniority_index,1)-1, 2
+          if (cfg_seniority_index(l) >= ncfgpersomo) then
+            ncfgpersomo = cfg_seniority_index(l)
+          endif
         enddo
+        !k = 0
+        !if ((i+2+k) < size(cfg_seniority_index,1)) then
+        !  do while(cfg_seniority_index(i+2+k) < ncfgpersomo)
+        !    k = k + 2
+        !    if ((i+2+k) >= size(cfg_seniority_index,1)) then
+        !      exit
+        !    endif
+        !    ncfgpersomo = cfg_seniority_index(i+2+k)
+        !  enddo
+        !endif
       endif
     endif
     ncfg = ncfgpersomo - ncfgprev
@@ -74,11 +87,24 @@
     if(cfg_seniority_index(i+2) > ncfgprev) then
       ncfgprev = cfg_seniority_index(i+2)
     else
-      k = 0
+      ! l = i+k+2 
-      do while(cfg_seniority_index(i+2+k) < ncfgprev)
+      ! Loop over l with a constraint to ensure that l <= size(cfg_seniority_index,1)-1
-        k = k + 2
+      ! Old version commented just below
-        ncfgprev = cfg_seniority_index(i+2+k)
+      do l = min(size(cfg_seniority_index,1)-1, i+2), size(cfg_seniority_index,1)-1, 2
+        if (cfg_seniority_index(l) >= ncfgprev) then
+          ncfgprev = cfg_seniority_index(l)
+        endif
       enddo
+      !k = 0
+      !if ((i+2+k) < size(cfg_seniority_index,1)) then
+      !  do while(cfg_seniority_index(i+2+k) < ncfgprev)
+      !    k = k + 2
+      !    if ((i+2+k) >= size(cfg_seniority_index,1)) then
+      !      exit
+      !    endif
+      !    ncfgprev = cfg_seniority_index(i+2+k)
+      !  enddo
+      !endif
     endif
   enddo
 END_PROVIDER

src/determinants/determinants.irp.f

@@ -77,28 +77,31 @@ BEGIN_PROVIDER [ integer, psi_det_size ]
   END_DOC
   PROVIDE ezfio_filename
   logical                        :: exists
-  if (mpi_master) then
+  psi_det_size = 1
-    call ezfio_has_determinants_n_det(exists)
+  PROVIDE mpi_master
-    if (exists) then
+  if (read_wf) then
-      call ezfio_get_determinants_n_det(psi_det_size)
+    if (mpi_master) then
-    else
+      call ezfio_has_determinants_n_det(exists)
-      psi_det_size = 1
+      if (exists) then
+        call ezfio_get_determinants_n_det(psi_det_size)
+      else
+        psi_det_size = 1
+      endif
+      call write_int(6,psi_det_size,'Dimension of the psi arrays')
     endif
-    call write_int(6,psi_det_size,'Dimension of the psi arrays')
+    IRP_IF MPI_DEBUG
+      print *,  irp_here, mpi_rank
+      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+    IRP_ENDIF
+    IRP_IF MPI
+      include 'mpif.h'
+      integer                        :: ierr
+      call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read psi_det_size with MPI'
+      endif
+    IRP_ENDIF
   endif
-  IRP_IF MPI_DEBUG
-    print *,  irp_here, mpi_rank
-    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-  IRP_ENDIF
-  IRP_IF MPI
-    include 'mpif.h'
-    integer                        :: ierr
-    call MPI_BCAST( psi_det_size, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    if (ierr /= MPI_SUCCESS) then
-      stop 'Unable to read psi_det_size with MPI'
-    endif
-  IRP_ENDIF
-
 
 END_PROVIDER