Merge remote-tracking branch 'origin/dev-stable' into dev-stable

src/cipsi/selection.irp.f

@@ -178,7 +178,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
   integer(bit_kind), allocatable        :: minilist(:, :, :), fullminilist(:, :, :)
   logical, allocatable                  :: banned(:,:,:), bannedOrb(:,:)
   double precision, allocatable         :: coef_fullminilist_rev(:,:)
-  double precision, allocatable         :: mat(:,:,:), hij_cache(:,:,:)
+  double precision, allocatable         :: mat(:,:,:)
 
 
   PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
@@ -205,7 +205,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
   ! Removed to avoid introducing determinants already presents in the wf
   !double precision, parameter :: norm_thr = 1.d-16
 
-  allocate (indices(N_det), hij_cache(mo_num,mo_num,2),          &
+  allocate (indices(N_det), &
       exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
 
   ! Pre-compute excitation degrees wrt alpha determinants
@@ -511,15 +511,10 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
 
           maskInd = maskInd + 1
           if(mod(maskInd, csubset) == (subset-1)) then
-            call get_mo_two_e_integrals_ij(h2,h1,mo_num,hij_cache(1,1,1),mo_integrals_map)
-            if (sp /= 3) then  ! AA or BB
-              call get_mo_two_e_integrals_ij(h1,h2,mo_num,hij_cache(1,1,2),mo_integrals_map)
-            endif
-
             call spot_isinwf(mask, fullminilist, i_generator, fullinteresting(0), banned, fullMatch, fullinteresting)
             if(fullMatch) cycle
 
-            call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting, hij_cache)
+            call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
 
             call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
           end if
@@ -535,7 +530,7 @@ subroutine select_singles_and_doubles(i_generator, hole_mask, particle_mask, foc
     enddo
   enddo
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
-  deallocate(banned, bannedOrb, mat, hij_cache)
+  deallocate(banned, bannedOrb, mat)
 end subroutine
 
 BEGIN_TEMPLATE
@@ -765,6 +760,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_
       enddo
 
       do_diag = sum(dabs(coef)) > 0.001d0 .and. N_states > 1
+      !do_diag = .False.
 
       double precision :: eigvalues(N_states+1)
       double precision :: work(1+6*(N_states+1)+2*(N_states+1)**2)
@@ -918,7 +914,7 @@ subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_
 
 END_TEMPLATE
 
-subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting, hij_cache)
+subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
   use bitmasks
   implicit none
   BEGIN_DOC
@@ -930,7 +926,6 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
   integer, intent(in)            :: sp, i_gen, N_sel
   integer, intent(in)            :: interesting(0:N_sel)
   integer(bit_kind),intent(in)   :: mask(N_int, 2), det(N_int, 2, N_sel)
-  double precision, intent(in)   :: hij_cache(mo_num, mo_num, 2)
   logical, intent(inout)         :: bannedOrb(mo_num, 2), banned(mo_num, mo_num, 2)
   double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
 
@@ -1000,9 +995,9 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
         if(nt == 4) then
           call get_d2(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         else if(nt == 3) then
-          call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache)
+          call get_d1(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         else
-          call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)), hij_cache)
+          call get_d0(det(1,1,i), phasemask, bannedOrb, banned, mat, mask, h, p, sp, psi_selectors_coef_transp(1, interesting(i)))
         end if
     else if(nt == 4) then
         call bitstring_to_list_in_selection(mobMask(1,1), p(1,1), p(0,1), N_int)
@@ -1203,7 +1198,7 @@ subroutine get_d2(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
 end
 
 
-subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache)
+subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   use bitmasks
   implicit none
 
@@ -1214,7 +1209,6 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
   double precision, intent(in)   :: coefs(N_states)
   double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
   integer, intent(in)            :: h(0:2,2), p(0:4,2), sp
-  double precision, intent(in)   :: hij_cache(mo_num, mo_num, 2)
 
   double precision, external     :: get_phase_bi, mo_two_e_integral
   logical                        :: ok
@@ -1256,11 +1250,13 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
     p1 = p(1,ma)
     p2 = p(2,ma)
     if(.not. bannedOrb(puti, mi)) then
+      call get_mo_two_e_integrals(hfix,p1,p2,mo_num,hij_cache1(1,1),mo_integrals_map)
+      call get_mo_two_e_integrals(hfix,p2,p1,mo_num,hij_cache1(1,2),mo_integrals_map)
       tmp_row = 0d0
       do putj=1, hfix-1
         if(lbanned(putj, ma)) cycle
         if(banned(putj, puti,bant)) cycle
-        hij = hij_cache(hfix,putj,1) - hij_cache(putj,hfix,1)
+        hij = hij_cache1(putj,1) - hij_cache1(putj,2)
         if (hij /= 0.d0) then
           hij = hij * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
@@ -1272,7 +1268,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
       do putj=hfix+1, mo_num
         if(lbanned(putj, ma)) cycle
         if(banned(putj, puti,bant)) cycle
-        hij = hij_cache(putj,hfix,1) - hij_cache(hfix,putj,1)
+        hij = hij_cache1(putj,2) - hij_cache1(putj,1)
         if (hij /= 0.d0) then
           hij = hij * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
           !DIR$ LOOP COUNT AVG(4)
