Solving Kohn-Sham density functional theory: I

The simplest way to use PEXSI to solve Kohn-Sham density functional theory is to use the PPEXSIDFTDriver routine. This routine uses built-in heuristics to obtain values of some parameters in PEXSI and provides a relatively small set of adjustable parameters for users to tune. This routine estimates the chemical potential self-consistently using a combined approach of inertia counting procedure and Newton's iteration through PEXSI. Some heuristic approach is also implemented in this routine for dynamic adjustment of the chemical potential and some stopping criterion.

An example routine is given in driver_ksdft.c, which solves a fake DFT problem by taking a Hamiltonian matrix from examples/lap2dr.matrix.

Here is the structure of the code using the simple driver routine.

#include  "c_pexsi_interface.h"
...
{
  /* Setup the input matrix in distributed compressed sparse column (CSC) format */ 
  ...;
  /* Initialize PEXSI. 
   * PPEXSIPlan is a handle communicating with the C++ internal data structure */
  PPEXSIPlan   plan;

  plan = PPEXSIPlanInitialize( 
      MPI_COMM_WORLD, 
      nprow,
      npcol,
      mpirank, 
      &info );

  /* Tuning parameters of PEXSI. See PPEXSIOption for explanation of the
   * parameters */
  PPEXSIOptions  options;
  PPEXSISetDefaultOptions( &options );

  options.numPole  = 60;
  options.temperature  = 0.019; // 3000K
  options.muPEXSISafeGuard  = 0.2; 
  options.numElectronPEXSITolerance = 0.001;

  /* Load the matrix into the internal data structure */
  PPEXSILoadRealSymmetricHSMatrix( 
      plan, 
      options,
      nrows,
      nnz,
      nnzLocal,
      numColLocal,
      colptrLocal,
      rowindLocal,
      HnzvalLocal,
      isSIdentity,
      SnzvalLocal,
      &info );

  /* Call the simple DFT driver using PEXSI */
  PPEXSIDFTDriver(
      plan,
      options,
      numElectronExact,
      &muPEXSI,                   
      &numElectronPEXSI,         
      &muMinInertia,              
      &muMaxInertia,             
      &numTotalInertiaIter,   
      &numTotalPEXSIIter,   
      &info );

  /* Retrieve the density matrix and other quantities from the plan */

  PPEXSIRetrieveRealSymmetricDFTMatrix(
      plan,
      DMnzvalLocal,
      EDMnzvalLocal,
      FDMnzvalLocal,
      &totalEnergyH,
      &totalEnergyS,
      &totalFreeEnergy,
      &info );

  /* Clean up */
  PPEXSIPlanFinalize(
      plan,
      &info );
}