Merged QMC modules (#208)

* Fixed mmap

* Truncated wf a la QMC=Chem

* Merged QmcChem and qmcpack modules

plugins/QMC/truncate_wf_spin.irp.f

@@ -0,0 +1,104 @@
+program e_curve
+  use bitmasks
+  implicit none
+  integer :: i,j,k, kk, nab, m, l
+  double precision :: norm, E, hij, num, ci, cj
+  integer, allocatable :: iorder(:)
+  double precision , allocatable :: norm_sort(:)
+  double precision :: e_0(N_states)
+  if (.not.read_wf) then
+    stop 'Please set read_wf to true'
+  endif
+
+  PROVIDE mo_bielec_integrals_in_map H_apply_buffer_allocated
+
+  nab = n_det_alpha_unique+n_det_beta_unique
+  allocate ( norm_sort(0:nab), iorder(0:nab) )
+
+  double precision :: thresh
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  double precision, allocatable  :: u_t(:,:), v_t(:,:), s_t(:,:)
+  double precision, allocatable  :: u_0(:,:), v_0(:,:)
+  allocate(u_t(N_states,N_det),v_t(N_states,N_det),s_t(N_states,N_det))
+  allocate(u_0(N_states,N_det),v_0(N_states,N_det))
+
+  print *,  'Threshold?'
+  read(*,*) thresh
+
+  norm_sort(0) = 0.d0
+  iorder(0) = 0
+  do i=1,n_det_alpha_unique
+   norm_sort(i) = det_alpha_norm(i)
+   iorder(i) = i
+  enddo
+  
+  do i=1,n_det_beta_unique
+   norm_sort(i+n_det_alpha_unique) = det_beta_norm(i)
+   iorder(i+n_det_alpha_unique) = -i
+  enddo
+  
+  call dsort(norm_sort(1),iorder(1),nab)
+
+
+  PROVIDE psi_bilinear_matrix_values nuclear_repulsion 
+  print *,  ''
+  do j=0,nab
+    i = iorder(j)
+    if (i<0) then
+      do k=1,n_det
+        if (psi_bilinear_matrix_columns(k) == -i) then
+          psi_bilinear_matrix_values(k,1) = 0.d0
+        endif
+      enddo
+    else
+      do k=1,n_det
+        if (psi_bilinear_matrix_rows(k) ==  i) then
+          psi_bilinear_matrix_values(k,1) = 0.d0
+        endif
+      enddo
+    endif
+    if (thresh > norm_sort(j)) then
+      cycle
+    endif
+
+    u_0 = psi_bilinear_matrix_values(1:N_det,1:N_states)
+    v_t = 0.d0
+    s_t = 0.d0
+    call dtranspose(                                                   &
+        u_0,                                                           &
+        size(u_0, 1),                                                  &
+        u_t,                                                           &
+        size(u_t, 1),                                                  &
+        N_det, N_states)
+    call H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_states,N_det,1,N_det,0,1)
+    call dtranspose(                                                   &
+        v_t,                                                           &
+        size(v_t, 1),                                                  &
+        v_0,                                                           &
+        size(v_0, 1),                                                  &
+        N_states, N_det)
+    
+    double precision, external :: u_dot_u, u_dot_v
+    do i=1,N_states
+      e_0(i) = u_dot_v(v_t(1,i),u_0(1,i),N_det)/u_dot_u(u_0(1,i),N_det)
+    enddo
+
+    m = 0
+    do k=1,n_det
+     if (psi_bilinear_matrix_values(k,1) /= 0.d0) then
+      m = m+1
+     endif
+    enddo
+
+    E = E_0(1) + nuclear_repulsion
+    norm = u_dot_u(u_0(1,1),N_det)
+    print *,  'Number of determinants:', m
+    print *,  'Energy', E
+    exit
+  enddo
+  call wf_of_psi_bilinear_matrix()
+  call save_wavefunction
+
+  deallocate (iorder, norm_sort)
+end
+

src/Determinants/spindeterminants.irp.f

@@ -1204,3 +1204,36 @@ N_int;;
 END_TEMPLATE
 
 
