Examples

1. Accessing sparse quantities
- 1.1. Fortran

1 Accessing sparse quantities

1.1 Fortran

program print_energy
  use trexio
  implicit none

  character*(128)  :: filename   ! Name of the input file
  integer          :: rc         ! Return code for error checking
  integer(8)       :: f          ! TREXIO file handle
  character*(128)  :: err_msg    ! Error message

This program computes the energy as:

\[ E = E_{\text{NN}} + \sum_{ij} D_{ij}\, \langle i | h | j \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle i j | k l \rangle\; \textrm{ with } \; 0 < i,j,k,l \le n \]

One needs to read from the TREXIO file:

\(n\): The number of molecular orbitals
\(E_{\text{NN}}\): The nuclear repulsion energy
\(D_{ij}\): The one-body reduced density matrix
\(\langle i |h| j \rangle\): The one-electron Hamiltonian integrals
\(\Gamma_{ijkl}\): The two-body reduced density matrix
\(\langle i j | k l \rangle\): The electron repulsion integrals

integer                       :: n
double precision              :: E, E_nn
double precision, allocatable :: D(:,:), h0(:,:)
double precision, allocatable :: G(:,:,:,:), W(:,:,:,:)

1.1.1 Declare Temporary variables

integer                       :: i, j, k, l, m
integer(8), parameter         :: BUFSIZE = 100000_8
integer(8)                    :: offset, icount, size_max
integer                       :: buffer_index(4,BUFSIZE)
double precision              :: buffer_values(BUFSIZE)

double precision, external    :: ddot   ! BLAS dot product

1.1.2 Obtain the name of the TREXIO file from the command line, and open it for reading

call getarg(1, filename)

f = trexio_open (filename, 'r', TREXIO_HDF5, rc)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error opening TREXIO file: '//trim(err_msg)
   stop
end if

1.1.3 Read the nuclear repulsion energy

rc = trexio_read_nucleus_repulsion(f, E_nn)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error reading nuclear repulsion: '//trim(err_msg)
   stop
end if

1.1.4 Read the number of molecular orbitals

rc = trexio_read_mo_num(f, n)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error reading number of MOs: '//trim(err_msg)
   stop
end if

1.1.5 Allocate memory

allocate( D(n,n), h0(n,n) )
allocate( G(n,n,n,n), W(n,n,n,n) )
G(:,:,:,:) = 0.d0
W(:,:,:,:) = 0.d0

1.1.6 Read one-electron quantities

rc = trexio_has_mo_1e_int_core_hamiltonian(f)
if (rc /= TREXIO_SUCCESS) then
   stop 'No core hamiltonian in file'
end if

rc = trexio_read_mo_1e_int_core_hamiltonian(f, h0)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error reading core Hamiltonian: '//trim(err_msg)
   stop
end if


rc = trexio_has_rdm_1e(f)
if (rc /= TREXIO_SUCCESS) then
   stop 'No 1e RDM in file'
end if

rc = trexio_read_rdm_1e(f, D)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error reading one-body RDM: '//trim(err_msg)
   stop
end if

1.1.7 Read two-electron quantities

Reading is done with OpenMP. Each thread reads its own buffer, and the buffers are then processed in parallel.

Reading the file requires a lock, so it is done in a critical section. The offset variable is shared, and it is incremented in the critical section. For each read, the function returns in icount the number of read integrals, so this variable needs also to be protected in the critical section when modified.

1.1.7.1 Electron repulsion integrals

rc = trexio_has_mo_2e_int_eri(f)
if (rc /= TREXIO_SUCCESS) then
   stop 'No electron repulsion integrals in file'
end if

rc = trexio_read_mo_2e_int_eri_size (f, size_max)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error reading number of ERIs: '//trim(err_msg)
   stop
end if

offset = 0_8
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(icount, i, j, k, l, &
!$OMP   buffer_index, buffer_values, m)
icount = BUFSIZE
do while (icount == BUFSIZE)
  !$OMP CRITICAL
  if (offset < size_max) then
    rc = trexio_read_mo_2e_int_eri(f, offset, icount, buffer_index, buffer_values)
    offset = offset + icount
  else
    icount = 0
  end if
  !$OMP END CRITICAL
  do m=1,icount
    i = buffer_index(1,m)
    j = buffer_index(2,m)
    k = buffer_index(3,m)
    l = buffer_index(4,m)
    W(i,j,k,l) = buffer_values(m)
    W(k,j,i,l) = buffer_values(m)
    W(i,l,k,j) = buffer_values(m)
    W(k,l,i,j) = buffer_values(m)
    W(j,i,l,k) = buffer_values(m)
    W(j,k,l,i) = buffer_values(m)
    W(l,i,j,k) = buffer_values(m)
    W(l,k,j,i) = buffer_values(m)
  end do
end do
!$OMP END PARALLEL

1.1.7.2 Reduced density matrix

rc = trexio_has_rdm_2e(f)
if (rc /= TREXIO_SUCCESS) then
   stop 'No two-body density matrix in file'
end if

rc = trexio_read_rdm_2e_size (f, size_max)
if (rc /= TREXIO_SUCCESS) then
   call trexio_string_of_error(rc, err_msg)
   print *, 'Error reading number of 2-RDM elements: '//trim(err_msg)
   stop
end if

offset = 0_8
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(icount, i, j, k, l, &
!$OMP   buffer_index, buffer_values, m)
icount = bufsize
do while (offset < size_max)
  !$OMP CRITICAL
  if (offset < size_max) then
    rc = trexio_read_rdm_2e(f, offset, icount, buffer_index, buffer_values)
    offset = offset + icount
  else
    icount = 0
  end if
  !$OMP END CRITICAL
  do m=1,icount
    i = buffer_index(1,m)
    j = buffer_index(2,m)
    k = buffer_index(3,m)
    l = buffer_index(4,m)
    G(i,j,k,l) = buffer_values(m)
  end do
end do
!$OMP END PARALLEL

1.1.8 Compute the energy

As \((n,m)\) 2D arrays are stored in memory as \((\n times m)\) 1D arrays, we could pass the matrices to the ddot BLAS function to perform the summations in a single call for the 1-electron quantities. Instead, we prefer to interleave the 1-electron (negative) and 2-electron (positive) summations to have a better cancellation of numerical errors.

Here \(n^4\) can be larger than the largest possible 32-bit integer, so it is not safe to pass \(n^4\) to the ddot BLAS function. Hence, we perform \(n^2\) loops, using vectors of size \(n^2\).

E = 0.d0
do l=1,n
  E = E + ddot( n, D(1,l), 1, h0(1,l), 1 ) 
  do k=1,n
     E = E + 0.5d0 * ddot( n*n, G(1,1,k,l), 1, W(1,1,k,l),  1 )
  end do
end do
E = E + E_nn

print *, 'Energy: ', E

1.1.9 Terminate

  deallocate( D, h0, G, W )

end program