Compute the H elements for QMC

devel/variance/rpt2_match.irp.f

@@ -88,7 +88,7 @@ program variance
     call print_summary(psi_energy_with_nucl_rep,pt2,error,var,norm,N_det,N_occ_pattern,N_states,psi_s2)
 
     do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
 
     call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)

devel/variance/variance.irp.f

@@ -88,7 +88,7 @@ program variance
     call print_summary(psi_energy_with_nucl_rep,pt2,error,var,norm,N_det,N_occ_pattern,N_states,psi_s2)
 
     do k=1,N_states
-      rpt2(:) = pt2(:)/(1.d0 + norm(k))
+      rpt2(k) = pt2(k)/(1.d0 + norm(k))
     enddo
 
     call save_iterations(psi_energy_with_nucl_rep(1:N_states),rpt2,N_det)

stable/champ/save_for_champ.irp.f

@@ -4,7 +4,8 @@ program qmcpack
 ! Generates a file for CHAMP
   END_DOC
 
-  integer :: i,j
+  integer :: i,j,k,l, istate
+
   read_wf = .True.
   TOUCH read_wf
   do j=1,ao_prim_num_max
@@ -22,4 +23,71 @@ program qmcpack
   call system('rm '//trim(ezfio_filename)//'/mo_basis/ao_md5')
   call system('$QP_ROOT/src/champ/qp_convert.py '//trim(ezfio_filename))
 
+  integer :: iunit
+  integer, external :: getUnitAndOpen
+  iunit = getUnitAndOpen(trim(ezfio_filename)//'.H','w')
+
+  double precision, external :: diag_h_mat_elem
+  write(iunit,*) N_states
+  do istate=1,N_states
+    write(iunit,*) istate, psi_energy_with_nucl_rep(istate)
+  enddo
+  write(iunit,*) N_det
+  do k=1,N_det
+    write(iunit,'(I10,X,F22.15)') k, diag_h_mat_elem(psi_det(1,1,k),N_int) + nuclear_repulsion
+  enddo
+
+
+  double precision :: F(N_states)
+  integer(bit_kind), allocatable :: det(:,:,:)
+  double precision , allocatable :: coef(:,:)
+  integer :: ispin
+  double precision :: norm(N_states), hij
+  allocate(det(N_int,2,N_det), coef(N_det,N_states))
+  do j=1,mo_num
+    do i=1,j-1
+      do ispin=1,2
+        call build_singly_excited_wavefunction(j,i,1,det,coef)
+        F = 0.d0
+        do istate=1,N_states
+          norm(istate) = 0.d0
+          do k=1,N_det
+            norm(istate) = norm(istate) + coef(k,istate) * coef(k,istate)
+          enddo
+          if (norm(istate) > 0.d0) then
+            norm(istate) = (1.d0/dsqrt(norm(istate)))
+          endif
+        enddo
+        if (sum(norm(:)) > 0.d0) then
+          
+          do istate = 1,N_states
+            coef(:,istate) = coef(:,istate) * norm(istate)
+          enddo
+          !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,istate,hij) REDUCTION(+:F)
+          do k=1,N_det
+            call i_H_j(det(1,1,k), det(1,1,k), N_int, hij)
+            do istate=1,N_states
+              F(istate) = F(istate) + hij*coef(k,istate)*coef(k,istate)
+            enddo
+            do l=1,k-1
+              call i_H_j(det(1,1,k), det(1,1,l), N_int, hij)
+              do istate=1,N_states
+                F(istate) = F(istate) + 2.d0*hij*coef(k,istate)*coef(l,istate)
+              enddo
+            enddo
+          enddo
+          !$OMP END PARALLEL DO
+          F(:) = F(:) - psi_energy(:)
+        endif
+        do istate=1,N_states
+          write(iunit,'(I4,X,I4,X,I1,X,I3,X,F22.15)') i, j, ispin, istate, F(istate) 
+        enddo
+      enddo
+    enddo
+  enddo
+
+  deallocate(det,coef)
+
+  close(iunit)
+
 end