From 233fa3e9bc0f8c103cd2f47edc6b44ff04cfb013 Mon Sep 17 00:00:00 2001
From: Emmanuel Giner <giner.emmanuel@gmail.com>
Date: Thu, 1 Jul 2021 16:14:58 +0200
Subject: [PATCH] added  fcidump_pyscf.irp.f

---
 src/tools/fcidump_pyscf.irp.f | 78 +++++++++++++++++++++++++++++++++++
 1 file changed, 78 insertions(+)
 create mode 100644 src/tools/fcidump_pyscf.irp.f

diff --git a/src/tools/fcidump_pyscf.irp.f b/src/tools/fcidump_pyscf.irp.f
new file mode 100644
index 00000000..9cbf733a
--- /dev/null
+++ b/src/tools/fcidump_pyscf.irp.f
@@ -0,0 +1,78 @@
+program fcidump
+  implicit none
+  BEGIN_DOC
+! Produce a regular `FCIDUMP` file from the |MOs| stored in the |EZFIO|
+! directory.
+!
+! To specify an active space, the class of the |MOs| have to set in the
+! |EZFIO| directory (see :ref:`qp_set_mo_class`).
+!
+! The :ref:`fcidump` program supports 3 types of |MO| classes :
+!
+! * the *core* orbitals which are always doubly occupied in the
+!   calculation
+!
+! * the *deleted* orbitals that are never occupied in the calculation
+!
+! * the *active* orbitals that are occupied with a varying number of
+!   electrons
+!
+  END_DOC
+  character*(128) :: output
+  integer :: i_unit_output,getUnitAndOpen
+  output=trim(ezfio_filename)//'.FCIDUMP'
+  i_unit_output = getUnitAndOpen(output,'w')
+
+  integer :: i,j,k,l
+  integer :: i1,j1,k1,l1
+  integer :: i2,j2,k2,l2
+  integer*8 :: m
+  character*(2), allocatable :: A(:)
+
+  write(i_unit_output,*) '&FCI NORB=', n_act_orb, ', NELEC=', elec_num-n_core_orb*2, &
+   ', MS2=', (elec_alpha_num-elec_beta_num), ','
+  allocate (A(n_act_orb))
+  A = '1,'
+  write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb)
+  write(i_unit_output,*) 'ISYM=0,'
+  write(i_unit_output,*) '&end'
+  deallocate(A)
+
+  integer(key_kind), allocatable :: keys(:)
+  double precision, allocatable  :: values(:)
+  integer(cache_map_size_kind)   :: n_elements, n_elements_max
+  PROVIDE mo_two_e_integrals_in_map
+
+  double precision :: get_two_e_integral, integral
+
+  do l=1,n_act_orb
+   l1 = list_act(l)
+   do k=1,n_act_orb
+    k1 = list_act(k)
+    do j=l,n_act_orb
+     j1 = list_act(j)
+     do i=k,n_act_orb
+      i1 = list_act(i)
+       if (i1>=j1) then
+          integral = get_two_e_integral(i1,j1,k1,l1,mo_integrals_map)
+          if (dabs(integral) > mo_integrals_threshold) then
+            write(i_unit_output,*) integral, i,k,j,l
+          endif
+       end if
+     enddo
+    enddo
+   enddo
+  enddo
+
+  do j=1,n_act_orb
+   j1 = list_act(j)
+   do i=j,n_act_orb
+    i1 = list_act(i)
+      integral = mo_one_e_integrals(i1,j1) + core_fock_operator(i1,j1)
+      if (dabs(integral) > mo_integrals_threshold) then
+        write(i_unit_output,*) integral, i,j,0,0
+      endif
+   enddo
+  enddo
+  write(i_unit_output,*) core_energy, 0, 0, 0, 0
+end