Update trexio plugin

devel/trexio/EZFIO.cfg

@@ -10,3 +10,21 @@ doc: Name of the exported TREXIO file
 interface: ezfio, ocaml, provider
 default: None
 
+[export_rdm]
+type: logical
+doc: If True, export two-body reduced density matrix
+interface: ezfio, ocaml, provider
+default: False
+
+[export_ao_ints]
+type: logical
+doc: If True, export two-electron integrals in AO basis
+interface: ezfio, ocaml, provider
+default: False
+
+[export_mo_ints]
+type: logical
+doc: If True, export two-electron integrals in MO basis
+interface: ezfio, ocaml, provider
+default: True
+

devel/trexio/LIB

@@ -1,2 +1,2 @@
--L/home/scemama/TREX/trexio/_install/lib -ltrexio
+-ltrexio
 

devel/trexio/NEED

@@ -4,3 +4,4 @@ mo_one_e_ints
 mo_two_e_ints
 ao_two_e_ints
 ao_one_e_ints
+two_body_rdm

devel/trexio/export_trexio.irp.f

@@ -1,4 +1,11 @@
 program export_trexio
+  implicit none
+  read_wf = .True.
+  SOFT_TOUCH read_wf
+  call run
+end
+
+subroutine run
   use trexio
   implicit none
   BEGIN_DOC
@@ -7,6 +14,7 @@ program export_trexio
 
   integer(8)                     :: f ! TREXIO file handle
   integer                        :: rc
+  double precision, allocatable  :: factor(:)
 
   print *, 'TREXIO file : '//trim(trexio_filename)
   print *, ''
@@ -26,6 +34,8 @@ program export_trexio
 ! Electrons
 ! ---------
 
+  print *, 'Electrons'
+
   rc = trexio_write_electron_up_num(f, elec_alpha_num)
   call check_success(rc)
 
@@ -36,6 +46,8 @@ program export_trexio
 ! Nuclei
 ! ------
 
+  print *, 'Nuclei'
+
   rc = trexio_write_nucleus_num(f, nucl_num)
   call check_success(rc)
 
@@ -48,10 +60,15 @@ program export_trexio
   rc = trexio_write_nucleus_label(f, nucl_label, 32)
   call check_success(rc)
 
+  rc = trexio_write_nucleus_repulsion(f, nuclear_repulsion)
+  call check_success(rc)
+
 
 ! Pseudo-potentials
 ! -----------------
 
+  print *, 'ECP'
+
   double precision, allocatable :: tmp_double(:,:)
   integer, allocatable          :: tmp_int(:,:)
 
@@ -90,39 +107,36 @@ program export_trexio
 ! Basis
 ! -----
 
+  print *, 'Basis'
+
+
   rc = trexio_write_basis_type(f, 'Gaussian', len('Gaussian'))
   call check_success(rc)
 
-  rc = trexio_write_basis_num(f, shell_num)
-  call check_success(rc)
-
-  rc = trexio_write_basis_nucleus_shell_num(f, nucleus_shell_num)
-  call check_success(rc)
-
-  rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
-  call check_success(rc)
-
-  rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
-  call check_success(rc)
-
   rc = trexio_write_basis_prim_num(f, prim_num)
   call check_success(rc)
 
-  rc = trexio_write_basis_shell_prim_num(f, shell_prim_num)
-  call check_success(rc)
+   rc = trexio_write_basis_shell_num(f, shell_num)
+   call check_success(rc)
 
-  double precision, allocatable :: factor(:)
-  allocate(factor(shell_num))
-  if (ao_normalized) then
-    factor(1:shell_num) = shell_normalization_factor(1:shell_num)
-  else
-    factor(1:shell_num) = 1.d0
-  endif
-  rc = trexio_write_basis_shell_factor(f, factor)
-  call check_success(rc)
-  deallocate(factor)
+   rc = trexio_write_basis_nucleus_index(f, basis_nucleus_index)
+   call check_success(rc)
 
-  rc = trexio_write_basis_shell_prim_index(f, shell_prim_index)
+   rc = trexio_write_basis_shell_ang_mom(f, shell_ang_mom)
+   call check_success(rc)
+
+   allocate(factor(shell_num))
+   if (ao_normalized) then
+     factor(1:shell_num) = shell_normalization_factor(1:shell_num)
+   else
+     factor(1:shell_num) = 1.d0
+   endif
+   rc = trexio_write_basis_shell_factor(f, factor)
+   call check_success(rc)
+
+   deallocate(factor)
+
+  rc = trexio_write_basis_shell_index(f, shell_index)
   call check_success(rc)
 
   rc = trexio_write_basis_exponent(f, prim_expo)
@@ -142,10 +156,11 @@ program export_trexio
   deallocate(factor)
 
 
-
 ! Atomic orbitals
 ! ---------------
 
+  print *, 'AOs'
+
   rc = trexio_write_ao_num(f, ao_num)
   call check_success(rc)
 
@@ -164,7 +179,6 @@ program export_trexio
   C_A(3) = 0.d0
 
   allocate(factor(ao_num))
-  print *, ao_first_of_shell
   if (ao_normalized) then
     do i=1,ao_num
       l = ao_first_of_shell(ao_shell(i))