@@ -1463,7 +1459,7 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
 
 
 
-subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs, hij_cache)
+subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
   use bitmasks
   implicit none
 
@@ -1474,7 +1470,6 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
   double precision, intent(in) :: coefs(N_states)
   double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
   integer, intent(in) :: h(0:2,2), p(0:4,2), sp
-  double precision, intent(in)   :: hij_cache(mo_num, mo_num, 2)
 
   integer :: i, j, k, s, h1, h2, p1, p2, puti, putj
   double precision :: hij, phase
@@ -1483,13 +1478,17 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
 
   integer, parameter :: bant=1
 
+  double precision, allocatable :: hij_cache1(:), hij_cache2(:)
+  allocate (hij_cache1(mo_num),hij_cache2(mo_num))
+
   PROVIDE mo_integrals_threshold
 
   if(sp == 3) then ! AB
     h1 = p(1,1)
     h2 = p(1,2)
     do p1=1, mo_num
       if(bannedOrb(p1, 1)) cycle
+      call get_mo_two_e_integrals(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map)
       do p2=1, mo_num
         if(bannedOrb(p2,2)) cycle
         if(banned(p1, p2, bant)) cycle ! rentable?
@@ -1498,7 +1497,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
           call i_h_j(gen, det, N_int, hij)
         else
           phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
-          hij = hij_cache(p2,p1,1) * phase
+          hij = hij_cache1(p2) * phase
         end if
         if (dabs(hij) < mo_integrals_threshold) cycle
         !DIR$ LOOP COUNT AVG(4)
@@ -1513,6 +1512,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
     p2 = p(2,sp)
     do puti=1, mo_num
       if (bannedOrb(puti, sp)) cycle
+      call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map)
+      call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map)
       do putj=puti+1, mo_num
         if(bannedOrb(putj, sp)) cycle
         if(banned(puti, putj, bant)) cycle ! rentable?
@@ -1521,7 +1522,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
           call i_h_j(gen, det, N_int, hij)
           if (dabs(hij) < mo_integrals_threshold) cycle
         else
-          hij = hij_cache(putj,puti,1) - hij_cache(putj,puti,2)
+          hij = hij_cache1(putj) - hij_cache2(putj)
           if (dabs(hij) < mo_integrals_threshold) cycle
           hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
         end if
@@ -1533,6 +1534,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs,
     end do
   end if
 
+  deallocate(hij_cache1,hij_cache2)
 end
 
 

src/cipsi_utils/pt2_stoch_routines.irp.f

@@ -894,9 +894,8 @@ integer function pt2_find_sample_lr(v, w, l_in, r_in)
    do t=1, pt2_N_teeth
      tooth_width = tilde_cW(pt2_n_0(t+1)) - tilde_cW(pt2_n_0(t))
      if (tooth_width == 0.d0) then
-       tooth_width = sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1)))
+       tooth_width = max(1.d-15,sum(tilde_w(pt2_n_0(t):pt2_n_0(t+1))))
      endif
-     ASSERT(tooth_width > 0.d0)
      do i=pt2_n_0(t)+1, pt2_n_0(t+1)
        pt2_w(i) = tilde_w(i) * pt2_W_T / tooth_width
      end do

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -70,6 +70,9 @@ subroutine mo_two_e_integrals_index(i,j,k,l,i1)
     else
       call add_integrals_to_map(full_ijkl_bitmask_4)
     endif
+
+    integer*8                      :: get_mo_map_size, mo_map_size
+    mo_map_size = get_mo_map_size()
     double precision, external     :: map_mb
     print*,'Molecular integrals provided:'
     print*,' Size of MO map           ', map_mb(mo_integrals_map) ,'MB'
@@ -79,9 +82,6 @@ subroutine mo_two_e_integrals_index(i,j,k,l,i1)
   call wall_time(wall_2)
   call cpu_time(cpu_2)
 
-  integer*8                      :: get_mo_map_size, mo_map_size
-  mo_map_size = get_mo_map_size()
-
   print*,' cpu  time :',cpu_2 - cpu_1, 's'
   print*,' wall time :',wall_2 - wall_1, 's  ( x ', (cpu_2-cpu_1)/(wall_2-wall_1), ')'
 

src/utils/mmap.f90

@@ -2,6 +2,8 @@ module mmap_module
 
   use iso_c_binding
 
+  character*(256)   :: mmap_prefix = '/tmp/'
+
   type mmap_type
     type(c_ptr)     :: ptr          ! Pointer to the data
     character*(128) :: filename     ! Name of the file
@@ -155,7 +157,15 @@ subroutine mmap_create(filename, shape, bytes, read_only, single_node, map)
         map%filename = filename
       else
         call getenv('EZFIO_FILE', map%filename)
-        map%filename = trim(map%filename) // '/work/tmpfile'
+        if (trim(map%filename) /= '') then
+          map%filename = trim(map%filename) // '/work/'
+        else
+          call getenv('TMPDIR', map%filename)
+          if (trim(map%filename) == '') then
+            map%filename = '/tmp/'
+          endif
+        endif
+        map%filename = trim(map%filename) // '/tmpfile'
       endif
 
       map%length = int(bytes,8)