Added svdwf module

devel/cc/spatial_to_spin_MO.irp.f

@@ -96,7 +96,7 @@ subroutine get_fock_matrix_alpha(det,F)
 
   integer :: i,j,k
 
-  F(:,:) = fock_op_cshell_ref_bitmask(:,:)
+  F(:,:) = fock_operator_closed_shell_ref_bitmask(:,:)
 
 end    
 
@@ -112,7 +112,7 @@ subroutine get_fock_matrix_beta(det,F)
 
   integer :: i,j,k
 
-  F(:,:) = fock_op_cshell_ref_bitmask(:,:)
+  F(:,:) = fock_operator_closed_shell_ref_bitmask(:,:)
 
 end    
 

devel/svdwf/.gitignore

@@ -0,0 +1,59 @@
+IRPF90_temp/
+IRPF90_man/
+build.ninja
+irpf90.make
+ezfio_interface.irp.f
+irpf90_entities
+tags
+Makefile
+ao_basis
+ao_one_e_ints
+ao_two_e_erf_ints
+ao_two_e_ints
+aux_quantities
+becke_numerical_grid
+bitmask
+cis
+cisd
+cipsi
+davidson
+davidson_dressed
+davidson_undressed
+density_for_dft
+determinants
+dft_keywords
+dft_utils_in_r
+dft_utils_one_e
+dft_utils_two_body
+dressing
+dummy
+electrons
+ezfio_files
+fci
+generators_cas
+generators_full
+hartree_fock
+iterations
+kohn_sham
+kohn_sham_rs
+mo_basis
+mo_guess
+mo_one_e_ints
+mo_two_e_erf_ints
+mo_two_e_ints
+mpi
+mrpt_utils
+nuclei
+perturbation
+pseudo
+psiref_cas
+psiref_utils
+scf_utils
+selectors_cassd
+selectors_full
+selectors_utils
+single_ref_method
+slave
+tools
+utils
+zmq

devel/svdwf/NEED

@@ -0,0 +1 @@
+determinants

devel/svdwf/README.rst

@@ -0,0 +1,4 @@
+=====
+svdwf
+=====
+

devel/svdwf/random_svd.irp.f

@@ -0,0 +1,101 @@
+program svdwf
+  implicit none
+  BEGIN_DOC
+! Make the SVD of the alpha-beta wave function and print singular values.
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run()
+end
+
+subroutine run
+  implicit none
+  include 'constants.include.F'
+  double precision, allocatable :: U(:,:), V(:,:), D(:), A(:,:)
+  integer :: i, j, k, l, q, r, m, n, iter
+  double precision,allocatable    :: Z(:,:), P(:,:), Yt(:,:), UYt(:,:)
+  double precision :: r1,r2
+
+  m = n_det_alpha_unique
+  n = n_det_beta_unique
+
+  r = min(1000,n)
+
+  allocate(Z(m,r))
+
+  ! Z(m,r) = A(m,n).P(n,r)                                                                    
+  Z(:,:) = 0.d0
+  do l=1,r
+   do k=1,N_det
+     i = psi_bilinear_matrix_rows(k)
+     j = psi_bilinear_matrix_columns(k)
+     call random_number(r1)
+     call random_number(r2)
+     r1 = dsqrt(-2.d0*dlog(r1))
+     r2 = dtwo_pi*r2
+     Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * r1*dcos(r2)
+   enddo
+  enddo
+
+  ! Power iterations
+  allocate(P(n,r))
+  do iter=1,20
+   ! P(n,r) = At(n,m).Z(m,r)
+   P(:,:) = 0.d0
+   do l=1,r
+    do k=1,N_det
+      i = psi_bilinear_matrix_rows(k)
+      j = psi_bilinear_matrix_columns(k)
+      P(j,l) = P(j,l) + psi_bilinear_matrix_values(k,1) * Z(i,l)
+    enddo
+   enddo
+   Z(:,:) = 0.d0
+   do l=1,r
+    do k=1,N_det
+      i = psi_bilinear_matrix_rows(k)
+      j = psi_bilinear_matrix_columns(k)
+      Z(i,l) = Z(i,l) + psi_bilinear_matrix_values(k,1) * P(j,l)
+    enddo
+   enddo
+   ! Compute QR
+   call ortho_qr(Z,size(Z,1),m,r)
+  enddo
+
+  ! Y(r,n) = Zt(r,m).A(m,n)
+  allocate(Yt(n,r))
+  Yt(:,:) = 0.d0
+  do l=1,r
+    do k=1,N_det
+     i = psi_bilinear_matrix_rows(k)
+     j = psi_bilinear_matrix_columns(k)
+     Yt(j,l) = Yt(j,l) + Z(i,l) * psi_bilinear_matrix_values(k,1)
+   enddo
+  enddo
+
+  allocate(D(r),V(n,r), UYt(r,r))
+  call svd(Yt,size(Yt,1),V,size(V,1),D,UYt,size(UYt,1),n,r)
+  deallocate(Yt)
+
+  ! U(m,r) = Z(m,r).UY(r,r)                                                                   
+  allocate(U(m,r))
+  call dgemm('N','T',m,r,r,1.d0,Z,size(Z,1),UYt,size(UYt,1),0.d0,U,size(U,1))
+  deallocate(UYt,Z)
+
+  do i=1,r
+    print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
+    if (D(i) < 1.d-15) then
+      k = i
+      exit
+    endif
+  enddo
+  print *, 'threshold: ', 2.858 * D(k/2)
+  do i=1,m
+    print '(I6,4(X,F12.8))', i, U(i,1:4)
+  enddo
+  print *, ''
+  do i=1,n
+    print '(I6,4(X,F12.8))', i, V(i,1:4)
+  enddo
+
+  deallocate(U,D,V)
+end

