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

bin/qp_set_frozen_core

@@ -7,12 +7,13 @@ setting all MOs as Active, except the n/2 first ones which are set as Core.
 If pseudo-potentials are used, all the MOs are set as Active.
 
 Usage:
-      qp_set_frozen_core [-q|--query] [(-l|-s|--large|--small)] EZFIO_DIR
+      qp_set_frozen_core [-q|--query] [(-l|-s|-u|--large|--small|--unset)] EZFIO_DIR
 
 Options:
     -q --query   Prints in the standard output the number of frozen MOs
     -l --large   Use a small core
     -s --small   Use a large core
+    -u --unset   Unset frozen core
 
 
 Default numbers of frozen electrons:
@@ -88,7 +89,9 @@ def main(arguments):
                 elif charge <=  54: n_frozen += 9
                 elif charge <=  86: n_frozen += 18
                 elif charge <= 118: n_frozen += 27
+        elif arguments["--unset"]:
 
+              n_frozen = 0
         else:  # default                                                                      
             for charge in ezfio.nuclei_nucl_charge:
                 if   charge <=   4: pass

external/qp2-dependencies

@@ -1 +1 @@
-Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0
+Subproject commit bc856147f6e626a6616b20344e5b8e3f30f44a92

src/cipsi/pt2_stoch_routines.irp.f

@@ -537,7 +537,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
 
           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, &
+            print '(I10, X, F12.6, X, G10.3, X, F10.6, X, G10.3, X, F10.6, X, G10.3, X, F10.1)', c, &
               pt2_data     % pt2(pt2_stoch_istate) +E, &
               pt2_data_err % pt2(pt2_stoch_istate), &
               pt2_data     % variance(pt2_stoch_istate), &

src/cipsi/run_pt2_slave.irp.f

@@ -116,10 +116,10 @@ subroutine run_pt2_slave_small(thread,iproc,energy)
     do k=1,n_tasks
       call pt2_alloc(pt2_data(k),N_states)
       b%cur = 0
-      double precision :: time2
+!      double precision :: time2
-      call wall_time(time2)
+!      call wall_time(time2)
       call select_connected(i_generator(k),energy,pt2_data(k),b,subset(k),pt2_F(i_generator(k)))
-      call wall_time(time1)
+!      call wall_time(time1)
 !      print *,  i_generator(1), time1-time2, n_tasks, pt2_F(i_generator(1))
     enddo
     call wall_time(time1)
@@ -190,8 +190,12 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
 
   integer :: bsize ! Size of selection buffers
   logical :: sending
+  double precision :: time_shift
+
   PROVIDE global_selection_buffer global_selection_buffer_lock
 
+  call random_number(time_shift)
+  time_shift = time_shift*15.d0
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
 
@@ -209,6 +213,9 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
 
   sending = .False.
   done = .False.
+  double precision :: time0, time1
+  call wall_time(time0)
+  time0 = time0+time_shift
   do while (.not.done)
 
     integer, external :: get_tasks_from_taskserver
@@ -244,19 +251,24 @@ subroutine run_pt2_slave_large(thread,iproc,energy)
       done = .true.
     endif
     call sort_selection_buffer(b)
+
+    call wall_time(time1)
+!    if (time1-time0 > 15.d0) then
       call omp_set_lock(global_selection_buffer_lock)
       global_selection_buffer%mini = b%mini
       call merge_selection_buffers(b,global_selection_buffer)
       b%cur=0
       call omp_unset_lock(global_selection_buffer_lock)
+      call wall_time(time0)
+!    endif
+
+    call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending)
     if ( iproc == 1 .or. i_generator < 100 .or. done) then
       call omp_set_lock(global_selection_buffer_lock)
-      call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending)
       call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), global_selection_buffer, (/task_id/), 1,sending)
       global_selection_buffer%cur = 0
       call omp_unset_lock(global_selection_buffer_lock)
     else
-      call push_pt2_results_async_recv(zmq_socket_push,b%mini,sending)
       call push_pt2_results_async_send(zmq_socket_push, (/i_generator/), (/pt2_data/), b, (/task_id/), 1,sending)
     endif
 

