9.2 KiB
Examples
- Accessing sparse quantities (integrals)
- Reading determinants
Accessing sparse quantities (integrals)
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(:,:,:,:)
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
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
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
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
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
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
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.
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
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
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
Terminate
deallocate( D, h0, G, W )
end program
Reading determinants
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
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)
print *, 'ndet', ndet
rc = trexio_get_int64_num(trex_determinant_file, int64_num)
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
end