devel/svdwf/svdwf.irp.f

@@ -0,0 +1,47 @@
+program svdwf
+  implicit none
+  BEGIN_DOC
+! TODO : Make the SVD of the alpha-beta wave function and print singular values.
+  END_DOC
+  read_wf = .True.
+  TOUCH read_wf
+  call run()
+end
+
+subroutine run
+  implicit none
+  double precision, allocatable :: U(:,:), Vt(:,:), D(:), A(:,:)
+  integer :: i, j, k, p, q
+
+  allocate( A (n_det_alpha_unique, n_det_beta_unique), &
+            U (n_det_alpha_unique, n_det_alpha_unique), &
+            Vt(n_det_beta_unique, n_det_beta_unique),  &
+            D(max(n_det_beta_unique,n_det_alpha_unique)) )
+
+  A = 0.D0
+  do k=1,N_det
+    i = psi_bilinear_matrix_rows(k)
+    j = psi_bilinear_matrix_columns(k)
+    A(i,j) = psi_bilinear_matrix_values(k,1)
+  enddo
+
+  call randomized_svd(A, size(A,1), &
+           U, size(U,1), D, Vt, size(Vt,1), n_det_alpha_unique, n_det_beta_unique, &
+           6,1000)
+
+  do i=1,n_det_beta_unique
+    print *, i, real(D(i)), real(D(i)**2), real(sum(D(1:i)**2))
+    if (D(i) < 1.d-15) then
+      k = i
+      exit
+    endif
+  enddo
+  print *, 'threshold: ', 2.858 * D(k/2)
+  do i=1,n_det_alpha_unique
+    print '(I6,4(X,F12.8))', i, U(i,1:4)
+  enddo
+  print *, ''
+  do i=1,n_det_beta_unique
+    print '(I6,4(X,F12.8))', i, Vt(1:4,i)
+  enddo
+end

stable/champ/qp_convert.py

@@ -2,11 +2,11 @@
 #
 # Modified from the QMCPACK interface developed by @tapplencourt and @abenali
 
-print "#QP -> CHAMP"
+print("#QP -> CHAMP")
 
-# ___           
-#  |  ._  o _|_ 
-# _|_ | | |  |_ 
+# ___
+#  |  ._  o _|_
+# _|_ | | |  |_
 #
 
 from ezfio import ezfio
@@ -19,7 +19,7 @@ ezfio.set_file(ezfio_path)
 
 do_pseudo = ezfio.get_pseudo_do_pseudo()
 if do_pseudo:
