GauXC in C

GauXC in C#

This guide shows how to use the GauXC library in C applications. We will cover

  • setting up a CMake project to include GauXC as dependency

  • initializing the GauXC runtime environment

  • reading molecule, basis set, and density matrix from an HDF5 input file

  • setting up the integration grid, load balancer, and exchange-correlation integrator

  • performing the exchange-correlation evaluation and outputting the result

Setting up CMake#

GauXC can be most conveniently used via CMake, therefore we will setup a minimal CMake project for a command line driver to use GauXC. Next to GauXC we will use other dependencies as needed.

The directory structure for the project will be

├── CMakeLists.txt
├── app
│   └── main.c
└── cmake
    ├── skala-argtable3.cmake
    ├── skala-dep-versions.cmake
    └── skala-gauxc.cmake

First we create the main CMakeLists.txt to define our project, include our dependencies, and declare our executable.

CMakeLists.txt#
cmake_minimum_required(VERSION 3.20 FATAL_ERROR)

project("Skala" LANGUAGES C CXX)

include(GNUInstallDirs)
include(FetchContent)
set(FETCHCONTENT_UPDATES_DISCONNECTED ON CACHE BOOL "Disable FC Updates")

list(PREPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/cmake)

option(Skala_GauXC_ENABLE_OPENMP "Enable OpenMP support in GauXC" ON)
option(Skala_GauXC_ENABLE_MPI "Enable MPI support in GauXC" OFF)
option(Skala_GauXC_ENABLE_CUDA "Enable CUDA support in GauXC" OFF)

include(skala-dep-versions)
include(skala-gauxc)
include(skala-argtable3)

add_executable(
  "${PROJECT_NAME}"
  "app/main.c"
)
target_link_libraries(
  "${PROJECT_NAME}"
  PRIVATE
  gauxc::gauxc
  argtable3::argtable3
)
install(
  TARGETS
  "${PROJECT_NAME}"
  DESTINATION
  "${CMAKE_INSTALL_BINDIR}"
)

list(REMOVE_AT CMAKE_MODULE_PATH 0)

For handling the dependencies, we create a separate file in the cmake/ subdirectory to include the path and checksums for all our dependencies.

cmake/skala-dep-versions.cmake#
set(Skala_GauXC_URL "https://github.com/microsoft/skala/releases/download/v1.1.1/gauxc-skala.tar.gz")
set(Skala_GauXC_SHA256 "f0535e79006b6e75d0b95aa26950d45738d4fd6f75fed975b75cf20eb993a15f")

set(Skala_Argtable3_URL "https://github.com/argtable/argtable3/releases/download/v3.3.1/argtable-v3.3.1.tar.gz")
set(Skala_Argtable3_SHA256 "c08bca4b88ddb9234726b75455b3b1670d7c864d8daf198eaa7a3b4d41addf2c")

For each dependency we will create a separate CMake include file for finding and making it available. First, we define how we will include GauXC our main dependency. For this we can rely in most cases to discover the GauXC config file, however we also provide a fallback to download and build GauXC in case it is not available in the environment. The options we defined in the main CMake file will be passed through to GauXC to ensure the library provides the requested features. Furthermore, after having GauXC available, we double check whether our requirements for GauXC are satisfied, this is especially necessary for the Skala implementation, which requires the GAUXC_HAS_ONEDFT feature flag.

cmake/skala-gauxc.cmake#
if(NOT DEFINED Skala_GauXC_URL)
  include(skala-dep-versions)
endif()
find_package(gauxc QUIET CONFIG)
if(NOT gauxc_FOUND)
  include(FetchContent)

  message(STATUS "Could not find GauXC... Building GauXC from source")
  message(STATUS "GAUXC URL: ${Skala_GauXC_URL}")

  set(GAUXC_ENABLE_ONEDFT ON CACHE BOOL "" FORCE)
  set(GAUXC_ENABLE_C ON CACHE BOOL "" FORCE)
  set(GAUXC_ENABLE_TESTS OFF CACHE BOOL "" FORCE)
  set(GAUXC_ENABLE_OPENMP ${Skala_GauXC_ENABLE_OPENMP} CACHE BOOL "" FORCE)
  set(GAUXC_ENABLE_MPI ${Skala_GauXC_ENABLE_MPI} CACHE BOOL "" FORCE)
  set(GAUXC_ENABLE_CUDA ${Skala_GauXC_ENABLE_CUDA} CACHE BOOL "" FORCE)

  FetchContent_Declare(
    gauxc
    URL ${Skala_GauXC_URL}
    # URL_HASH SHA256=${Skala_GauXC_SHA256}
    DOWNLOAD_EXTRACT_TIMESTAMP ON
  )
  FetchContent_MakeAvailable(gauxc)

else()
  if(NOT ${GAUXC_HAS_ONEDFT})
    message(FATAL_ERROR "GauXC found but without Skala support enabled")
  endif()
  if(NOT ${GAUXC_HAS_C})
    message(FATAL_ERROR "GauXC found but without C API support enabled")
  endif()
  if(${Skala_GauXC_ENABLE_OPENMP} AND NOT ${GAUXC_HAS_OPENMP})
    message(FATAL_ERROR "GauXC found without OpenMP support but Skala_GauXC_ENABLE_OPENMP is ON")
  endif()
  if(${Skala_GauXC_ENABLE_MPI} AND NOT ${GAUXC_HAS_MPI})
    message(FATAL_ERROR "GauXC found without MPI support but Skala_GauXC_ENABLE_MPI is ON")
  endif()
  if(${Skala_GauXC_ENABLE_CUDA} AND NOT ${GAUXC_HAS_CUDA})
    message(FATAL_ERROR "GauXC found without CUDA support but Skala_GauXC_ENABLE_CUDA is ON")
  endif()