src/cipsi/selection.irp.f

@@ -474,17 +474,11 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 !      endif
 
       do i=1,fullinteresting(0)
-        do k=1,N_int
+        fullminilist(:,:,i) = psi_det_sorted(:,:,fullinteresting(i))
-          fullminilist(k,1,i) = psi_det_sorted(k,1,fullinteresting(i))
-          fullminilist(k,2,i) = psi_det_sorted(k,2,fullinteresting(i))
-        enddo
       enddo
 
       do i=1,interesting(0)
-        do k=1,N_int
+        minilist(:,:,i) = psi_det_sorted(:,:,interesting(i))
-          minilist(k,1,i) = psi_det_sorted(k,1,interesting(i))
-          minilist(k,2,i) = psi_det_sorted(k,2,interesting(i))
-        enddo
       enddo
 
       do s2=s1,2

src/cipsi/selection_buffer.irp.f

@@ -92,38 +92,51 @@ subroutine merge_selection_buffers(b1, b2)
   allocate(val(sze), detmp(N_int, 2, sze))
   i1=1
   i2=1
+
+  select case (N_int)
+BEGIN_TEMPLATE 
+  case $case
     do i=1,nmwen
       if ( (i1 > b1%cur).and.(i2 > b2%cur) ) then
         exit
       else if (i1 > b1%cur) then
           val(i) = b2%val(i2)
-        detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2)
+          detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2)
-        detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2)
+          detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2)
           i2=i2+1
       else if (i2 > b2%cur) then
           val(i) = b1%val(i1)
-        detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1)
+          detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1)
-        detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1)
+          detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1)
           i1=i1+1
       else
         if (b1%val(i1) <= b2%val(i2)) then
           val(i) = b1%val(i1)
-        detmp(1:N_int,1,i) = b1%det(1:N_int,1,i1)
+          detmp(1:$N_int,1,i) = b1%det(1:$N_int,1,i1)
-        detmp(1:N_int,2,i) = b1%det(1:N_int,2,i1)
+          detmp(1:$N_int,2,i) = b1%det(1:$N_int,2,i1)
           i1=i1+1
         else
           val(i) = b2%val(i2)
-        detmp(1:N_int,1,i) = b2%det(1:N_int,1,i2)
+          detmp(1:$N_int,1,i) = b2%det(1:$N_int,1,i2)
-        detmp(1:N_int,2,i) = b2%det(1:N_int,2,i2)
+          detmp(1:$N_int,2,i) = b2%det(1:$N_int,2,i2)
           i2=i2+1
         endif
       endif
     enddo
-  deallocate(b2%det, b2%val)
     do i=nmwen+1,b2%N
       val(i) = 0.d0
-    detmp(1:N_int,1:2,i) = 0_bit_kind
+!      detmp(1:$N_int,1,i) = 0_bit_kind
+!      detmp(1:$N_int,2,i) = 0_bit_kind
     enddo
+SUBST [ case, N_int ]
+(1); 1;;
+(2); 2;;
+(3); 3;;
+(4); 4;;
+default; N_int;;
+END_TEMPLATE
+  end select
+  deallocate(b2%det, b2%val)
   b2%det => detmp
   b2%val => val
   b2%mini = min(b2%mini,b2%val(b2%N))

src/csf/configurations.irp.f

@@ -779,6 +779,7 @@ subroutine binary_search_cfg(cfgInp,addcfg)
 end subroutine
 
  BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det, (2,N_configuration) ]
+&BEGIN_PROVIDER [ integer, psi_configuration_n_det, (N_configuration) ]
 &BEGIN_PROVIDER [ integer, psi_configuration_to_psi_det_data, (N_det) ]
 
  implicit none
@@ -867,6 +868,29 @@ end subroutine
  enddo
 
  deallocate(dets, old_order)
+ integer :: ndet_conf
+ do i = 1, N_configuration
+  ndet_conf = psi_configuration_to_psi_det(2,i) - psi_configuration_to_psi_det(1,i) + 1
+  psi_configuration_n_det(i) = ndet_conf
+ enddo
 
 END_PROVIDER
 
