OpenMP lock bug

plugins/Full_CI/micro_pt2.irp.f

@@ -35,6 +35,8 @@ subroutine run_wf
   call provide_everything
   integer :: rc, i
 
+  print *,  'Contribution to PT2 running'
+
   !$OMP PARALLEL PRIVATE(i)
   i = omp_get_thread_num()
   call H_apply_FCI_PT2_slave_tcp(i)

plugins/QmcChem/e_curve_qmc.irp.f

@@ -0,0 +1,102 @@
+program e_curve
+  use bitmasks
+  implicit none
+  integer :: i,j,k, nab, m, l
+  double precision :: norm, E, hij, num, ci, cj
+  integer, allocatable :: iorder(:)
+  double precision , allocatable :: norm_sort(:)
+  nab = n_det_alpha_unique+n_det_beta_unique
+  allocate ( norm_sort(0:nab), iorder(0:nab) )
+
+
+  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)
+
+  if (.not.read_wf) then
+    stop 'Please set read_wf to true'
+  endif
+
+  PROVIDE psi_bilinear_matrix_values nuclear_repulsion 
+  print *,  ''
+  print *,  '=============================='
+  print *,  'Energies at different cut-offs'
+  print *,  '=============================='
+  print *,  ''
+  print *,  '=========================================================='
+  print '(A8,2X,A8,2X,A12,2X,A10,2X,A12)',  'Thresh.', 'Ndet', 'Cost', 'Norm', 'E'
+  print *,  '=========================================================='
+  double precision :: thresh
+  integer(bit_kind), allocatable :: det_i(:,:), det_j(:,:)
+  thresh = 1.d-10
+  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
+    num = 0.d0
+    norm = 0.d0
+    m = 0
+    !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(k,l,det_i,det_j,ci,cj,hij) REDUCTION(+:norm,m,num)
+    allocate( det_i(N_int,2), det_j(N_int,2))
+    !$OMP DO SCHEDULE(guided)
+    do k=1,n_det
+      if (psi_bilinear_matrix_values(k,1) == 0.d0) then
+        cycle
+      endif
+      ci = psi_bilinear_matrix_values(k,1)
+      det_i(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(k))
+      det_i(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(k))
+      do l=1,n_det
+        if (psi_bilinear_matrix_values(l,1) == 0.d0) then
+          cycle
+        endif
+        cj = psi_bilinear_matrix_values(l,1)
+        det_j(:,1) = psi_det_alpha_unique(:,psi_bilinear_matrix_rows(l))
+        det_j(:,2) = psi_det_beta_unique(:,psi_bilinear_matrix_columns(l))
+        call i_h_j(det_i, det_j, N_int, hij)
+        num = num + ci*cj*hij
+      enddo
+     norm = norm + ci*ci
+     m = m+1
+    enddo
+    !$OMP END DO
+    deallocate (det_i,det_j)
+    !$OMP END PARALLEL
+    if (m == 0) then
+       exit
+    endif
+    E = num / norm + nuclear_repulsion
+    print '(E9.1,2X,I8,2X,F10.2,2X,F10.8,2X,F12.6)',  thresh, m, &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * (nab-1) ) / &
+      dble( elec_alpha_num**3 + elec_alpha_num**2 * (j-1)), norm, E
+    thresh = thresh * 2.d0
+  enddo
+  print *,  '=========================================================='
+
+  deallocate (iorder, norm_sort)
+end
+

plugins/QmcChem/save_for_qmcchem.irp.f

@@ -1,9 +1,46 @@
 program save_for_qmc
- read_wf = .True.
- TOUCH read_wf
- print *,  "N_det = ", N_det
- call write_spindeterminants 
- if (do_pseudo) then
-   call write_pseudopotential
- endif
+
+  integer                        :: iunit
+  integer, external              :: get_unit_and_open
+  logical                        :: exists
+  double precision               :: e_ref
+  
+  ! Determinants
+  read_wf = .True.
+  TOUCH read_wf
+  print *,  "N_det = ", N_det
+  call write_spindeterminants
+  
+  ! Reference Energy
+  if (do_pseudo) then
+    call write_pseudopotential
+  endif
+  call system( &
+   'mkdir -p '//trim(ezfio_filename)//'/simulation ;' // &
+   'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/simulation/.version ; ' // &
+   'mkdir -p '//trim(ezfio_filename)//'/properties ;' // &
+   'cp '//trim(ezfio_filename)//'/.version '//trim(ezfio_filename)//'/properties/.version ; ' // &
+   'echo T > '//trim(ezfio_filename)//'/properties/e_loc' &
+  )
+  iunit = 13
+  open(unit=iunit,file=trim(ezfio_filename)//'/simulation/e_ref',action='write')
+  call ezfio_has_full_ci_energy_pt2(exists)
+  if (exists) then
+    call ezfio_get_full_ci_energy_pt2(e_ref)
+  else
+    call ezfio_has_full_ci_energy(exists)
+    if (exists) then
+      call ezfio_get_full_ci_energy(e_ref)
+    else
+      call ezfio_has_hartree_fock_energy(exists)
+      if (exists) then
+        call ezfio_get_hartree_fock_energy(e_ref)
+      else
+        e_ref = 0.d0
+      endif
+    endif
+  endif
+  write(iunit,*) e_ref
+  close(iunit)
+
 end

src/Determinants/H_apply.irp.f

@@ -214,8 +214,13 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
     duplicate(i) = .False.
   enddo
 
-  do i=1,N_det-1
+  found_duplicates = .False.
+  i=0
+  j=0
+  do while (i<N_det-1)
+    i = max(i+1,j)
     if (duplicate(i)) then
+      found_duplicates = .True.
       cycle
     endif
     j = i+1
@@ -239,14 +244,6 @@ subroutine remove_duplicates_in_psi_det(found_duplicates)
     enddo
   enddo
 
-  found_duplicates = .False.
-  do i=1,N_det
-    if (duplicate(i)) then
-      found_duplicates = .True.
-      exit
-    endif
-  enddo
-
   if (found_duplicates) then
     call write_bool(output_determinants,found_duplicates,'Found duplicate determinants')
     k=0

src/Determinants/H_apply.template.f

@@ -167,7 +167,8 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
   double precision               :: diag_H_mat_elem
   integer                        :: iproc
   integer                        :: jtest_vvvv
-  integer(omp_lock_kind), save   :: lck, ifirst=0
+  integer(omp_lock_kind), save   :: lck
+  integer, save                  :: ifirst=0
   if (ifirst == 0) then
 !$    call omp_init_lock(lck)
     ifirst=1
@@ -417,7 +418,8 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
   integer, allocatable           :: ia_ja_pairs(:,:,:)
   logical, allocatable           :: array_pairs(:,:)
   double precision               :: diag_H_mat_elem
-  integer(omp_lock_kind), save   :: lck, ifirst=0
+  integer(omp_lock_kind), save   :: lck
+  integer, save                  :: ifirst=0
   integer                        :: iproc
   
   integer(bit_kind)              :: key_mask(N_int, 2)
@@ -428,6 +430,11 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
   logical :: is_a_1p
   logical :: is_a_2p
   
+  if (ifirst == 0) then
+    ifirst=1
+!$    call omp_init_lock(lck)
+  endif
+  
   do k=1,N_int
     key_mask(k,1) = 0_bit_kind
     key_mask(k,2) = 0_bit_kind
@@ -439,11 +446,6 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,fock_diag_tmp,i_generato
   $check_double_excitation
   
 
-  if (ifirst == 0) then
-    ifirst=1
-!$    call omp_init_lock(lck)
-  endif
-  
   $initialization
   
   $omp_parallel