Adding New Functionality (Example)

Overview

On this Wiki page, we will sketch out how one might go about adding new functionality to NTPoly. The example we give will be based on trying to add a McWeeny purification module. Clone a copy of NTPoly from the repository, and create a new branch to get started.

git clone https://github.com/william-dawson/NTPoly.git

git checkout -b McWeeny

Theory

McWeeny purification can be used to compute the density matrix from a Hamiltonian using a very tight recursive scheme:

P_{k+1} = 3P_k^2 - 2P_k^3

P_0 = \frac{\lambda}{2}(\mu I - H) + \frac{1}{2}I

Here is a potential python implementation:

p_0 = (lambda_val/2.0) * (chemical_potential*identity_mat - hamiltonian) + 0.5*identity_mat
for i in range(0,100):
  p2 = numpy.dot(p_0,p_0)
  p3 = numpy.dot(p_0,p2)
  p0 = 3.0*p2 - 2.0*p3

Getting Started

First, it is important to build an outline of the routine so that it will be easy to test. Let's begin by adding a new subroutine in Source/Fortran/DensityMatrixSolversModule.F90.

Fortran Code Outline

subroutine SolveMcWeeny(Hamiltonian, InverseSquareRoot, Density, &
  & chemical_potential, energy_value_out solver_parameters_in)
  type(Matrix_ps), intent(in) :: Hamiltonian
  type(Matrix_ps), intent(in) :: InverseSquareRoot
  integer, intent(in) :: nel
  type(Matrix_ps), intent(inout) :: Density
  real(ntreal), intent(out), optional :: energy_value_out
  real(ntreal), intent(out) :: chemical_potential
  type(SolverParameters_t), intent(in), optional :: solver_parameters_in
  !! Handling Optional Parameters
  type(SolverParameters_t) :: solver_parameters

  !! Optional Parameters
  if (present(solver_parameters_in)) then
    solver_parameters = solver_parameters_in
  else
    solver_parameters = SolverParameters_t()
  end if

end subroutine SolveMcWeeny

This new module will take in the necessary matrices along with the solver parameters and the chemical potential. Before coding the actual solver, we need to wrap up this module so it an be called from C++. This time, we go to the Source/Wrappers/ directory, and add to the DensityMatrixSolvers_wrp.f90 file.

Wrapper Code Outline

subroutine SolveMcWeeny_wrp(ih_Hamiltonian, ih_InverseSquareRoot, &
  & ih_Density, chemical_potential, energy_value_out, ih_solver_parameters) &
  & bind(c,name="SolveMcWeeny")
  integer(kind=c_int), intent(in) :: ih_Hamiltonian(SIZE_wrp)
  integer(kind=c_int), intent(in) :: ih_InverseSquareRoot(SIZE_wrp)
  integer(kind=c_int), intent(inout) :: ih_Density(SIZE_wrp)
  real(ntreal), intent(out) :: energy_value_out
  real(ntreal), intent(in) :: chemical_potential
  integer(kind=c_int), intent(in) :: ih_solver_parameters(SIZE_wrp)
  type(Matrix_ps_wrp) :: h_Hamiltonian
  type(Matrix_ps_wrp) :: h_InverseSquareRoot
  type(Matrix_ps_wrp) :: h_Density
  type(Matrix_ps_wrp) :: h_solver_parameters

  h_Hamiltonian = transfer(ih_Hamiltonian,h_Hamiltonian)
  h_InverseSquareRoot = transfer(ih_InverseSquareRoot,h_InverseSquareRoot)
  h_Density = transfer(ih_Density,h_Density)
  h_solver_parameters = transfer(ih_solver_parameters, h_solver_parameters)

  call SolveMcWeeny(h_Hamiltonian%data, h_InverseSquareRoot%data, &
       & h_Density%data, chemical_potential, h_solver_parameters%data)
end subroutine SolveMcWeeny_wrp

This wrapper is necessary to hide the distributed datatypes from C. Writing a subroutine like this is mostly just copy and pasting. For all the derived types, you need to pass an integer handle to it, and then use the transfer function to convert handle to the real object. We also need to add a signature to the C header file Source/C/DensityMatrixSolvers_c.h:

Next we add an object oriented C++ layer in Source/CPlusPlus.

C++ Header File

static void McWeeny(const Matrix_ps &Hamiltonian,
               const Matrix_ps &InverseSquareRoot, int nel,
               Matrix_ps &Density, double &chemical_potential,
               double &energy_value_out,
               const SolverParameters &solver_parameters);

C++ Source File

void DensityMatrixSolvers::McWeeny(const Matrix_ps &Hamiltonian,
                            const Matrix_ps &Overlap, int nel,
                            Matrix_ps &Density,
                            double &chemical_potential_out,
                            double &energy_value_out,
                            const SolverParameters &solver_parameters) {
  McWeeny_wrp(GetIH(Hamiltonian), GetIH(Overlap), &nel, GetIH(Density),
         &chemical_potential, &energy_value_out, GetIH(solver_parameters));
}

This is also pretty much just copy and pasting.

Compilation and Testing

Now let's add a new test that calls and verifies our McWeeny routine. In the UnitTests/Tests/ folder, there is test_chemistry.py.in file. Here if we add a new test, Python will automatically add it to our unit testing setup. The Chemistry solvers have been specially set up to compare to dense python routines.

def test_mcweeny(self):
      '''Test routines to compute the density matrix with TRS4.'''
      self.basic_solver(nt.DensityMatrixSolvers.McWeeny)

Now when one runs make test, this test will be automatically executed. Note that this is a .py.in file. Any changes here require you to recompile the code, even though it is python.

Running every test every time you make a change can get very time consuming. That is why I've provided a per branch testing facility. Go to the UnitTests directory, and copy the file ExampleBranchTest.sh to the exact name of your git branch:

cp ExampleBranchTest.sh McWeeny.sh

Now we'll set the branch test to our new McWeeny test:

BRANCHTEST="test_chemistry.TestChemistry_r.test_mcweeny"

Then you can run just this test by typing:

ctest -R CurrentTest