+
+BEGIN_PROVIDER [ integer, n_elec_alpha_for_psi_configuration, (N_configuration)]
+ implicit none
+ integer :: i,j,k,l
+ integer(bit_kind) :: det_tmp(N_int,2),det_alpha(N_int)
+ n_elec_alpha_for_psi_configuration = 0
+ do i = 1, N_configuration
+  j = psi_configuration_to_psi_det(2,i) 
+  det_tmp(:,:) = psi_det(:,:,j)
+  k = 0
+  do l = 1, N_int
+   det_alpha(N_int) = iand(det_tmp(l,1),psi_configuration(l,1,i))
+   k += popcnt(det_alpha(l))
+  enddo
+  n_elec_alpha_for_psi_configuration(i) = k
+ enddo
+
+END_PROVIDER 

src/dav_general_mat/dav_diag_dressed_ext_rout.irp.f

@@ -1,5 +1,5 @@
 
-subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
+subroutine davidson_general_ext_rout_diag_dressed(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_diag_in,converged,hcalc)
   use mmap_module
   implicit none
   BEGIN_DOC
@@ -412,36 +412,6 @@ subroutine davidson_general_ext_rout(u_in,H_jj,Dress_jj,energies,sze,N_st,N_st_d
   FREE nthreads_davidson
 end
 
-subroutine hcalc_template(v,u,N_st,sze)
-  use bitmasks
-  implicit none
-  BEGIN_DOC
-  ! Template of routine for the application of H
-  !
-  ! Here, it is done with the Hamiltonian matrix 
-  !
-  ! on the set of determinants of psi_det 
-  !
-  ! Computes $v = H | u \rangle$ 
-  !
-  END_DOC
-  integer, intent(in)              :: N_st,sze
-  double precision, intent(in)     :: u(sze,N_st)
-  double precision, intent(inout)  :: v(sze,N_st)
-  integer :: i,j,istate
-  v = 0.d0
-  do istate = 1, N_st
-   do i = 1, sze
-    do j = 1, sze
-      v(i,istate) += H_matrix_all_dets(j,i) * u(j,istate)
-    enddo
-   enddo
-   do i = 1, sze
-    v(i,istate) += u(i,istate) * nuclear_repulsion
-   enddo
-  enddo
-end
-
 subroutine dressing_diag_uv(v,u,dress_diag,N_st,sze)
   implicit none
   BEGIN_DOC

src/davidson/u0_hs2_u0.irp.f

@@ -203,7 +203,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   integer, allocatable           :: doubles(:)
   integer, allocatable           :: singles_a(:)
   integer, allocatable           :: singles_b(:)
-  integer, allocatable           :: idx(:), idx0(:)
+  integer, allocatable           :: idx(:), buffer_lrow(:), idx0(:)
   integer                        :: maxab, n_singles_a, n_singles_b, kcol_prev
   integer*8                      :: k8
   logical                        :: compute_singles
@@ -253,7 +253,7 @@ compute_singles=.True.
       !$OMP   PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i,     &
       !$OMP          lcol, lrow, l_a, l_b, utl, kk, u_is_sparse,     &
       !$OMP          buffer, doubles, n_doubles, umax,               &
-      !$OMP          tmp_det2, hij, sij, idx, l, kcol_prev,          &
+      !$OMP          tmp_det2, hij, sij, idx, buffer_lrow, l, kcol_prev,          &
       !$OMP          singles_a, n_singles_a, singles_b, ratio,       &
       !$OMP          n_singles_b, k8, last_found,left,right,right_max)
 
@@ -264,7 +264,7 @@ compute_singles=.True.
       singles_a(maxab),                                              &
       singles_b(maxab),                                              &
       doubles(maxab),                                                &
-      idx(maxab), utl(N_st,block_size))
+      idx(maxab), buffer_lrow(maxab), utl(N_st,block_size))
 
   kcol_prev=-1
 
