UP | HOME

Examples

Table of Contents

1. Writing nuclear coordinates

This example demonstrates how to write the nuclear coordinates of a water molecule to a file using TREXIO. It covers the basic steps involved in opening a file, writing data, and closing it, as well as the necessary TREXIO functions to perform these actions.

1.1. C

#include <stdio.h>
#include <trexio.h>

int main() {
  int num = 3;  // Number of atoms
  double coord[][3] = {
    // xyz coordinates in atomic units
    0.  ,  0.        , -0.24962655,
    0.  ,  2.70519714,  1.85136466,
    0.  , -2.70519714,  1.85136466 };

  trexio_exit_code rc;

  // Open the TREXIO file
  trexio_t* f = trexio_open("water.trexio", 'w', TREXIO_HDF5, &rc);
  if (rc != TREXIO_SUCCESS) {
    fprintf(stderr, "Error: %s\n", trexio_string_of_error(rc));
    return -1;
  }

  // Write the number of nuclei
  rc = trexio_write_nucleus_num (f, num);
  if (rc != TREXIO_SUCCESS) {
    fprintf(stderr, "Error: %s\n", trexio_string_of_error(rc));
    return -1;
  }

  // Write the nuclear coordinates
  rc = trexio_write_nucleus_coord (f, &coord[0][0]);
  if (rc != TREXIO_SUCCESS) {
    fprintf(stderr, "Error: %s\n", trexio_string_of_error(rc));
    return -1;
  }

  // Close the TREXIO file
  rc = trexio_close(f);
  if (rc != TREXIO_SUCCESS) {
    fprintf(stderr, "Error: %s\n", trexio_string_of_error(rc));
    return -1;
  }
  return 0;
}

1.2. Python

This code uses the TREXIO Python binding to create a new TREXIO file named water.trexio, and write the nuclear coordinates of a water molecule.

The coord variable is a list of three lists, each containing the x, y, and z coordinates of the water molecule's nuclei.

The with statement is used to ensure the file is properly closed after the write is complete.

The trexio.write_nucleus_num function is used to write the number of nuclei in the system.

The trexio.write_nucleus_coord function is used to write the nuclear coordinates of the system. 

import trexio
coord = [    # xyz coordinates in atomic units
    [0. , 0., -0.24962655],
    [0. , 2.70519714, 1.85136466],
    [0. , -2.70519714, 1.85136466]
]
# The Python API calls can raise `trexio.Error`
# exceptions to be handled via try/except clauses
# in the user application
with trexio.File("water.trexio", 'w',
                 back_end=trexio.TREXIO_HDF5) as f:
    trexio.write_nucleus_num(f, len(coord))
    trexio.write_nucleus_coord(f, coord)

1.3. Fortran

In Fortran, you will need to use the trexio module. Due to potential issues with compiler versions and .mod files, it is recommended to include the source version of this module directly in your source code. The corresponding file is installed in a standard location when the library is installed using make install, in the include directory.

To avoid duplicating this file every time you update the TREXIO library, create a file named trexio_module.F90 (note the capital F in the suffix) containing: 

#include <trexio_f.f90>

Compile this new file and add the trexio_module.o file at link time. 

program trexio_water
  use trexio

  integer, parameter        :: num=3       ! Number of nuclei
  double precision          :: coord(3,3)  ! Array of atom coordinates

  integer(trexio_t)         :: f        ! The TREXIO file handle
  integer(trexio_exit_code) :: rc       ! TREXIO return code
  character*(128)           :: err_msg  ! String holding the error message

  coord(:,:) = reshape( (/ 0.d0  ,  0.d0        , -0.24962655d0, &
                           0.d0  ,  2.70519714d0,  1.85136466d0, &
                           0.d0  , -2.70519714d0,  1.85136466d0  /), &
                         shape(coord) )

  ! Open the TREXIO file
  f = trexio_open ('water.trexio', 'w', TREXIO_HDF5, rc)
  if (rc /= TREXIO_SUCCESS) then
     call trexio_string_of_error(rc, err_msg)
     print *, 'Error: '//trim(err_msg)
     call exit(-1)
  end if

  ! Write the number of nuclei
  rc = trexio_write_nucleus_num (f, num)
  if (rc /= TREXIO_SUCCESS) then
     call trexio_string_of_error(rc, err_msg)
     print *, 'Error: '//trim(err_msg)
     call exit(-1)
  end if

  ! Write the nuclear coordinates
  rc = trexio_write_nucleus_coord (f, coord)
  if (rc /= TREXIO_SUCCESS) then
     call trexio_string_of_error(rc, err_msg)
     print *, 'Error: '//trim(err_msg)
     call exit(-1)
  end if

  ! Close the TREXIO file
  rc = trexio_close(f)
  if (rc /= TREXIO_SUCCESS) then
     call trexio_string_of_error(rc, err_msg)
     print *, 'Error: '//trim(err_msg)
     call exit(-1)
  end if