@@ -180,6 +194,8 @@ program export_trexio
 ! One-e AO integrals
 ! ------------------
 
+  print *, 'AO integrals'
+
   rc = trexio_write_ao_1e_int_overlap(f,ao_overlap)
   call check_success(rc)
 
@@ -200,10 +216,59 @@ program export_trexio
   rc = trexio_write_ao_1e_int_core_hamiltonian(f,ao_one_e_integrals)
   call check_success(rc)
 
+! Two-e AO integrals
+! ------------------
+
+  if (export_ao_ints) then
+    PROVIDE ao_two_e_integrals_in_map
+
+    integer(8), parameter :: BUFSIZE=10000_8
+    double precision :: eri_buffer(BUFSIZE), integral
+    integer(4) :: eri_index(4,BUFSIZE)
+    integer(8) :: icount, offset
+
+    double precision, external :: get_ao_two_e_integral
+
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,ao_num
+      do k=1,ao_num
+        do j=l,ao_num
+          do i=k,ao_num
+            if (i==j .and. k<l) cycle
+            if (i<j) cycle
+            integral = get_ao_two_e_integral(i,j,k,l,ao_integrals_map)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+              call check_success(rc)
+              icount = 0_8
+              offset += icount
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_ao_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+      call check_success(rc)
+    end if
+  end if
+
 
 ! Molecular orbitals
 ! ------------------
 
+  print *, 'MOs'
+
 !  rc = trexio_write_mo_type(f, mo_label)
 !  call check_success(rc)
 
@@ -217,6 +282,8 @@ program export_trexio
 ! One-e MO integrals
 ! ------------------
 
+  print *, 'MO integrals'
+
   rc = trexio_write_mo_1e_int_kinetic(f,mo_kinetic_integrals)
   call check_success(rc)
 
@@ -230,22 +297,103 @@ program export_trexio
     rc = trexio_write_mo_1e_int_ecp_non_local(f,mo_pseudo_integrals_non_local)
     call check_success(rc)
   endif
-!
+
   rc = trexio_write_mo_1e_int_core_hamiltonian(f,one_e_dm_mo)
   call check_success(rc)
 
+! Two-e MO integrals
+! ------------------
 
-! RDM
-! ----
+  if (export_mo_ints) then
+    PROVIDE mo_two_e_integrals_in_map
 
-!  rc = trexio_write_rdm_one_e(f,one_e_dm_mo)
-!  call check_success(rc)
-!
-!  rc = trexio_write_rdm_one_e_up(f,one_e_dm_mo_alpha_average)
-!  call check_success(rc)
-!
-!  rc = trexio_write_rdm_one_e_dn(f,one_e_dm_mo_beta_average)
-!  call check_success(rc)
+    double precision, external :: mo_two_e_integral
+
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,mo_num
+      do k=1,mo_num
+        do j=l,mo_num
+          do i=k,mo_num
+            if (i==j .and. k<l) cycle
+            if (i<j) cycle
+            integral = mo_two_e_integral(i,j,k,l)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+              call check_success(rc)
+              icount = 0_8
+              offset += icount
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_mo_2e_int_eri(f, offset, icount, eri_index, eri_buffer)
+      call check_success(rc)
+    end if
+  end if
+
+
+
+! One-e RDM
+! ---------
+
+  rc = trexio_write_rdm_1e(f,one_e_dm_mo)
+  call check_success(rc)
+
+  rc = trexio_write_rdm_1e_up(f,one_e_dm_mo_alpha_average)
+  call check_success(rc)
+
+  rc = trexio_write_rdm_1e_dn(f,one_e_dm_mo_beta_average)
+  call check_success(rc)
+
+
+! Two-e RDM
+! ---------
+
+  if (export_rdm) then
+    PROVIDE two_e_dm_mo
+
+    icount = 0_8
+    offset = 0_8
+    do l=1,mo_num
+      do k=1,mo_num
+        do j=1,mo_num
+          do i=1,mo_num
+            integral = two_e_dm_mo(i,j,k,l)
+            if (integral == 0.d0) cycle
+            icount += 1_8
+            eri_buffer(icount) = integral
+            eri_index(1,icount) = i
+            eri_index(2,icount) = j
+            eri_index(3,icount) = k
+            eri_index(4,icount) = l
+            if (icount == BUFSIZE) then
+              rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
+              call check_success(rc)
+              icount = 0_8
+              offset += icount
+            end if
+          end do
+        end do
+      end do
+    end do
+
+    if (icount >= 0_8) then
+      rc = trexio_write_rdm_2e(f, offset, icount, eri_index, eri_buffer)
+      call check_success(rc)
+    end if
+  end if
 
 
 ! ------------------------------------------------------------------------------

stable/utilities/truncate_wf.irp.f

@@ -59,7 +59,7 @@ subroutine routine_s2
  do i=1,N_det
    print *,  i, real(weight_configuration(det_to_configuration(i),:)), real(sum(weight_configuration(det_to_configuration(i),:)))
  enddo
- print*, 'Min weight of the occupation pattern ?'
+ print*, 'Min weight of the configuration?'
  read(5,*) wmin
 
  ndet_max = 0