+subroutine wf_of_psi_bilinear_matrix(truncate)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+! Generate a wave function containing all possible products 
+! of alpha and beta determinants
+  END_DOC
+  logical, intent(in)            :: truncate
+  integer                        :: i,j,k
+  integer(bit_kind)              :: tmp_det(N_int,2)
+  integer                        :: idx
+  integer, external              :: get_index_in_psi_det_sorted_bit
+  double precision               :: norm(N_states)
+  PROVIDE psi_bilinear_matrix
+
+  do k=1,N_det
+   i = psi_bilinear_matrix_rows(k)
+   j = psi_bilinear_matrix_columns(k)
+   psi_det(1:N_int,1,k) = psi_det_alpha_unique(1:N_int,i)
+   psi_det(1:N_int,2,k) = psi_det_beta_unique (1:N_int,j)
+  enddo
+  psi_coef(1:N_det,1:N_states) = psi_bilinear_matrix_values(1:N_det,1:N_states)
+  TOUCH psi_det psi_coef 
+
+  psi_det   = psi_det_sorted
+  psi_coef  = psi_coef_sorted
+  do while (sum( dabs(psi_coef(N_det,1:N_states)) ) == 0.d0)
+    N_det -= 1
+  enddo
+  SOFT_TOUCH psi_det psi_coef N_det
+
+end
+

src/Utils/fortran_mmap.c

@@ -70,3 +70,12 @@ void munmap_fortran(size_t bytes, int fd, void* map)
     }
     close(fd);
 }
+
+
+void msync_fortran(size_t bytes, int fd, void* map)
+{
+    if (msync(map, bytes, MS_SYNC) == -1) {
+        perror("Error syncing the mmap file");
+    }
+}
+

src/Utils/map_functions.irp.f

@@ -52,18 +52,14 @@ subroutine map_save_to_disk(filename,map)
   map % consolidated_idx (map % map_size + 2_8) = k
   map % consolidated = .True.
 
+  integer*8 :: n_elements
+  n_elements = int(map % n_elements,8)
 
-  call munmap( (/ map % map_size + 2_8 /), 8, fd(1), c_pointer(1))
+  print *,  'Writing data to disk...'
-  call mmap(trim(filename)//'_consolidated_idx', (/ map % map_size + 2_8 /), 8, fd(1), .True., c_pointer(1))
+  call msync ( (/ map % map_size + 2_8 /),    8, fd(1), c_pointer(1))
-  call c_f_pointer(c_pointer(1),map % consolidated_idx, (/ map % map_size +2_8/))
+  call msync ( (/ n_elements /), cache_key_kind, fd(2), c_pointer(2))
-
+  call msync ( (/ n_elements /), integral_kind , fd(3), c_pointer(3))
-  call munmap( (/ map % n_elements /), cache_key_kind, fd(2), c_pointer(2))
+  print *,  'Done'
-  call mmap(trim(filename)//'_consolidated_key', (/ map % n_elements /), cache_key_kind, fd(2), .True., c_pointer(2))
-  call c_f_pointer(c_pointer(2),map % consolidated_key, (/ map % n_elements /))
-
-  call munmap( (/ map % n_elements /), integral_kind, fd(3), c_pointer(3))
-  call mmap(trim(filename)//'_consolidated_value', (/ map % n_elements /), integral_kind, fd(3), .True., c_pointer(3))
-  call c_f_pointer(c_pointer(3),map % consolidated_value, (/ map % n_elements /))
 
 end
 
@@ -79,8 +75,6 @@ subroutine map_load_from_disk(filename,map)
   integer*8                      :: i,k,l
   integer*4                      :: j,n_elements
 
-
-
   if (map % consolidated) then
     stop 'map already consolidated'
   endif

src/Utils/mmap.f90

@@ -15,7 +15,14 @@ module mmap_module
       integer(c_int), intent(in), value    :: read_only
     end function
 
-    subroutine c_munmap(length, fd, map) bind(c,name='munmap_fortran')
+    subroutine c_munmap_fortran(length, fd, map) bind(c,name='munmap_fortran')
+      use iso_c_binding
+      integer(c_size_t), intent(in), value :: length
+      integer(c_int), intent(in), value :: fd
+      type(c_ptr), intent(in), value    :: map
+    end subroutine
+
+    subroutine c_msync_fortran(length, fd, map) bind(c,name='msync_fortran')
       use iso_c_binding
       integer(c_size_t), intent(in), value :: length
       integer(c_int), intent(in), value :: fd
@@ -61,7 +68,23 @@ module mmap_module
 
       length = PRODUCT( shape(:) ) * bytes
       fd_ = fd
-      call c_munmap( length, fd_, map)
+      call c_munmap_fortran( length, fd_, map)
+  end subroutine
+
+  subroutine msync(shape, bytes, fd, map)
+      use iso_c_binding
+      implicit none
+      integer*8, intent(in)          :: shape(:)  ! Shape of the array to map
+      integer, intent(in)            :: bytes     ! Number of bytes per element
+      integer, intent(in)            :: fd        ! File descriptor
+      type(c_ptr), intent(in)        :: map       ! C pointer
+
+      integer(c_size_t)              :: length
+      integer(c_int)                 :: fd_
+
+      length = PRODUCT( shape(:) ) * bytes
+      fd_ = fd
+      call c_msync_fortran( length, fd_, map)
   end subroutine
 
 end module mmap_module