endif()

if(GAUXC_HAS_GAU2GRID AND NOT TARGET gau2grid::gg)
  find_package(gau2grid CONFIG REQUIRED)
endif()

For our command line driver, we will be using Argtable3 to create the command line interface.

cmake/skala-argtable3.cmake#
if(NOT DEFINED Skala_Argtable3_URL)
  include(skala-dep-versions)
endif()
find_package(Argtable3 QUIET CONFIG)
if(NOT Argtable3_FOUND)
  include(FetchContent)

  message(STATUS "Could not find Argtable3... Building Argtable3 from source")
  message(STATUS "Argtable3 URL: ${Skala_Argtable3_URL}")

  FetchContent_Declare(
    argtable3
    URL ${Skala_Argtable3_URL}
    URL_HASH SHA256=${Skala_Argtable3_SHA256}
    DOWNLOAD_EXTRACT_TIMESTAMP ON
  )
  FetchContent_MakeAvailable(Argtable3)
endif()

With this we have the full CMake setup we need for creating our command line driver.

Setting up headers#

For our main driver program we include the relevant headers from GauXC, next to the ones needed for the HDF5 I/O and command line interface.

app/main.c (header includes)#
// Standard C libraries
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

// For GauXC core functionality
#include <gauxc/status.h>
#include <gauxc/molecule.h>
#include <gauxc/basisset.h>
#include <gauxc/molgrid.h>
#include <gauxc/runtime_environment.h>
#include <gauxc/load_balancer.h>
#include <gauxc/molecular_weights.h>
#include <gauxc/functional.h>
#include <gauxc/xc_integrator.h>
#include <gauxc/matrix.h>

// For HDF5 I/O
#include <gauxc/external/hdf5.h>

// For command line interface
#include <argtable3.h>

For each of the GauXC components we will be using, we include the respective header file from GauXC.

gauxc/status.h

For handling GauXC status codes and errors.

gauxc/molecule.h

For handling molecular data structures, like atomic numbers and coordinates.

gauxc/basisset.h

For handling basis set data, like basis function definitions.

gauxc/molgrid.h

For setting up and managing the integration grid.

gauxc/runtime_environment.h

For interacting with MPI and device runtime environments.

gauxc/load_balancer.h

For setting up and managing the load balancer for distributing grid points.

gauxc/molecular_weights.h

For computing molecular weights needed for grid generation.

gauxc/functional.h

For handling exchange-correlation functionals.

gauxc/xc_integrator.h

For setting up and managing the exchange-correlation integrator.

gauxc/matrix.h

For handling matrix data structures, like density matrices and potentials.

We start our main driver program with a guarded call to the MPI initialization.

app/main.c (MPI initialize)#
int
main(int argc, char** argv)
{
#ifdef GAUXC_HAS_MPI
  MPI_Init(NULL, NULL);
#endif

For the finalization of the MPI environment we also add a guarded call to the MPI finalize function at the end of our main program.

app/main.c (MPI finalize)#
#ifdef GAUXC_HAS_MPI
  MPI_Finalize();
#endif
}

GauXC provides a macro GAUXC_HAS_MPI to inform users whether GauXC was built with MPI support. We can use this macro to conditionally compile the MPI initialization and finalization code only when GauXC was built with MPI support.

Creating the Command Line Driver#

We will provide our input data mainly via an HDF5 input file. This file will contain the molecular structure, basis set, and density matrix. Additionally, we have parameters for defining the integration grid and also the parallelization strategy. Our main parameters for the input will therefore be

input_file

an HDF5 input file to provide the molecule, basis and density matrix

model

the model checkpoint we want to evaluate

grid_spec

the grid size specification which defines the number of angular and radial integration points

rad_quad_spec

the radial quadrature scheme which defines the spacing of radial integration points

prune_spec

the pruning scheme for combining the atomic grids to a molecular one

For the command line driver we will use Argtable3 to handle the command line arguments concisely. We define the command line arguments for the input HDF5 file, model type, and other parameters using Argtable3.

app/main.c (command line arguments)#
  struct arg_lit *help; 
  struct arg_file *input_file_;
  struct arg_int *batch_size_;
  struct arg_dbl *basis_tol_;
  struct arg_str *model_, *grid_spec_, *rad_quad_spec_, *prune_spec_;
  struct arg_str *lb_exec_space_, *int_exec_space_;
  struct arg_end *end;

  void* argtable[] = {
    input_file_    = arg_filen(NULL, NULL, "<file>", 1, 1, "Input file containing molecular geometry and density matrix"),
    model_         = arg_strn(NULL, "model", "<str>", 1, 1, "OneDFT model to use, can be a path to a checkpoint"),
    grid_spec_     = arg_strn(NULL, "grid-spec", "<str>", 0, 1, "Atomic grid size specification (default: Fine)"),
                     arg_rem(NULL, "Possible values are: Fine, UltraFine, SuperFine, GM3, GM5"),
    rad_quad_spec_ = arg_strn(NULL, "radial-quad", "<str>", 0, 1, "Radial quadrature scheme (default: MuraKnowles)"),
                     arg_rem(NULL, "Possible values are: Becke, MuraKnowles, TreutlerAhlrichs, MurrayHandyLaming"),
    prune_spec_    = arg_strn(NULL, "prune-scheme", "<str>", 0, 1, "Pruning scheme (default: Robust)"),
                     arg_rem(NULL, "Possible values are: Unpruned, Robust, Treutler"),
    lb_exec_space_ = arg_strn(NULL, "lb-exec-space", "<str>", 0, 1, "Load balancer execution space"),
                     arg_rem(NULL, "Possible values are: Host, Device"),
    int_exec_space_= arg_strn(NULL, "int-exec-space", "<str>", 0, 1, "Integrator execution space"),
                     arg_rem(NULL, "Possible values are: Host, Device"),
    batch_size_    = arg_intn(NULL, "batch-size", "<int>", 0, 1, "Batch size for grid point processing (default: 512)"),
    basis_tol_     = arg_dbln(NULL, "basis-tol", "<double>", 0, 1, "Basis function evaluation tolerance (default: 1e-10)"),
    help           = arg_litn(NULL, "help", 0, 1, "Print this help and exit"),
    end = arg_end(20)
  };