-    print "do_pseudo True"
+    print("do_pseudo True")
     from qp_path import QP_ROOT
 
     l_ele_path = os.path.join(QP_ROOT,"data","list_element.txt")
@@ -28,9 +28,9 @@ if do_pseudo:
 
     l_element_raw = data_raw.split("\n")
     l_element = [element_raw.split() for element_raw in l_element_raw]
-    d_z = dict((abr, z) for (z, abr, ele, _) in filter(lambda x: x != [], l_element) )
+    d_z = dict((abr, z) for (z, abr, ele, _) in [x for x in l_element if x != []] )
 else:
-    print "do_pseudo False"
+    print("do_pseudo False")
 
 try:
     n_det = ezfio.get_determinants_n_det()
@@ -38,13 +38,13 @@ except IOError:
     n_det = 1
 
 if n_det == 1:
-    print "multi_det False"
+    print("multi_det False")
 else:
-    print "multi_det True"
+    print("multi_det True")
 
-#             
-# |\/| o  _  _ 
-# |  | | _> (_ 
+#
+# |\/| o  _  _
+# |  | | _> (_
 #
 
 def list_to_string(l):
@@ -52,17 +52,17 @@ def list_to_string(l):
 
 
 ao_num = ezfio.get_ao_basis_ao_num()
-print "ao_num", ao_num
+print("ao_num", ao_num)
 
 mo_num = ezfio.get_mo_basis_mo_num()
-print "mo_num", mo_num
+print("mo_num", mo_num)
 
 alpha = ezfio.get_electrons_elec_alpha_num()
 beta = ezfio.get_electrons_elec_beta_num()
-print "elec_alpha_num", alpha
-print "elec_beta_num", beta
-print "elec_tot_num", alpha + beta
-print "spin_multiplicity", (alpha - beta) + 1
+print("elec_alpha_num", alpha)
+print("elec_beta_num", beta)
+print("elec_tot_num", alpha + beta)
+print("spin_multiplicity", (alpha - beta) + 1)
 
 l_label = ezfio.get_nuclei_nucl_label()
 l_charge = ezfio.get_nuclei_nucl_charge()
@@ -70,19 +70,19 @@ l_coord = ezfio.get_nuclei_nucl_coord()
 
 l_coord_str = [list_to_string(i) for i in zip(*l_coord)]
 
-print "nucl_num", len(l_label)
+print("nucl_num", len(l_label))
 
-#  _               
-# /   _   _  ._ _| 
-# \_ (_) (_) | (_| 
+#  _
+# /   _   _  ._ _|
+# \_ (_) (_) | (_|
 #
-print "Atomic coord in Bohr"
+print("Atomic coord in Bohr")
 
 for i, t in enumerate(zip(l_label, l_charge, l_coord_str)):
     t_1 = d_z[t[0]] if do_pseudo else t[1]
-    
+
     t_new = [t[0],t_1,t[2]]
-    print list_to_string(t_new)
+    print(list_to_string(t_new))
 
 #
 # Call externet process to get the sysmetry
@@ -93,14 +93,22 @@ process = subprocess.Popen(
     stdout=subprocess.PIPE)
 out, err = process.communicate()
 
-basis_raw, sym_raw, _ = out.split("\n\n\n")
+basis_raw, sym_raw, _ = out.decode().split("\n\n\n")
+basis_split = basis_raw.split('\n')
 
-#  _                 __        
-# |_)  _.  _ o  _   (_   _ _|_ 
-# |_) (_| _> | _>   __) (/_ |_ 
+#  _                 __
+# |_)  _.  _ o  _   (_   _ _|_
+# |_) (_| _> | _>   __) (/_ |_
 #
 
-basis_without_header = "\n".join(basis_raw.split("\n")[11:])
+beginning = 0
+for x in basis_split:
+  if x.startswith("Basis set"):
+    break
+  beginning += 1
+
+beginning+=2
+basis_without_header = "\n".join(basis_split[beginning:])
 
 import re
 l_basis_raw = re.split('\n\s*\n', basis_without_header)
@@ -118,37 +126,33 @@ for i, (a,b) in enumerate(zip(l_label,l_basis_raw)):
     else:
         continue
 