end program

2. Accessing sparse quantities (integrals)

2.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} \gamma_{ij}\, \langle j | h | i \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle k l | i j \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
\(\gamma_{ij}\)
The one-body reduced density matrix
\(\langle j |h| i \rangle\)
The one-electron Hamiltonian integrals
\(\Gamma_{ijkl}\)
The two-body reduced density matrix
\(\langle k l | i j \rangle\)
The electron repulsion integrals
integer                       :: n
double precision              :: E, E_nn
double precision, allocatable :: D(:,:), h0(:,:)
double precision, allocatable :: G(:,:,:,:), W(:,:,:,:)

2.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

2.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_AUTO, 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

2.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

2.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

2.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

2.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

2.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.

2.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
2.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

2.1.8. Compute the energy

When the orbitals are real, we can use

\begin{eqnarray*} E &=& E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle j | h | i \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle k l | i j \rangle \\ &=& E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle i | h | j \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle i j | k l \rangle \\ \end{eqnarray*}

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

2.1.9. Terminate

  deallocate( D, h0, G, W )

end program

2.2. Python

import sys
import trexio
import numpy as np

BUFSIZE = 100000

This program computes the energy as:

\[ E = E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle j | h | i \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle k l | i j \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
\(\gamma_{ij}\)
The one-body reduced density matrix
\(\langle j |h| i \rangle\)
The one-electron Hamiltonian integrals
\(\Gamma_{ijkl}\)
The two-body reduced density matrix
\(\langle k l | i j \rangle\)
The electron repulsion integrals

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

filename = sys.argv[1]
f = trexio.File(filename, 'r', trexio.TREXIO_AUTO)

2.2.2. Read the nuclear repulsion energy

E_nn = trexio.read_nucleus_repulsion(f)

2.2.3. Read the number of molecular orbitals

n = trexio.read_mo_num(f)

2.2.4. Read one-electron quantities

if not trexio.has_mo_1e_int_core_hamiltonian(f):
  print("No core hamiltonian in file")
  sys.exit(-1)

h0 = trexio.read_mo_1e_int_core_hamiltonian(f)

if not trexio.has_rdm_1e(f):
  print("No 1e RDM in file")
  sys.exit(-1)

D = trexio.read_rdm_1e(f)

2.2.5. Read two-electron quantities

2.2.5.1. Electron repulsion integrals
if not trexio.has_mo_2e_int_eri(f):
  print("No electron repulsion integrals in file")
  sys.exit(-1)

size_max = trexio.read_mo_2e_int_eri_size(f)

offset = 0
icount = BUFSIZE
feof = False
W = np.zeros( (n,n,n,n) )
while not feof:
  buffer_index, buffer_values, icount, feof = trexio.read_mo_2e_int_eri(f, offset, icount)
  offset += icount
  for m in range(icount):
    i, j, k, l = buffer_index[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]
2.2.5.2. Reduced density matrix
if not trexio.has_rdm_2e(f):
  print("No two-body density matrix in file")

offset = 0
icount = BUFSIZE
feof = False
G = np.zeros( (n,n,n,n) )
while not feof:
  buffer_index, buffer_values, icount, feof = trexio.read_rdm_2e(f, offset, icount)
  offset += icount
  for m in range(icount):
    i, j, k, l = buffer_index[m]
    G[i,j,k,l] = buffer_values[m]

2.2.6. Compute the energy

When the orbitals are real, we can use

\begin{eqnarray*} E &=& E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle j | h | i \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle k l | i j \rangle \\ &=& E_{\text{NN}} + \sum_{ij} \gamma_{ij}\, \langle i | h | j \rangle\, +\, \frac{1}{2} \sum_{ijkl} \Gamma_{ijkl}\, \langle i j | k l \rangle \\ \end{eqnarray*} 
G = np.reshape(G, (n*n, n*n) )
W = np.reshape(W, (n*n, n*n) )
E = E_nn
E += 0.5*sum( [ np.dot(G[:,l], W[:,l]) for l in range(n*n) ] )
E += sum( [ np.dot(D[:,l], h0[:,l]) for l in range(n) ] )

print (f"Energy: {E}")

3. Reading determinants

3.1. Fortran

program test

  use trexio
  implicit none

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


  integer*8, allocatable :: buffer(:,:,:)
  integer(8) :: offset, icount, BUFSIZE
  integer :: ndet, int64_num, m

  integer :: occ_num_up, occ_num_dn
  integer, allocatable :: orb_list_up(:), orb_list_dn(:)

  call getarg(1, filename)

  trex_determinant_file = trexio_open(filename, 'r', TREXIO_AUTO, 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

  rc = trexio_read_determinant_num(trex_determinant_file, ndet)
  if (rc /= TREXIO_SUCCESS) then
     call trexio_string_of_error(rc, err_msg)
     print *, 'Error reading determinant_num: '//trim(err_msg)
     stop
  end if
  print *, 'ndet', ndet

  rc = trexio_get_int64_num(trex_determinant_file, int64_num)
  if (rc /= TREXIO_SUCCESS) then
     call trexio_string_of_error(rc, err_msg)
     print *, 'Error reading int64_num: '//trim(err_msg)
     stop
  end if
  print *, 'int64_num', int64_num

  BUFSIZE = 1000_8
  allocate(buffer(int64_num, 2, BUFSIZE))
  allocate(orb_list_up(int64_num*64), orb_list_dn(int64_num*64))

  offset = 0_8
  icount = BUFSIZE
  do while (icount == BUFSIZE)
    if (offset < ndet) then
      rc = trexio_read_determinant_list(trex_determinant_file, offset, icount, buffer)
      offset = offset + icount
    else
      icount = 0
    end if
    print *, '---'
    do m=1,icount
      rc = trexio_to_orbital_list_up_dn(int64_num, buffer(1,1,m), &
           orb_list_up, orb_list_dn, occ_num_up, occ_num_dn)
      print '(100(I3,X))', (orb_list_up(1:occ_num_up)), (orb_list_dn(1:occ_num_dn))
      print *, ''
    end do
  end do

  deallocate(buffer, orb_list_dn, orb_list_up)

end

4. Reordering basis functions

Here is a simple recipe to reorder your atomic orbitals in TREXIO ordering, and to apply this re-ordering to your MOs.

4.1. Cartesian AOs

  1. Specify the ordering used by your software in a list, for each angular momentum. Here is a simple example for Molden ordering. 

    conv = [ [ 's' ],
         ['x', 'y', 'z'],
         ['xx', 'yy', 'zz', 'xy', 'xz', 'yz'],
         ['xxx', 'yyy', 'zzz', 'xyy', 'xxy', 'xxz', 'xzz', 'yzz', 'yyz', 'xyz'],
         ['xxxx', 'yyyy', 'zzzz', 'xxxy', 'xxxz', 'xyyy', 'yyyz', 'xzzz', 'yzzz', 'xxyy', 'xxzz', 'yyzz', 'xxyz', 'xyyz', 'xyzz'] ]

    conv[l] returns the ordering for angular momentum l.

  2. Iterate over the shells, and append the corresponding predefined lists 

    ao = []
for l in shell_ang_mom:
       ao.append(conv[l])

  3. For each shell, assign a numerical index to each AO, then reorder them lexicographically 

    ao_ordering = []
j = 0
for k,l in enumerate(ao):
      accu = [ x, i+j, for i,x in enumerate(l) ]
      accu.sort()
      ao_ordering += accu
      j += len(l)
ao_ordering = [ i for (_,i) in ao_ordering ]

4.2. Spherical AOs

For spherical AOs, you can adopt a similar approach as for Cartesian AOs, but you will need to define a function that generates a string to ensure lexicographic sorting functions correctly.

  1. Specify the ordering used by your software in a list, for each angular momentum. Here is a simple example for Molden ordering. 

    conv = [ ['s'],
         ['p+1', 'p-1', 'p+0'],
         ['d+0', 'd+1', 'd-1', 'd+2', 'd-2'],
         ['f+0', 'f+1', 'f-1', 'f+2', 'f-2', 'f+3', 'f-3'],
         ['g+0', 'g+1', 'g-1', 'g+2', 'g-2', 'g+3', 'g-3', 'g+4', 'g-4'] ]

  2. Iterate over the shells, and append the corresponding predefined lists 

    ao = []
for l in shell_ang_mom:
       ao.append(conv[l])

  3. Create a sorting function 

    def f_sort(x):
    m = int(x[1:])
    if m>=0:
        return 2*m
    else:
        return -2*m+1

  1. For each shell, assign a numerical index to each AO, then reorder them lexicographically 

    ao_ordering = []
j = 0
for k,l in enumerate(ao):
      accu = [ f_sort(x), i+j, for i,x in enumerate(l) ]
      accu.sort()
      ao_ordering += accu
      j += len(l)

ao_ordering = [ i for (_,i) in ao_ordering ]

4.3. Reordering the MOs

Apply the ordering to each MO, and concatenate them in a 1D array to pass it to trexio: 

MoMatrix = []
for mo in mo_coef:
      vector = np.zeros(ao_num)
      for i, x in mo:
         vector[i] = x
      for i in ao_ordering:
         MoMatrix.append(vector[i])

trexio.write_mo_coefficient(trexio_file, MoMatrix)

Author: TREX-CoE

Created: 2025-08-02 Sat 12:30

Validate