Note

The settings for the molecular grid (grid_spec, rad_quad_spec, prune_spec) are defined in more detail in Molecular grid reference.

With this we can parse the command line and handle potential errors.

app/main.c (parse command line arguments)#
  int nerrors = arg_parse(argc, argv, argtable);
  if(help->count > 0) {
    printf("Usage: %s", argv[0]);
    arg_print_syntax(stdout, argtable, "\n\n");
    printf("Options:\n");
    arg_print_glossary(stdout, argtable, "  %-25s %s\n");
    return EXIT_SUCCESS;
  }

  if (nerrors > 0 || input_file_->count == 0 || model_->count == 0) {
    printf("Usage: %s", argv[0]);
    arg_print_syntax(stdout, argtable, "\n");
    arg_print_errors(stderr, end, argv[0]);
    return EXIT_FAILURE;
  }

Finally, we can extract the values of the command line arguments and store them in variables for later use. For this purpose we will define two helper functions, one for copying the values from the Argtable3 structs and a normalization function to ensure all inputs are lowercase.

app/main.c (helper functions for command line arguments)#
inline static char * arg_strcpy(struct arg_str * arg, char * default_str)
{
  if(arg->count == 0 && default_str != NULL) {
    return strcpy(malloc(strlen(default_str) + 1), default_str);
  }
  if(arg->count > 0) {
    const char* sval = arg->sval[arg->count - 1];
    return strcpy(malloc(strlen(sval) + 1), sval);
  }
  return NULL;
}

inline static char * lowercase(char * str)
{
  for(char * p = str; *p; ++p) *p = tolower(*p);
  return str;
}

With this we can extract the command line argument values.

app/main.c (extract command line argument values)#
  char* input_file = strcpy(malloc(strlen(input_file_->filename[0]) + 1), input_file_->filename[0]);
  char* model = arg_strcpy(model_, NULL);
  char* grid_spec = lowercase(arg_strcpy(grid_spec_, "Fine"));
  char* rad_quad_spec = lowercase(arg_strcpy(rad_quad_spec_, "MuraKnowles"));
  char* prune_spec = lowercase(arg_strcpy(prune_spec_, "Robust"));
  char* lb_exec_space_str = lowercase(arg_strcpy(lb_exec_space_, "Host"));
  char* int_exec_space_str = lowercase(arg_strcpy(int_exec_space_, "Host"));
  int batch_size = batch_size_->count > 0 ? batch_size_->ival[0] : 512;
  double basis_tol = basis_tol_->count > 0 ? basis_tol_->dval[0] : 1e-10;

At this point we can already free the Argtable3 structures as we do not need them anymore.

app/main.c (free Argtable3 structures)#
  // Clean up argtable memory
  arg_freetable(argtable, sizeof(argtable) / sizeof(argtable[0]));

Also, we want to ensure our input string variables will get freed at the end of our program. We add the respective free() calls at the end of our main program, before the MPI finalization.

app/main.c (free input strings)#
  // Clean up memory
  free(input_file);
  free(model);
  free(grid_spec);
  free(rad_quad_spec);
  free(prune_spec);
  free(lb_exec_space_str);
  free(int_exec_space_str);

Before adding any further implementation, we will create a first build of the project. To configure the CMake build run

cmake -B build -G Ninja -S .
cmake --build build

After we build our project successfully, we can run our driver directly from the build directory with

./build/Skala --help

As output you should see the help page generated by our command-line interface

Usage: ./_build/Skala <file> --model=<str> [--grid-spec=<str>] [--radial-quad=<str>] [--prune-scheme=<str>] [--lb-exec-space=<str>] [--int-exec-space=<str>] [--batch-size=<int>] [--basis-tol=<double>] [--help]

Options:
  <file>                    Input file containing molecular geometry and density matrix
  --model=<str>             OneDFT model to use, can be a path to a checkpoint
  --grid-spec=<str>         Atomic grid size specification (default: Fine)
                            Possible values are: Fine, UltraFine, SuperFine, GM3, GM5
  --radial-quad=<str>       Radial quadrature scheme (default: MuraKnowles)
                            Possible values are: Becke, MuraKnowles, TreutlerAhlrichs, MurrayHandyLaming
  --prune-scheme=<str>      Pruning scheme (default: Robust)
                            Possible values are: Unpruned, Robust, Treutler
  --lb-exec-space=<str>     Load balancer execution space
                            Possible values are: Host, Device
  --int-exec-space=<str>    Integrator execution space
                            Possible values are: Host, Device
  --batch-size=<int>        Batch size for grid point processing (default: 512)
  --basis-tol=<double>      Basis function evaluation tolerance (default: 1e-10)
  --help                    Print this help and exit

With this we are able to change our configuration conveniently from the command line.

Tip

If you encounter any issues when running the driver, like segmentation faults, rerun your binary with a debugger, like gdb.

gdb --args ./build/Skala --help

This way you can inspect the stack trace and find the source of the error more easily.

As first step for any interaction with GauXC, we need to initialize the GauXC status object and runtime environment.

app/main.c (GauXC status and runtime environment)#
  // Add handler for events like exceptions
  GauXCStatus status;

  // Create runtime
  GauXCRuntimeEnvironment rt = gauxc_runtime_environment_new(&status GAUXC_MPI_CODE(, MPI_COMM_WORLD));
  int world_rank = gauxc_runtime_environment_comm_rank(&status, rt);
  int world_size = gauxc_runtime_environment_comm_size(&status, rt);