-print "BEGIN_BASIS_SET\n"
-print "\n\n".join(l_basis_clean)
-print "END_BASIS_SET"
+print("BEGIN_BASIS_SET\n")
+print("\n\n".join(l_basis_clean))
+print("END_BASIS_SET")
 
-#       _     
-# |\/| / \  _ 
-# |  | \_/ _> 
+#       _
+# |\/| / \  _
+# |  | \_/ _>
 #
 
 
 #
 # Function
 #
-d_gms_order ={ 0:["s"],
-     1:[ "x", "y", "z" ],
-     2:[ "xx", "yy", "zz", "xy", "xz", "yz" ],
-     3:[ "xxx", "yyy", "zzz", "xxy", "xxz", "yyx", "yyz", "zzx", "zzy", "xyz"],
-     4:[ "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy"] }
+d_gms_order = ["s",
+      "x", "y", "z" ,
+      "xx", "yy", "zz", "xy", "xz", "yz" ,
+      "xxx", "yyy", "zzz", "xxy", "xxz", "yyx", "yyz", "zzx", "zzy", "xyz",
+      "xxxx", "yyyy", "zzzz", "xxxy", "xxxz", "yyyx", "yyyz", "zzzx", "zzzy", "xxyy", "xxzz", "yyzz", "xxyz", "yyxz", "zzxy"] 
 
-def compare_gamess_style(item1, item2):
-    n1,n2 = map(len,(item1,item2))
-    assert (n1 == n2)
+def key_gamess_style(item1):
     try:
-        l = d_gms_order[n1]
-    except KeyError:
+        l = d_gms_order.index(item1)
+    except ValueError:
         return 0
-#       raise (KeyError, "We dont handle L than 4")
     else:
-        a = l.index(item1)
-        b = l.index(item2)
-        return cmp( a, b )
+        result = d_gms_order.index(item1)
+    return result
 
 def expend_sym_str(str_):
     #Expend x2 -> xx
@@ -184,20 +188,23 @@ def get_nb_permutation(str_):
     if (str_) == 's': return 1
     else: return n_orbital(len(str_))
 
+import functools
 def order_l_l_sym(l_l_sym):
     n = 1
-    iter_ = range(len(l_l_sym))
+    iter_ = list(range(len(l_l_sym)))
     for i in iter_:
         if n != 1:
             n += -1
-            continue 
+            continue
 
         l = l_l_sym[i]
         n = get_nb_permutation(l[2])
 
+        def local_key(x):
+          return key_gamess_style(x[2])
+
         l_l_sym[i:i + n] = sorted(l_l_sym[i:i + n],
-                                  key=lambda x: x[2],
-                                  cmp=compare_gamess_style)
+                                  key=local_key)
 
 
     return l_l_sym
@@ -237,21 +244,21 @@ def order_phase(mo_coef):
 def chunked(l, chunks_size):
     l_block = []
     for i in l:
-        chunks = [i[x:x + chunks_size] for x in xrange(0, len(i), chunks_size)]
+        chunks = [i[x:x + chunks_size] for x in range(0, len(i), chunks_size)]
         l_block.append(chunks)
     return l_block
 
 
 def print_mo_coef(mo_coef_block, l_l_sym):
-    print ""
-    print "BEGIN_MO"
-    print ""
+    print("")
+    print("BEGIN_MO")
+    print("")
     len_block_curent = 0
     nb_block = len(mo_coef_block[0])
     for i_block in range(0, nb_block):
         a = [i[i_block] for i in mo_coef_block]
-        r_ = range(len_block_curent, len_block_curent + len(a[0]))
-        print " ".join([str(i + 1) for i in r_])
+        r_ = list(range(len_block_curent, len_block_curent + len(a[0])))
+        print(" ".join([str(i + 1) for i in r_]))
 
         len_block_curent += len(a[0])
 
@@ -260,27 +267,27 @@ def print_mo_coef(mo_coef_block, l_l_sym):
             i_a = int(l[1]) - 1
             sym = l[2]
 
-            print l_label[i_a], sym, " ".join('%20.15e'%i
-                                              for i in a[i])
+            print(l_label[i_a], sym, " ".join('%20.15e'%i
+                                              for i in a[i]))
 
         if i_block != nb_block - 1:
