From 3a23f4933e4b47eae2a36a505949894e60083ad6 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Mon, 15 Oct 2018 11:18:45 +0200
Subject: [PATCH] FCIdump module updated from @eginer

---
 plugins/FCIdump/NEEDED_CHILDREN_MODULES      |  2 +-
 plugins/FCIdump/fcidump.irp.f                | 48 ++++++++++++--------
 scripts/make_binary.sh                       | 13 +++++-
 src/core_integrals/NEEDED_CHILDREN_MODULES   |  1 +
 src/core_integrals/core_integrals.main.irp.f |  7 +++
 src/core_integrals/core_quantities.irp.f     | 35 ++++++++++++++
 6 files changed, 85 insertions(+), 21 deletions(-)
 create mode 100644 src/core_integrals/NEEDED_CHILDREN_MODULES
 create mode 100644 src/core_integrals/core_integrals.main.irp.f
 create mode 100644 src/core_integrals/core_quantities.irp.f

diff --git a/plugins/FCIdump/NEEDED_CHILDREN_MODULES b/plugins/FCIdump/NEEDED_CHILDREN_MODULES
index 34de8ddb..491f2016 100644
--- a/plugins/FCIdump/NEEDED_CHILDREN_MODULES
+++ b/plugins/FCIdump/NEEDED_CHILDREN_MODULES
@@ -1 +1 @@
-Determinants Davidson
+Determinants DavidsonUndressed core_integrals
diff --git a/plugins/FCIdump/fcidump.irp.f b/plugins/FCIdump/fcidump.irp.f
index f93c1128..8d334fc5 100644
--- a/plugins/FCIdump/fcidump.irp.f
+++ b/plugins/FCIdump/fcidump.irp.f
@@ -1,21 +1,25 @@
 program fcidump
   implicit none
+  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 :: ii(8), jj(8), kk(8),ll(8)
+  integer :: i1,j1,k1,l1
+  integer :: i2,j2,k2,l2
   integer*8 :: m
   character*(2), allocatable :: A(:)
 
-  print *, '&FCI NORB=', mo_tot_num, ', NELEC=', elec_num, &
+  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(mo_tot_num))
+  allocate (A(n_act_orb))
   A = '1,'
-  print *, 'ORBSYM=', (A(i), i=1,mo_tot_num) 
-  print *,'ISYM=0,'
-  print *,'/'
+  write(i_unit_output,*) 'ORBSYM=', (A(i), i=1,n_act_orb) 
+  write(i_unit_output,*) 'ISYM=0,'
+  write(i_unit_output,*) '/'
   deallocate(A)
   
-  integer*8 :: i8, k1
   integer(key_kind), allocatable :: keys(:)
   double precision, allocatable  :: values(:)
   integer(cache_map_size_kind)   :: n_elements, n_elements_max
@@ -23,14 +27,18 @@ program fcidump
 
   double precision :: get_mo_bielec_integral, integral
 
-  do l=1,mo_tot_num
-   do k=1,mo_tot_num
-    do j=l,mo_tot_num
-     do i=k,mo_tot_num
-       if (i>=j) then
-          integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
+  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_mo_bielec_integral(i1,j1,k1,l1,mo_integrals_map)
           if (dabs(integral) > mo_integrals_threshold) then 
-            print *, integral, i,k,j,l
+            write(i_unit_output,*) integral, i,k,j,l
           endif
        end if
      enddo
@@ -38,13 +46,15 @@ program fcidump
    enddo
   enddo
 
-  do j=1,mo_tot_num
-   do i=j,mo_tot_num
-      integral = mo_mono_elec_integral(i,j)
+  do j=1,n_act_orb
+   j1 = list_act(j)
+   do i=j,n_act_orb
+    i1 = list_act(i)
+      integral = mo_mono_elec_integral(i1,j1) + core_fock_operator(i1,j1)
       if (dabs(integral) > mo_integrals_threshold) then 
-        print *, integral, i,j,0,0
+        write(i_unit_output,*) integral, i,j,0,0
       endif
    enddo
   enddo
-  print *, 0.d0, 0, 0, 0, 0
+  write(i_unit_output,*) core_energy, 0, 0, 0, 0
 end
diff --git a/scripts/make_binary.sh b/scripts/make_binary.sh
index d18e0d04..0028fbbd 100755
--- a/scripts/make_binary.sh
+++ b/scripts/make_binary.sh
@@ -35,7 +35,18 @@ fi
 
 
 # Build all sources
-for dir in ${QP_ROOT}/{src,ocaml}
+for dir in ${QP_ROOT}/{src}
+do
+  pushd $dir
+  ninja
+  if [[ $? -ne 0 ]]
+  then
+    echo "Error building ${dir}"
+  fi
+  popd
+done
+
+for dir in ${QP_ROOT}/{ocaml}
 do
   make -C ${dir}
   if [[ $? -ne 0 ]]
diff --git a/src/core_integrals/NEEDED_CHILDREN_MODULES b/src/core_integrals/NEEDED_CHILDREN_MODULES
new file mode 100644
index 00000000..6a4d0040
--- /dev/null
+++ b/src/core_integrals/NEEDED_CHILDREN_MODULES
@@ -0,0 +1 @@
+Integrals_Monoelec Integrals_Bielec Bitmask
diff --git a/src/core_integrals/core_integrals.main.irp.f b/src/core_integrals/core_integrals.main.irp.f
new file mode 100644
index 00000000..f5e9fd1b
--- /dev/null
+++ b/src/core_integrals/core_integrals.main.irp.f
@@ -0,0 +1,7 @@
+program core_integrals
+  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  print*,'core energy = ',core_energy 
+end
diff --git a/src/core_integrals/core_quantities.irp.f b/src/core_integrals/core_quantities.irp.f
new file mode 100644
index 00000000..a120cf78
--- /dev/null
+++ b/src/core_integrals/core_quantities.irp.f
@@ -0,0 +1,35 @@
+BEGIN_PROVIDER [double precision, core_energy]
+ implicit none
+ integer :: i,j,k,l
+ core_energy = 0.d0
+ do i = 1, n_core_orb
+  j = list_core(i)
+  core_energy += 2.d0 * mo_mono_elec_integral(j,j) + mo_bielec_integral_jj(j,j)
+  do k = i+1, n_core_orb
+   l = list_core(k)
+   core_energy += 2.d0 * (2.d0 * mo_bielec_integral_jj(j,l) - mo_bielec_integral_jj_exchange(j,l))
+  enddo
+ enddo
+ core_energy += nuclear_repulsion
+
+END_PROVIDER
+
+BEGIN_PROVIDER [double precision, core_fock_operator, (mo_tot_num,mo_tot_num)]
+ implicit none
+ integer :: i,j,k,l,m,n
+ double precision :: get_mo_bielec_integral
+ BEGIN_DOC 
+! this is the contribution to the Fock operator from the core electrons
+ END_DOC
+ core_fock_operator = 0.d0
+ do i = 1, n_act_orb
+  j = list_act(i)
+  do k = 1, n_act_orb
+   l = list_act(k)
+   do m = 1, n_core_orb
+    n = list_core(m)
+    core_fock_operator(j,l) += 2.d0 * get_mo_bielec_integral(j,n,l,n,mo_integrals_map) - get_mo_bielec_integral(j,n,n,l,mo_integrals_map)
+   enddo
+  enddo
+ enddo
+END_PROVIDER