The world_size and world_rank variables describe our MPI environment and contain dummy values if we do not use MPI.

At this point we want to show the configuration we are using before proceeding further. We can use the world_rank variable to ensure that only the root process outputs the configuration in case of MPI parallel execution.

app/main.c (configuration summary)#
  if(!world_rank && !status.code) {
    printf("Configuration\n");
    printf("-> Input file        : %s\n", input_file);
    printf("-> Model             : %s\n", model);
    printf("-> Grid              : %s\n", grid_spec);
    printf("-> Radial quadrature : %s\n", rad_quad_spec);
    printf("-> Pruning scheme    : %s\n", prune_spec);
    printf("\n");
  }

Finally, we also add calls to free the runtime environment at the end of our main program. To ease freeing memory we will define a type generic macro for calling the respective free functions.

app/main.c (free macro)#
// Generic delete macro
#define gauxc_delete(status, ptr) _Generic( (ptr), \
                 GauXCMolecule : gauxc_molecule_delete, \
                 GauXCBasisSet : gauxc_basisset_delete, \
                  GauXCMolGrid : gauxc_molgrid_delete, \
       GauXCRuntimeEnvironment : gauxc_runtime_environment_delete, \
             GauXCLoadBalancer : gauxc_load_balancer_delete, \
      GauXCLoadBalancerFactory : gauxc_load_balancer_factory_delete, \
         GauXCMolecularWeights : gauxc_molecular_weights_delete, \
  GauXCMolecularWeightsFactory : gauxc_molecular_weights_factory_delete, \
               GauXCFunctional : gauxc_functional_delete, \
               GauXCIntegrator : gauxc_integrator_delete, \
        GauXCIntegratorFactory : gauxc_integrator_factory_delete, \
                   GauXCMatrix : gauxc_matrix_delete \
                               )(status, &(ptr))

With this we can free the runtime environment at the end of our main program.

app/main.c (free runtime)#
  gauxc_delete(&status, rt);

Molecule data#

For reading the molecule we will use GauXC’s built-in functionality to read from an HDF5 dataset.

Note

GauXC stores the information about the atomic numbers and their coordinates as an array of structs, defined in C++ as

struct GauXC::Atom {
  int64_t Z; ///< Atomic number
  double x;  ///< X coordinate (bohr)
  double y;  ///< Y coordinate (bohr)
  double z;  ///< Z coordinate (bohr)
};
class GauXC::Molecule : public std::vector<GauXC::Atom> {
  ...
};

The HDF5 wrapper directly maps this struct representation to an HDF5 dataset.

We use gauxc_molecule_read_hdf5_record function which implements the reading of the molecule data.

app/main.c (read molecule)#
  // Get molecule (atomic numbers and cartesian coordinates)
  GauXCMolecule mol = gauxc_molecule_new(&status);
  // Load molecule from HDF5 dataset
  gauxc_molecule_read_hdf5_record(&status, mol, input_file, "/MOLECULE");

To ensure proper memory management we also add a call to free the molecule object at the end of our main program.

app/main.c (free molecule)#
  gauxc_delete(&status, mol);

Basis set data#

For the basis set we will use the same approach as for the molecule and use GauXC’s built-in HDF5 reading functionality.

Note

Similar to the molecule the basis set object is built as an array of shell objects. The shell object contains the information about the primitives, contraction coefficients, angular momentum, and center of the shell.

template <typename F>
class alignas(256) GauXC::Shell {
  std::array<F, 32> alpha; ///< exponents of primitives
  std::array<F, 32> coeff; ///< contraction coefficients
  std::array<double, 3> O; ///< origin of the shell
  int32_t nprim; ///< number of primitives
  int32_t l; ///< angular moment of the shell
  int32_t pure; ///< pure=1: spherical Gaussianss; pure=0: cartesian Gaussianss
};

template <typename F>
struct GauXC::BasisSet : public std::vector<GauXC::Shell<F>> {
  ...
};

Again, this allows to directly map the object’s representation to an HDF5 dataset.

With GauXC’s gauxc_basisset_read_hdf5_record function we can read the basis set data conveniently from the HDF5 file.

app/main.c (read basisset)#
  // Get basis set
  GauXCBasisSet basis = gauxc_basisset_new(&status);
  // Load basis set from HDF5 dataset
  gauxc_basisset_read_hdf5_record(&status, basis, input_file, "/BASIS");

To avoid memory leaks we also add a call to free the basis set object at the end of our main program.

app/main.c (free basisset)#
  gauxc_delete(&status, basis);

Integration grid#

To setup the integration grid, which is the part of the input to GauXC for computing the exchange-correlation functional, we create a molecular grid. We have three main input parameters here, the grid size which defines the density of angular and radial points, the radial quadrature scheme which defines the spacing of the radial points, and the pruning scheme which defines how atomic grids are combined to a molecular grid. In GauXC these are defined as enumerators and we add a number of helper functions for turning the input strings from the command-line to the respective enumerator values.

app/main.c (enumerator conversion functions)#
enum GauXC_RadialQuad
read_radial_quad(GauXCStatus* status, const char* rad_quad_spec)
{
  status->code = 0;
  if(strcmp(rad_quad_spec, "becke") == 0)
    return GauXC_RadialQuad_Becke;
  if(strcmp(rad_quad_spec, "muraknowles") == 0)
    return GauXC_RadialQuad_MuraKnowles;
  if(strcmp(rad_quad_spec, "treutlerahlrichs") == 0)
    return GauXC_RadialQuad_TreutlerAhlrichs;
  if(strcmp(rad_quad_spec, "murrayhandylaming") == 0)
    return GauXC_RadialQuad_MurrayHandyLaming;
  status->code = 1;
}