-            print ""
+            print("")
         else:
-            print "END_MO"
+            print("END_MO")
 
 
 mo_coef = ezfio.get_mo_basis_mo_coef()
-mo_coef_transp = zip(*mo_coef)
+mo_coef_transp = list(zip(*mo_coef))
 mo_coef_block = chunked(mo_coef_transp, 4)
 print_mo_coef(mo_coef_block, l_l_sym_ordered)
 
-#  _                    
-# |_) _  _       _|  _  
-# |  _> (/_ |_| (_| (_) 
+#  _
+# |_) _  _       _|  _
+# |  _> (/_ |_| (_| (_)
 #
 if do_pseudo:
-    print ""
-    print "BEGIN_PSEUDO"
+    print("")
+    print("BEGIN_PSEUDO")
     klocmax = ezfio.get_pseudo_pseudo_klocmax()
     kmax = ezfio.get_pseudo_pseudo_kmax()
     lmax = ezfio.get_pseudo_pseudo_lmax()
@@ -318,21 +325,21 @@ if do_pseudo:
                 l_str.append(l_dump)
 
         str_ = "PARAMETERS FOR {0} ON ATOM {1} WITH ZCORE {2} AND LMAX {3} ARE"
-        print str_.format(a, i + 1, int(d_z[a])-int(l_charge[i]), int(len(l_str) - 1))
+        print(str_.format(a, i + 1, int(d_z[a])-int(l_charge[i]), int(len(l_str) - 1)))
 
         for i, l in enumerate(l_str):
             str_ = "FOR L= {0} COEFF N ZETA"
-            print str_.format(int(len(l_str) - i - 1))
+            print(str_.format(int(len(l_str) - i - 1)))
             for ii, ll in enumerate(l):
-                print " ", ii + 1, ll
+                print(" ", ii + 1, ll)
 
     str_ = "THE ECP RUN REMOVES {0} CORE ELECTRONS, AND THE SAME NUMBER OF PROTONS."
-    print str_.format(sum([int(d_z[a])-int(l_charge[i]) for i,a in enumerate(l_label)]))
-    print "END_PSEUDO"
+    print(str_.format(sum([int(d_z[a])-int(l_charge[i]) for i,a in enumerate(l_label)])))
+    print("END_PSEUDO")
 
-#  _         
-# | \  _ _|_ 
-# |_/ (/_ |_ 
+#  _
+# | \  _ _|_
+# |_/ (/_ |_
 #
 
 
@@ -343,12 +350,12 @@ bit_kind = ezfio.get_determinants_bit_kind()
 
 nexcitedstate = ezfio.get_determinants_n_states()
 
-print ""
-print "BEGIN_DET"
-print ""
-print "mo_num", mo_num
-print "det_num", n_det
-print ""
+print("")
+print("BEGIN_DET")
+print("")
+print("mo_num", mo_num)
+print("det_num", n_det)
+print("")
 
 if "QP_STATE" in os.environ:
   state = int(os.environ["QP_STATE"])-1
@@ -372,13 +379,13 @@ MultiDetAlpha = []
 MultiDetBeta = []
 for coef, (det_a, det_b) in zip(psi_coef_small, psi_det):
 
-        print coef
+        print(coef)
         MyDetA=decode(det_a)
         MyDetB=decode(det_b)
-        print MyDetA 
-        print MyDetB
-        print ''
+        print(MyDetA)
+        print(MyDetB)
+        print('')
         MultiDetAlpha.append( det_a  )
-	MultiDetBeta.append( det_b  )
-print "END_DET"
+        MultiDetBeta.append( det_b  )
+print("END_DET")
 

stable/champ/save_for_champ.irp.f

@@ -23,73 +23,73 @@ 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
-            if (sum(coef(k,:)*coef(k,:))==0.d0) cycle
-            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
-              if (sum(coef(l,:)*coef(l,:))==0.d0) cycle
-              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)
+!  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
+!            if (sum(coef(k,:)*coef(k,:))==0.d0) cycle
+!            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
+!              if (sum(coef(l,:)*coef(l,:))==0.d0) cycle
+!              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