@@ -332,18 +332,20 @@ compute_singles=.True.
         l_a = psi_bilinear_matrix_columns_loc(lcol)
         ASSERT (l_a <= N_det)
 
-        !DIR$ UNROLL(8)
-        !DIR$ LOOP COUNT avg(50000)
         do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
           lrow = psi_bilinear_matrix_rows(l_a)
           ASSERT (lrow <= N_det_alpha_unique)
 
-          buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, lrow)  ! hot spot
+          buffer_lrow(j) = lrow
 
           ASSERT (l_a <= N_det)
           idx(j) = l_a
           l_a = l_a+1
         enddo
+
+        do j=1,psi_bilinear_matrix_columns_loc(lcol+1) - psi_bilinear_matrix_columns_loc(lcol)
+          buffer(1:$N_int,j) = psi_det_alpha_unique(1:$N_int, buffer_lrow(j))  ! hot spot
+        enddo
         j = j-1
 
         call get_all_spin_singles_$N_int(                              &
@@ -789,7 +791,7 @@ compute_singles=.True.
 
   end do
   !$OMP END DO
-  deallocate(buffer, singles_a, singles_b, doubles, idx, utl)
+  deallocate(buffer, singles_a, singles_b, doubles, idx, buffer_lrow, utl)
   !$OMP END PARALLEL
 
 end

src/determinants/s2.irp.f

@@ -103,13 +103,17 @@ BEGIN_PROVIDER [ double precision, expected_s2]
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, s2_values, (N_states) ]
+ BEGIN_PROVIDER [ double precision, s2_values, (N_states) ]
+&BEGIN_PROVIDER [ double precision, s_values, (N_states) ]
  implicit none
  BEGIN_DOC
 ! array of the averaged values of the S^2 operator on the various states
  END_DOC
  integer :: i
  call u_0_S2_u_0(s2_values,psi_coef,n_det,psi_det,N_int,N_states,psi_det_size)
+ do i = 1, N_states
+  s_values(i) = 0.5d0 *(-1.d0 + dsqrt(1.d0 + 4 * s2_values(i)))
+ enddo
 
 END_PROVIDER
 

src/dft_utils_in_r/ao_in_r.irp.f

@@ -91,7 +91,19 @@
  enddo
  END_PROVIDER
 
- BEGIN_PROVIDER[double precision, aos_lapl_in_r_array, (ao_num,n_points_final_grid,3)]
+ BEGIN_PROVIDER [double precision, aos_lapl_in_r_array_transp, (ao_num, n_points_final_grid,3)]
+ implicit none
+ integer :: i,j,m
+ do i = 1, n_points_final_grid
+  do j = 1, ao_num
+   do m = 1, 3
+    aos_lapl_in_r_array_transp(j,i,m) =  aos_lapl_in_r_array(m,j,i)
+   enddo
+  enddo
+ enddo
+ END_PROVIDER 
+
+ BEGIN_PROVIDER [double precision, aos_lapl_in_r_array, (3,ao_num,n_points_final_grid)]
  implicit none
  BEGIN_DOC
  ! aos_lapl_in_r_array(i,j,k)   = value of the kth component of the laplacian of jth ao on the ith grid point
@@ -100,20 +112,20 @@
  END_DOC
  integer :: i,j,m
  double precision :: aos_array(ao_num), r(3)
- double precision :: aos_grad_array(ao_num,3)
+ double precision :: aos_grad_array(3,ao_num)
- double precision :: aos_lapl_array(ao_num,3)
+ double precision :: aos_lapl_array(3,ao_num)
  !$OMP PARALLEL DO &
  !$OMP DEFAULT (NONE)  &
  !$OMP PRIVATE (i,r,aos_array,aos_grad_array,aos_lapl_array,j,m) & 
  !$OMP SHARED(aos_lapl_in_r_array,n_points_final_grid,ao_num,final_grid_points)
- do m = 1, 3
  do i = 1, n_points_final_grid
   r(1) = final_grid_points(1,i)
   r(2) = final_grid_points(2,i)
   r(3) = final_grid_points(3,i)
   call give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_lapl_array)
   do j = 1, ao_num