enum GauXC_AtomicGridSizeDefault
read_atomic_grid_size(GauXCStatus* status, const char* spec)
{
    status->code = 0;
    if(strcmp(spec, "fine") == 0)
      return GauXC_AtomicGridSizeDefault_FineGrid;
    if(strcmp(spec, "ultrafine") == 0)
      return GauXC_AtomicGridSizeDefault_UltraFineGrid;
    if(strcmp(spec, "superfine") == 0)
      return GauXC_AtomicGridSizeDefault_SuperFineGrid;
    if(strcmp(spec, "gm3") == 0)
      return GauXC_AtomicGridSizeDefault_GM3;
    if(strcmp(spec, "gm5") == 0)
      return GauXC_AtomicGridSizeDefault_GM5;
    status->code = 1;
}

enum GauXC_PruningScheme
read_pruning_scheme(GauXCStatus* status, const char* spec)
{
  status->code = 0;
  if(strcmp(spec, "Unpruned") == 0)
    return GauXC_PruningScheme_Unpruned;
  if(strcmp(spec, "Robust") == 0)
    return GauXC_PruningScheme_Robust;
  if(strcmp(spec, "Treutler") == 0)
    return GauXC_PruningScheme_Treutler;
  status->code = 1;
}

For the main program we can now create the molecular grid based on our input parameters. We also have to define the batch size for the grid, the default is 512 points per batch, however larger values up around 10000 are recommended for better performance.

app/main.c (grid setup)#
  // Define molecular grid from grid size, radial quadrature and pruning scheme
  enum GauXC_AtomicGridSizeDefault grid_type = read_atomic_grid_size(&status, grid_spec);
  enum GauXC_RadialQuad radial_quad = read_radial_quad(&status, rad_quad_spec);
  enum GauXC_PruningScheme pruning_scheme = read_pruning_scheme(&status, prune_spec);
  GauXCMolGrid grid = gauxc_molgrid_new_default(
    &status,
    mol,
    pruning_scheme,
    batch_size,
    radial_quad,
    grid_type);

We free the molecular grid object at the end of our main program.

app/main.c (free grid)#
  gauxc_delete(&status, grid);

Exchange-correlation integrator#

To distribute the work of evaluating the exchange-correlation functional on the grid points, we create a load balancer. The load balancer will take care of distributing the grid points to the available resources, either host or device, based on the execuation space we provide. Again we have a helper function to convert the input string to the respective enumerator value.

app/main.c (execuation space enumerator)#
enum GauXC_ExecutionSpace
read_execution_space(GauXCStatus* status, const char* exec_space_str)
{
  status->code = 0;
  if(strcmp(exec_space_str, "host") == 0)
    return GauXC_ExecutionSpace_Host;
  if(strcmp(exec_space_str, "device") == 0)
    return GauXC_ExecutionSpace_Device;
  status->code = 1;
}

Note that the load balancer will provide access to the molecule, basis and grid data for all further usage in GauXC. We can now create the load balancer based on our input parameters.

app/main.c (load balancer setup)#
  // Choose whether we run on host or device
  enum GauXC_ExecutionSpace lb_exec_space, int_exec_space;
  lb_exec_space = read_execution_space(&status, lb_exec_space_str);
  int_exec_space = read_execution_space(&status, int_exec_space_str);

  // Setup load balancer based on molecule, grid and basis set
  GauXCLoadBalancerFactory lb_factory = gauxc_load_balancer_factory_new(&status, lb_exec_space, "Replicated");
  GauXCLoadBalancer lb = gauxc_load_balancer_factory_get_shared_instance(&status, lb_factory, rt, mol, grid, basis);

  // Apply partitioning weights to the molecule grid
  GauXCMolecularWeightsSettings settings = {GauXC_XCWeightAlg_SSF, false};
  GauXCMolecularWeightsFactory mw_factory = gauxc_molecular_weights_factory_new(&status, int_exec_space,
    "Default", settings);
  GauXCMolecularWeights mw = gauxc_molecular_weights_factory_get_instance(&status, mw_factory);
  gauxc_molecular_weights_modify_weights(&status, mw, lb);

Finally, we can create the main GauXC integrator, for this we setup the exchange-correlation integrator factory for producing an instance of the integrator. To configure the integrator we create an additional settings object which holds the model checkpoint we want to evaluate.

app/main.c (integrator setup)#
  // Setup exchange-correlation integrator
  GauXCFunctional func = gauxc_functional_from_string(&status, "PBE", true);
  GauXCIntegratorFactory integrator_factory = gauxc_integrator_factory_new(&status, int_exec_space,
    "Replicated", "Default", "Default", "Default");
  GauXCIntegrator integrator = gauxc_integrator_factory_get_instance(&status, integrator_factory, func, lb);

We free the integrator and its associated objects at the end of our main program.

app/main.c (free integrator)#
  gauxc_delete(&status, lb_factory);
  gauxc_delete(&status, lb);
  gauxc_delete(&status, mw_factory);
  gauxc_delete(&status, mw);
  gauxc_delete(&status, func);
  gauxc_delete(&status, integrator_factory);
  gauxc_delete(&status, integrator);

Density matrix#

The final input we need to provide to GauXC is the density matrix. Similar to the molecule and basis set we will read it from our HDF5 input file.

app/main.c (read density matrix)#
  // Load density matrix from input
  GauXCMatrix P_s = gauxc_matrix_empty(&status);
  GauXCMatrix P_z = gauxc_matrix_empty(&status);
  gauxc_matrix_read_hdf5_record(&status, P_s, input_file, "/DENSITY_SCALAR");
  gauxc_matrix_read_hdf5_record(&status, P_z, input_file, "/DENSITY_Z");