-    aos_lapl_in_r_array(j,i,m) = aos_lapl_array(j,m)
+   do m = 1, 3
+    aos_lapl_in_r_array(m,j,i) = aos_lapl_array(m,j)
    enddo
   enddo
  enddo

src/dft_utils_in_r/mo_in_r.irp.f

@@ -138,7 +138,7 @@
  integer :: m
  mos_lapl_in_r_array = 0.d0
  do m=1,3
-  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num)
+  call dgemm('N','N',mo_num,n_points_final_grid,ao_num,1.d0,mo_coef_transp,mo_num,aos_lapl_in_r_array_transp(1,1,m),ao_num,0.d0,mos_lapl_in_r_array(1,1,m),mo_num)
  enddo
  END_PROVIDER
 

src/ezfio_files/output.irp.f

@@ -25,7 +25,7 @@ subroutine write_time(iunit)
   ct = ct - output_cpu_time_0
   call wall_time(wt)
   wt = wt - output_wall_time_0
-  write(6,'(A,F14.6,A,F14.6,A)') &
+  write(6,'(A,F14.2,A,F14.2,A)') &
     '.. >>>>> [ WALL TIME: ', wt, '  s ] [ CPU  TIME: ', ct, '  s ] <<<<< ..'
   write(6,*)
 end

src/mo_guess/h_core_guess_routine.irp.f

@@ -7,7 +7,7 @@ subroutine hcore_guess
   label = 'Guess'
   call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals,          &
                                      size(mo_one_e_integrals,1),  &
-                                     size(mo_one_e_integrals,2),label,1,.false.)
+                                     size(mo_one_e_integrals,2),label,1,.true.)
   call nullify_small_elements(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-12 )
   call save_mos
   TOUCH mo_coef mo_label

src/tools/print_wf.irp.f

@@ -32,8 +32,9 @@ subroutine routine
  double precision :: norm_mono_a,norm_mono_b
  double precision :: norm_mono_a_2,norm_mono_b_2
  double precision :: norm_mono_a_pert_2,norm_mono_b_pert_2
- double precision :: norm_mono_a_pert,norm_mono_b_pert
+ double precision :: norm_mono_a_pert,norm_mono_b_pert,norm_double_1
  double precision :: delta_e,coef_2_2
+
  norm_mono_a = 0.d0
  norm_mono_b = 0.d0
  norm_mono_a_2 = 0.d0
@@ -42,6 +43,7 @@ subroutine routine
  norm_mono_b_pert = 0.d0
  norm_mono_a_pert_2 = 0.d0
  norm_mono_b_pert_2 = 0.d0
+ norm_double_1 = 0.d0
  do i = 1, min(N_det_print_wf,N_det)
   print*,''
   print*,'i = ',i
@@ -93,6 +95,7 @@ subroutine routine
     print*,'h1,p1 = ',h1,p1
     print*,'s2',s2
     print*,'h2,p2 = ',h2,p2
+    norm_double_1 += dabs(psi_coef_sorted(i,1)/psi_coef_sorted(1,1))
    endif
 
    print*,'<Ref| H |D_I> = ',hij
@@ -109,6 +112,7 @@ subroutine routine
  print*,''
  print*,'L1 norm of mono alpha = ',norm_mono_a
  print*,'L1 norm of mono beta  = ',norm_mono_b
+ print*,'L1 norm of double exc = ',norm_double_1
  print*, '---'
  print*,'L2 norm of mono alpha = ',norm_mono_a_2
  print*,'L2 norm of mono beta  = ',norm_mono_b_2

src/utils/memory.irp.f

@@ -114,7 +114,7 @@ subroutine print_memory_usage()
   call resident_memory(rss)
   call total_memory(mem)
 
-  write(*,'(A,F14.6,A,F14.6,A)') &
+  write(*,'(A,F14.3,A,F14.3,A)') &
     '.. >>>>> [ RES  MEM : ', rss , &
         ' GB ] [ VIRT MEM : ', mem, ' GB ] <<<<< ..'
 end