Also for the density matrix we add a call to free the matrix object at the end of our main program.

app/main.c (free density matrix)#
  gauxc_delete(&status, P_s);
  gauxc_delete(&status, P_z);

Exchange-correlation evaluation#

With all inputs provided we can now perform the exchange-correlation evaluation.

app/main.c (exchange-correlation evaluation)#
  // Integrate exchange correlation energy
  double EXC;
  GauXCMatrix VXC_s = gauxc_matrix_empty(&status);
  GauXCMatrix VXC_z = gauxc_matrix_empty(&status);

  gauxc_integrator_eval_exc_vxc_onedft_uks(&status, integrator, P_s, P_z, model, &EXC, &VXC_s, &VXC_z);

After the evaluation we can output the computed exchange-correlation energy.

app/main.c (exchange-correlation output)#
  if(!world_rank && !status.code) {
    printf("Results\n");
    printf("-> EXC : %.10f\n", EXC);
    printf("\n");
  }

Since we have allocated new matrices for the exchange-correlation potential, we also need to free them at the end of our main program.

app/main.c (free exchange-correlation potential)#
  gauxc_delete(&status, VXC_s);
  gauxc_delete(&status, VXC_z);

Now we can rebuild our project with

cmake --build build

After we build our project successfully, we run the driver again from the build directory

./build/Skala He_def2-svp.h5 --model PBE

Note

The He_def2-svp.h5 input file can be created with the skala package.

from pyscf import gto

from skala.gauxc.export import write_gauxc_h5_from_pyscf
from skala.pyscf import SkalaRKS

mol = gto.M(atom="He 0 0 0", basis="def2-svp", unit="Bohr", spin=0)
ks = SkalaRKS(mol, xc="pbe")
ks.kernel()

dm = ks.make_rdm1()
exc = ks.scf_summary["exc"]
_, _, vxc = ks._numint.nr_rks(ks.mol, ks.grids, ks.xc, dm)

write_gauxc_h5_from_pyscf("He_def2-svp.h5", mol, dm=dm, exc=exc, vxc=vxc)

As output we can see the results for the PBE functional

Configuration
-> Input file        : He_def2-svp.h5
-> Model             : PBE
-> Grid              : fine
-> Radial quadrature : muraknowles
-> Pruning scheme    : robust

Results
-> EXC : -1.0540318683

Download checkpoint from HuggingFace#

To evaluate Skala we first need to download the model checkpoint from HuggingFace. For this we can use the hf command line tool from the huggingface_hub Python package. After downloading the model checkpoint we can run our driver again with the new model

hf download microsoft/skala skala-1.0.fun --local-dir .
./build/Skala He_def2-svp.h5 --model ./skala-1.0.fun

In the output we can see the results for the Skala functional

Configuration
-> Input file        : He_def2-svp.h5
-> Model             : ./skala-1.0.fun
-> Grid              : fine
-> Radial quadrature : muraknowles
-> Pruning scheme    : robust

Results
-> EXC : -1.0712560886

Full source code#

Full source code of the main driver
app/main.c#
// Standard C libraries
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>

// For GauXC core functionality
#include <gauxc/status.h>
#include <gauxc/molecule.h>
#include <gauxc/basisset.h>
#include <gauxc/molgrid.h>
#include <gauxc/runtime_environment.h>
#include <gauxc/load_balancer.h>
#include <gauxc/molecular_weights.h>
#include <gauxc/functional.h>
#include <gauxc/xc_integrator.h>
#include <gauxc/matrix.h>

// For HDF5 I/O
#include <gauxc/external/hdf5.h>

// For command line interface
#include <argtable3.h>

// Generic delete macro
#define gauxc_delete(status, ptr) _Generic( (ptr), \
                 GauXCMolecule : gauxc_molecule_delete, \
                 GauXCBasisSet : gauxc_basisset_delete, \
                  GauXCMolGrid : gauxc_molgrid_delete, \
       GauXCRuntimeEnvironment : gauxc_runtime_environment_delete, \
             GauXCLoadBalancer : gauxc_load_balancer_delete, \
      GauXCLoadBalancerFactory : gauxc_load_balancer_factory_delete, \
         GauXCMolecularWeights : gauxc_molecular_weights_delete, \
  GauXCMolecularWeightsFactory : gauxc_molecular_weights_factory_delete, \
               GauXCFunctional : gauxc_functional_delete, \
               GauXCIntegrator : gauxc_integrator_delete, \
        GauXCIntegratorFactory : gauxc_integrator_factory_delete, \
                   GauXCMatrix : gauxc_matrix_delete \
                               )(status, &(ptr))

enum GauXC_ExecutionSpace
read_execution_space(GauXCStatus* status, const char* exec_space_str)
{
  status->code = 0;
  if(strcmp(exec_space_str, "host") == 0)
    return GauXC_ExecutionSpace_Host;
  if(strcmp(exec_space_str, "device") == 0)
    return GauXC_ExecutionSpace_Device;
  status->code = 1;
}

enum GauXC_RadialQuad
read_radial_quad(GauXCStatus* status, const char* rad_quad_spec)
{
  status->code = 0;
  if(strcmp(rad_quad_spec, "becke") == 0)
    return GauXC_RadialQuad_Becke;
  if(strcmp(rad_quad_spec, "muraknowles") == 0)
    return GauXC_RadialQuad_MuraKnowles;
  if(strcmp(rad_quad_spec, "treutlerahlrichs") == 0)
    return GauXC_RadialQuad_TreutlerAhlrichs;
  if(strcmp(rad_quad_spec, "murrayhandylaming") == 0)
    return GauXC_RadialQuad_MurrayHandyLaming;
  status->code = 1;
}

enum GauXC_AtomicGridSizeDefault
read_atomic_grid_size(GauXCStatus* status, const char* spec)
{
    status->code = 0;
    if(strcmp(spec, "fine") == 0)
      return GauXC_AtomicGridSizeDefault_FineGrid;
    if(strcmp(spec, "ultrafine") == 0)
      return GauXC_AtomicGridSizeDefault_UltraFineGrid;
    if(strcmp(spec, "superfine") == 0)
      return GauXC_AtomicGridSizeDefault_SuperFineGrid;
    if(strcmp(spec, "gm3") == 0)
      return GauXC_AtomicGridSizeDefault_GM3;
    if(strcmp(spec, "gm5") == 0)
      return GauXC_AtomicGridSizeDefault_GM5;
    status->code = 1;
}

enum GauXC_PruningScheme
read_pruning_scheme(GauXCStatus* status, const char* spec)
{
  status->code = 0;
  if(strcmp(spec, "Unpruned") == 0)
    return GauXC_PruningScheme_Unpruned;
  if(strcmp(spec, "Robust") == 0)
    return GauXC_PruningScheme_Robust;
  if(strcmp(spec, "Treutler") == 0)
    return GauXC_PruningScheme_Treutler;
  status->code = 1;
}

inline static char * arg_strcpy(struct arg_str * arg, char * default_str)
{
  if(arg->count == 0 && default_str != NULL) {
    return strcpy(malloc(strlen(default_str) + 1), default_str);
  }
  if(arg->count > 0) {
    const char* sval = arg->sval[arg->count - 1];
    return strcpy(malloc(strlen(sval) + 1), sval);
  }
  return NULL;
}

inline static char * lowercase(char * str)
{
  for(char * p = str; *p; ++p) *p = tolower(*p);
  return str;
}

int
main(int argc, char** argv)
{
#ifdef GAUXC_HAS_MPI
  MPI_Init(NULL, NULL);
#endif
  struct arg_lit *help; 
  struct arg_file *input_file_;
  struct arg_int *batch_size_;
  struct arg_dbl *basis_tol_;
  struct arg_str *model_, *grid_spec_, *rad_quad_spec_, *prune_spec_;
  struct arg_str *lb_exec_space_, *int_exec_space_;
  struct arg_end *end;

  void* argtable[] = {
    input_file_    = arg_filen(NULL, NULL, "<file>", 1, 1, "Input file containing molecular geometry and density matrix"),
    model_         = arg_strn(NULL, "model", "<str>", 1, 1, "OneDFT model to use, can be a path to a checkpoint"),
    grid_spec_     = arg_strn(NULL, "grid-spec", "<str>", 0, 1, "Atomic grid size specification (default: Fine)"),
                     arg_rem(NULL, "Possible values are: Fine, UltraFine, SuperFine, GM3, GM5"),
    rad_quad_spec_ = arg_strn(NULL, "radial-quad", "<str>", 0, 1, "Radial quadrature scheme (default: MuraKnowles)"),
                     arg_rem(NULL, "Possible values are: Becke, MuraKnowles, TreutlerAhlrichs, MurrayHandyLaming"),
    prune_spec_    = arg_strn(NULL, "prune-scheme", "<str>", 0, 1, "Pruning scheme (default: Robust)"),
                     arg_rem(NULL, "Possible values are: Unpruned, Robust, Treutler"),
    lb_exec_space_ = arg_strn(NULL, "lb-exec-space", "<str>", 0, 1, "Load balancer execution space"),
                     arg_rem(NULL, "Possible values are: Host, Device"),
    int_exec_space_= arg_strn(NULL, "int-exec-space", "<str>", 0, 1, "Integrator execution space"),
                     arg_rem(NULL, "Possible values are: Host, Device"),
    batch_size_    = arg_intn(NULL, "batch-size", "<int>", 0, 1, "Batch size for grid point processing (default: 512)"),
    basis_tol_     = arg_dbln(NULL, "basis-tol", "<double>", 0, 1, "Basis function evaluation tolerance (default: 1e-10)"),
    help           = arg_litn(NULL, "help", 0, 1, "Print this help and exit"),
    end = arg_end(20)
  };

  int nerrors = arg_parse(argc, argv, argtable);
  if(help->count > 0) {
    printf("Usage: %s", argv[0]);
    arg_print_syntax(stdout, argtable, "\n\n");
    printf("Options:\n");
    arg_print_glossary(stdout, argtable, "  %-25s %s\n");
    return EXIT_SUCCESS;
  }

  if (nerrors > 0 || input_file_->count == 0 || model_->count == 0) {
    printf("Usage: %s", argv[0]);
    arg_print_syntax(stdout, argtable, "\n");
    arg_print_errors(stderr, end, argv[0]);
    return EXIT_FAILURE;
  }

  char* input_file = strcpy(malloc(strlen(input_file_->filename[0]) + 1), input_file_->filename[0]);
  char* model = arg_strcpy(model_, NULL);
  char* grid_spec = lowercase(arg_strcpy(grid_spec_, "Fine"));
  char* rad_quad_spec = lowercase(arg_strcpy(rad_quad_spec_, "MuraKnowles"));
  char* prune_spec = lowercase(arg_strcpy(prune_spec_, "Robust"));
  char* lb_exec_space_str = lowercase(arg_strcpy(lb_exec_space_, "Host"));
  char* int_exec_space_str = lowercase(arg_strcpy(int_exec_space_, "Host"));
  int batch_size = batch_size_->count > 0 ? batch_size_->ival[0] : 512;
  double basis_tol = basis_tol_->count > 0 ? basis_tol_->dval[0] : 1e-10;

  // Clean up argtable memory
  arg_freetable(argtable, sizeof(argtable) / sizeof(argtable[0]));

  // Add handler for events like exceptions
  GauXCStatus status;

  // Create runtime
  GauXCRuntimeEnvironment rt = gauxc_runtime_environment_new(&status GAUXC_MPI_CODE(, MPI_COMM_WORLD));
  int world_rank = gauxc_runtime_environment_comm_rank(&status, rt);
  int world_size = gauxc_runtime_environment_comm_size(&status, rt);

  if(!world_rank && !status.code) {
    printf("Configuration\n");
    printf("-> Input file        : %s\n", input_file);
    printf("-> Model             : %s\n", model);
    printf("-> Grid              : %s\n", grid_spec);
    printf("-> Radial quadrature : %s\n", rad_quad_spec);
    printf("-> Pruning scheme    : %s\n", prune_spec);
    printf("\n");
  }

  // Get molecule (atomic numbers and cartesian coordinates)
  GauXCMolecule mol = gauxc_molecule_new(&status);
  // Load molecule from HDF5 dataset
  gauxc_molecule_read_hdf5_record(&status, mol, input_file, "/MOLECULE");

  // Get basis set
  GauXCBasisSet basis = gauxc_basisset_new(&status);
  // Load basis set from HDF5 dataset
  gauxc_basisset_read_hdf5_record(&status, basis, input_file, "/BASIS");

  // Define molecular grid from grid size, radial quadrature and pruning scheme
  enum GauXC_AtomicGridSizeDefault grid_type = read_atomic_grid_size(&status, grid_spec);
  enum GauXC_RadialQuad radial_quad = read_radial_quad(&status, rad_quad_spec);
  enum GauXC_PruningScheme pruning_scheme = read_pruning_scheme(&status, prune_spec);
  GauXCMolGrid grid = gauxc_molgrid_new_default(
    &status,
    mol,
    pruning_scheme,
    batch_size,
    radial_quad,
    grid_type);

  // Choose whether we run on host or device
  enum GauXC_ExecutionSpace lb_exec_space, int_exec_space;
  lb_exec_space = read_execution_space(&status, lb_exec_space_str);
  int_exec_space = read_execution_space(&status, int_exec_space_str);

  // Setup load balancer based on molecule, grid and basis set
  GauXCLoadBalancerFactory lb_factory = gauxc_load_balancer_factory_new(&status, lb_exec_space, "Replicated");
  GauXCLoadBalancer lb = gauxc_load_balancer_factory_get_shared_instance(&status, lb_factory, rt, mol, grid, basis);

  // Apply partitioning weights to the molecule grid
  GauXCMolecularWeightsSettings settings = {GauXC_XCWeightAlg_SSF, false};
  GauXCMolecularWeightsFactory mw_factory = gauxc_molecular_weights_factory_new(&status, int_exec_space,
    "Default", settings);
  GauXCMolecularWeights mw = gauxc_molecular_weights_factory_get_instance(&status, mw_factory);
  gauxc_molecular_weights_modify_weights(&status, mw, lb);

  // Setup exchange-correlation integrator
  GauXCFunctional func = gauxc_functional_from_string(&status, "PBE", true);
  GauXCIntegratorFactory integrator_factory = gauxc_integrator_factory_new(&status, int_exec_space,
    "Replicated", "Default", "Default", "Default");
  GauXCIntegrator integrator = gauxc_integrator_factory_get_instance(&status, integrator_factory, func, lb);

  // Configure model checkpoint
  // GauXCOneDFTSettings onedft_settings;
  // onedft_settings.model = model;

  // Load density matrix from input
  GauXCMatrix P_s = gauxc_matrix_empty(&status);
  GauXCMatrix P_z = gauxc_matrix_empty(&status);
  gauxc_matrix_read_hdf5_record(&status, P_s, input_file, "/DENSITY_SCALAR");
  gauxc_matrix_read_hdf5_record(&status, P_z, input_file, "/DENSITY_Z");

#ifdef GAUXC_HAS_MPI
  MPI_Barrier(MPI_COMM_WORLD);
#endif

  // Integrate exchange correlation energy
  double EXC;
  GauXCMatrix VXC_s = gauxc_matrix_empty(&status);
  GauXCMatrix VXC_z = gauxc_matrix_empty(&status);

  gauxc_integrator_eval_exc_vxc_onedft_uks(&status, integrator, P_s, P_z, model, &EXC, &VXC_s, &VXC_z);

#ifdef GAUXC_HAS_MPI
  MPI_Barrier(MPI_COMM_WORLD);
#endif

  if(!world_rank && !status.code) {
    printf("Results\n");
    printf("-> EXC : %.10f\n", EXC);
    printf("\n");
  }

  // Clean up memory
  free(input_file);
  free(model);
  free(grid_spec);
  free(rad_quad_spec);
  free(prune_spec);
  free(lb_exec_space_str);
  free(int_exec_space_str);
  gauxc_delete(&status, mol);
  gauxc_delete(&status, basis);
  gauxc_delete(&status, grid);
  gauxc_delete(&status, rt);
  gauxc_delete(&status, lb_factory);
  gauxc_delete(&status, lb);
  gauxc_delete(&status, mw_factory);
  gauxc_delete(&status, mw);
  gauxc_delete(&status, func);
  gauxc_delete(&status, integrator_factory);
  gauxc_delete(&status, integrator);
  gauxc_delete(&status, P_s);
  gauxc_delete(&status, P_z);
  gauxc_delete(&status, VXC_s);
  gauxc_delete(&status, VXC_z);

#ifdef GAUXC_HAS_MPI
  MPI_Finalize();
#endif
}

Summary#

Within this guide the usage of GauXC in C applications was demonstrated. A minimal CMake project was created to setup the build environment and include GauXC. We created a command line driver which reads the molecule, basis set, and density matrix from an HDF5 input file and sets up the integration grid, load balancer, and exchange-correlation integrator. Finally, we evaluated the exchange-correlation energy and potential and output the result.

Contents