Molecular Orbitals

1. Context
2. Computation

1. Context

The following arrays are stored in the context:

`mo_num`		Number of MOs
`coefficient`	`[mo_num][ao_num]`	MO coefficients
`coefficient_t`	`[ao_num][mo_num]`	Transposed of the Orbital coefficients
`r_cusp`	`[nucl_num]`	Radius of the functions for Cusp adjustments

Computed data:

`cusp_param`	`[nucl_num][4][mo_num]`	Parameters of the functions for Cusp adjustments
`mo_value`	`[point_num][mo_num]`	Value of the MOs at point positions
`mo_vgl`	`[point_num][5][mo_num]`	Value, gradients, Laplacian of the MOs at point positions

1.1. Data structure

The uninitialized integer contains one bit set to one for each initialization function which has not been called. It becomes equal to zero after all initialization functions have been called. The struct is then initialized and provided == true. Some values are initialized by default, and are not concerned by this mechanism.

1.2. Initialization functions

To set the basis set, all the following functions need to be called.

qmckl_exit_code  qmckl_set_mo_basis_mo_num           (qmckl_context context, const int64_t   mo_num);
qmckl_exit_code  qmckl_set_mo_basis_coefficient      (qmckl_context context, const double  * coefficient, const int64_t size_max);
qmckl_exit_code  qmckl_set_mo_basis_r_cusp           (qmckl_context context, const double  * r_cusp, const int64_t size_max);

When the basis set is completely entered, other data structures are computed to accelerate the calculations. The function qmckl_set_mo_basis_coefficients can be called multiple times to override the current coefficients in the context, but the function qmckl_set_mo_basis_mo_num can only be called once in a context.

1.3. Cusp adjsutment functions

To activate the cusp adjustment, the user must enter the radius of the fitting function for each atom.

This function requires the computation of the value and gradients of the \(s\) AOs at the distance equal to the radius, and the values of the non-\(s\) AOs at the center.

1.4. Access functions

When all the data for the AOs have been provided, the following function returns true.

bool qmckl_mo_basis_provided (const qmckl_context context);

1.4.1. Fortran interfaces

1.5. Update

It may be desirable to remove certain molecular orbitals (MOs) that do not significantly contribute to the wave function. In particular, in a single determinant calculation, the virtual MOs can be removed as they do not participate in the ground state configuration.

To select a subset of MOs that will be kept, an array of integers of size mo_num can be created. If the integer corresponding to an MO is zero, that MO is dropped and will not be included in the calculation. If the integer is non-zero, the MO will be kept.

qmckl_exit_code
qmckl_mo_basis_select_mo (const qmckl_context context,
                          const int32_t* keep,
                          const int64_t size_max);

1.5.1. Fortran interface

2. Computation

2.1. Parameters of the cusp-correction functions

2.2. Computation of MOs: values only

2.2.1. Get

qmckl_exit_code
qmckl_get_mo_basis_mo_value(qmckl_context context,
                            double* const mo_value,
                            const int64_t size_max);

Uses the given array to compute the values.

qmckl_exit_code
qmckl_get_mo_basis_mo_value_inplace (qmckl_context context,
                                     double* const mo_value,
                                     const int64_t size_max);

2.2.2. Provide

qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context);

qmckl_exit_code qmckl_provide_mo_basis_mo_value(qmckl_context context)
{

  qmckl_exit_code rc = QMCKL_SUCCESS;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_provide_mo_basis_mo_value",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  if (!ctx->mo_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_provide_mo_basis_mo_value",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->mo_basis.mo_value_date) {

    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    mem_info.size = ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);

    if (ctx->mo_basis.mo_value != NULL) {
      qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
      rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_value, &mem_info_test);

      /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
         memory was not allocated with qmckl_malloc */

      if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
        rc = qmckl_free(context, ctx->mo_basis.mo_value);
        assert (rc == QMCKL_SUCCESS);
        ctx->mo_basis.mo_value = NULL;
      }
    }

    /* Allocate array */
    if (ctx->mo_basis.mo_value == NULL) {

      double* mo_value = (double*) qmckl_malloc(context, mem_info);

      if (mo_value == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_mo_basis_mo_value",
                               NULL);
      }
      ctx->mo_basis.mo_value = mo_value;
    }

    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->mo_basis.mo_value_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}

2.2.3. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`ao_num`	`int64_t`	in	Number of AOs
`mo_num`	`int64_t`	in	Number of MOs
`point_num`	`int64_t`	in	Number of points
`coefficient_t`	`double[mo_num][ao_num]`	in	Transpose of the AO to MO transformation matrix
`ao_value`	`double[point_num][ao_num]`	in	Value of the AOs
`mo_value`	`double[point_num][mo_num]`	out	Value of the MOs

The matrix of AO values is very sparse, so we use a sparse-dense matrix multiplication instead of a dgemm, as exposed in https://dx.doi.org/10.1007/978-3-642-38718-0_14.

function qmckl_compute_mo_basis_mo_value_doc(context, &
     ao_num, mo_num, point_num, &
     coefficient_t, ao_value, mo_value) &
     result(info) bind(C)
  use qmckl_constants
  implicit none
  integer (qmckl_context) , intent(in)  , value :: context
  integer (c_int64_t) , intent(in)  , value :: ao_num
  integer (c_int64_t) , intent(in)  , value :: mo_num
  integer (c_int64_t) , intent(in)  , value :: point_num
  real    (c_double ) , intent(in)          :: coefficient_t(mo_num,ao_num)
  real    (c_double ) , intent(in)          :: ao_value(ao_num,point_num)
  real    (c_double ) , intent(out)         :: mo_value(mo_num,point_num)

  integer(qmckl_exit_code) :: info

  integer*8 :: j,k

  info = QMCKL_SUCCESS

  do j=1,point_num
     mo_value(:,j) = 0.d0
     do k=1,ao_num
        if (ao_value(k,j) == 0.d0) cycle
        mo_value(:,j) = mo_value(:,j) + coefficient_t(:,k) * ao_value(k,j)
     end do
  end do

end function qmckl_compute_mo_basis_mo_value_doc

2.3. Computation of MOs: values, gradient, Laplacian

2.3.1. Get

qmckl_exit_code
qmckl_get_mo_basis_mo_vgl(qmckl_context context,
                          double* const mo_vgl,
                          const int64_t size_max);

Uses the given array to compute the VGL.

qmckl_exit_code
qmckl_get_mo_basis_mo_vgl_inplace (qmckl_context context,
                                   double* const mo_vgl,
                                   const int64_t size_max);

2.3.2. Provide

qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context);

qmckl_exit_code qmckl_provide_mo_basis_mo_vgl(qmckl_context context)
{

  qmckl_exit_code rc = QMCKL_SUCCESS;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_provide_mo_basis_mo_vgl",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  if (!ctx->mo_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_provide_mo_basis_mo_vgl",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->mo_basis.mo_vgl_date) {

    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    mem_info.size = 5 * ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);

    if (ctx->mo_basis.mo_vgl != NULL) {
      qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
      rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_vgl, &mem_info_test);

      /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
         memory was not allocated with qmckl_malloc */

      if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
        rc = qmckl_free(context, ctx->mo_basis.mo_vgl);
        assert (rc == QMCKL_SUCCESS);
        ctx->mo_basis.mo_vgl = NULL;
      }
    }

    /* Allocate array */
    if (ctx->mo_basis.mo_vgl == NULL) {

      double* mo_vgl = (double*) qmckl_malloc(context, mem_info);

      if (mo_vgl == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_mo_basis_mo_vgl",
                               NULL);
      }
      ctx->mo_basis.mo_vgl = mo_vgl;
    }

    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->mo_basis.mo_vgl_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}

2.3.3. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`ao_num`	`int64_t`	in	Number of AOs
`mo_num`	`int64_t`	in	Number of MOs
`point_num`	`int64_t`	in	Number of points
`coefficient_t`	`double[ao_num][mo_num]`	in	Transpose of the AO to MO transformation matrix
`ao_vgl`	`double[point_num][5][ao_num]`	in	Value, gradients and Laplacian of the AOs
`mo_vgl`	`double[point_num][5][mo_num]`	out	Value, gradients and Laplacian of the MOs

The matrix of AO values is very sparse, so we use a sparse-dense matrix multiplication instead of a dgemm, as exposed in https://dx.doi.org/10.1007/978-3-642-38718-0_14.

function qmckl_compute_mo_basis_mo_vgl_doc(context, &
     ao_num, mo_num, point_num, &
     coefficient_t, ao_vgl, mo_vgl) &
     result(info) bind(C)
  use qmckl_constants
  implicit none
  integer(qmckl_context), intent(in), value  :: context
  integer (c_int64_t) , intent(in)  , value :: ao_num
  integer (c_int64_t) , intent(in)  , value :: mo_num
  integer (c_int64_t) , intent(in)  , value :: point_num
  real    (c_double ) , intent(in)          :: coefficient_t(mo_num,ao_num)
  real    (c_double ) , intent(in)          :: ao_vgl(ao_num,5,point_num)
  real    (c_double ) , intent(out)         :: mo_vgl(mo_num,5,point_num)

  integer(qmckl_exit_code) :: info

  integer*8 :: i,j,k
  double precision :: c1, c2, c3, c4, c5

  info = QMCKL_SUCCESS

  do j=1,point_num
     mo_vgl(:,:,j) = 0.d0
     do k=1,ao_num
        if (ao_vgl(k,1,j) /= 0.d0) then
           c1 = ao_vgl(k,1,j)
           c2 = ao_vgl(k,2,j)
           c3 = ao_vgl(k,3,j)
           c4 = ao_vgl(k,4,j)
           c5 = ao_vgl(k,5,j)
           do i=1,mo_num
              mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
              mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
              mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
              mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
              mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
           end do
        end if
     end do
  end do

end function qmckl_compute_mo_basis_mo_vgl_doc

2.4. Computation of MOs: Hessian

2.4.1. Get

qmckl_exit_code
qmckl_get_mo_basis_mo_hessian(qmckl_context context,
                          double* const mo_hessian,
                          const int64_t size_max);

2.4.2. Provide

qmckl_exit_code qmckl_provide_mo_basis_mo_hessian(qmckl_context context);

qmckl_exit_code qmckl_provide_mo_basis_mo_hessian(qmckl_context context)
{

  qmckl_exit_code rc = QMCKL_SUCCESS;

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return qmckl_failwith( context,
                           QMCKL_INVALID_CONTEXT,
                           "qmckl_provide_mo_basis_mo_hessian",
                           NULL);
  }

  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  if (!ctx->mo_basis.provided) {
    return qmckl_failwith( context,
                           QMCKL_NOT_PROVIDED,
                           "qmckl_provide_mo_basis_mo_hessian",
                           NULL);
  }

  /* Compute if necessary */
  if (ctx->point.date > ctx->mo_basis.mo_hessian_date) {

    qmckl_memory_info_struct mem_info = qmckl_memory_info_struct_zero;
    mem_info.size = 3 * 3 * ctx->mo_basis.mo_num * ctx->point.num * sizeof(double);

    if (ctx->mo_basis.mo_hessian != NULL) {
      qmckl_memory_info_struct mem_info_test = qmckl_memory_info_struct_zero;
      rc = qmckl_get_malloc_info(context, ctx->mo_basis.mo_hessian, &mem_info_test);

      /* if rc != QMCKL_SUCCESS, we are maybe in an _inplace function because the
         memory was not allocated with qmckl_malloc */

      if ((rc == QMCKL_SUCCESS) && (mem_info_test.size != mem_info.size)) {
        rc = qmckl_free(context, ctx->mo_basis.mo_hessian);
        assert (rc == QMCKL_SUCCESS);
        ctx->mo_basis.mo_hessian = NULL;
      }
    }

    /* Allocate array */
    if (ctx->mo_basis.mo_hessian == NULL) {

      double* mo_hessian = (double*) qmckl_malloc(context, mem_info);

      if (mo_hessian == NULL) {
        return qmckl_failwith( context,
                               QMCKL_ALLOCATION_FAILED,
                               "qmckl_mo_basis_mo_hessian",
                               NULL);
      }
      ctx->mo_basis.mo_hessian = mo_hessian;
    }

    rc = qmckl_provide_ao_basis_ao_hessian(context);
    if (rc != QMCKL_SUCCESS) {
      return qmckl_failwith( context,
                             QMCKL_NOT_PROVIDED,
                             "qmckl_ao_basis",
                             NULL);
    }

    if (ctx->mo_basis.r_cusp == NULL) {
      /* No cusp correction */
      rc = qmckl_compute_mo_basis_mo_hessian_doc_f(context,
                                         ctx->ao_basis.ao_num,
                                         ctx->mo_basis.mo_num,
                                         ctx->point.num,
                                         ctx->mo_basis.coefficient_t,
                                         ctx->ao_basis.ao_hessian,
                                         ctx->mo_basis.mo_hessian);
    } else {
      rc = qmckl_provide_en_distance(context);
      if (rc != QMCKL_SUCCESS) {
        return qmckl_failwith( context,
                               QMCKL_NOT_PROVIDED,
                               "qmckl_electron_en_distance",
                               NULL);
      }
      rc = qmckl_compute_mo_basis_mo_hessian_cusp(context,
                                              ctx->nucleus.num,
                                              ctx->ao_basis.ao_num,
                                              ctx->mo_basis.mo_num,
                                              ctx->point.num,
                                              ctx->ao_basis.ao_nucl,
                                              ctx->ao_basis.ao_ang_mom,
                                              ctx->electron.en_distance,
                                              ctx->nucleus.coord,
                                              ctx->point.coord,
                                              ctx->mo_basis.r_cusp,
                                              ctx->mo_basis.cusp_param,
                                              ctx->mo_basis.coefficient_t,
                                              ctx->ao_basis.ao_hessian,
                                              ctx->mo_basis.mo_hessian);
    }

    if (rc != QMCKL_SUCCESS) {
      return rc;
    }

    ctx->mo_basis.mo_hessian_date = ctx->date;
  }

  return QMCKL_SUCCESS;
}

2.4.3. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`ao_num`	`int64_t`	in	Number of AOs
`mo_num`	`int64_t`	in	Number of MOs
`point_num`	`int64_t`	in	Number of points
`coefficient_t`	`double[mo_num][ao_num]`	in	Transpose of the AO to MO transformation matrix
`ao_hessian`	`double[3][point_num][4][ao_num]`	in	Hessian of the AOs
`mo_hessian`	`double[3][point_num][3][mo_num]`	out	Hessian of the MOs

integer function qmckl_compute_mo_basis_mo_hessian_doc_f(context, &
     ao_num, mo_num, point_num, &
     coefficient_t, ao_hessian, mo_hessian) &
     bind(C) result(info)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer (qmckl_context), intent(in), value  :: context
  integer (c_int64_t)    , intent(in), value  :: ao_num, mo_num
  integer (c_int64_t)    , intent(in), value  :: point_num
  real    (c_double )    , intent(in)         :: ao_hessian(ao_num,4,point_num,3)
  real    (c_double )    , intent(in)         :: coefficient_t(mo_num,ao_num)
  real    (c_double )    , intent(out)        :: mo_hessian(mo_num,3,point_num,3)
  integer*8 :: i,j,k, kk, jj

  info = QMCKL_SUCCESS

  do jj=1,3
    do j=1,point_num
      mo_hessian(:,:,j,jj) = 0.d0
      do k=1,ao_num
          do kk=1,3
            do i=1,mo_num
              mo_hessian(i,kk,j,jj) = mo_hessian(i,kk,j,jj) + coefficient_t(i,k) * ao_hessian(k,kk,j,jj)
            end do
          end do
      end do
    end do
  end do

end function qmckl_compute_mo_basis_mo_hessian_doc_f

qmckl_exit_code qmckl_compute_mo_basis_mo_hessian_doc_f (
      const qmckl_context context,
      const int64_t ao_num,
      const int64_t mo_num,
      const int64_t point_num,
      const double* coefficient_t,
      const double* ao_hessian,
      double* const mo_hessian );

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`nucl_num`	`int64_t`	in	Number of nuclei
`ao_num`	`int64_t`	in	Number of AOs
`mo_num`	`int64_t`	in	Number of MOs
`point_num`	`int64_t`	in	Number of points
`ao_nucl`	`int64_t[ao_num]`	in	Nucleus on which the AO is centered
`ao_ang_mom`	`int32_t[ao_num]`	in	Angular momentum of the shell
`en_distance`	`double[point_num][nucl_num]`	in	Electron-nucleus distances
`nucl_coord`	`double[3][nucl_num]`	in	Nuclear coordinates
`point_coord`	`double[3][point_num]`	in	Electron coordinates
`r_cusp`	`double[nucl_num]`	in	Cusp-adjustment radius
`cusp_param`	`double[nucl_num][4][mo_num]`	in	Cusp-adjustment parameters
`coefficient_t`	`double[mo_num][ao_num]`	in	Transpose of the AO to MO transformation matrix
`ao_hessian`	`double[3][point_num][4][ao_num]`	in	Hessian of the AOs
`mo_hessian`	`double[3][point_num][3][mo_num]`	out	Hessian of the MOs

integer function qmckl_compute_mo_basis_mo_hessian_cusp_doc_f(context, &
     nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
     nucl_coord, point_coord, r_cusp, cusp_param, coefficient_t, ao_hessian, mo_hessian) &
     bind(C) result(info)
  use, intrinsic :: iso_c_binding
  use qmckl
  implicit none
  integer(qmckl_context) , intent(in), value  :: context
  integer (c_int64_t)    , intent(in), value  :: nucl_num, ao_num, mo_num, point_num
  integer (c_int64_t)    , intent(in)         :: ao_nucl(ao_num)
  integer (c_int32_t)    , intent(in)         :: ao_ang_mom(ao_num)
  real    (c_double )    , intent(in)         :: en_distance(nucl_num, point_num)
  real    (c_double )    , intent(in)         :: nucl_coord(nucl_num,3)
  real    (c_double )    , intent(in)         :: point_coord(point_num,3)
  real    (c_double )    , intent(in)         :: r_cusp(nucl_num)
  real    (c_double )    , intent(in)         :: cusp_param(mo_num,4,nucl_num)
  real    (c_double )    , intent(in)         :: coefficient_t(mo_num,ao_num)
  real    (c_double )    , intent(in)         :: ao_hessian(ao_num,4,point_num,3)
  real    (c_double )    , intent(out)        :: mo_hessian(mo_num,3,point_num,3)
  integer*8 :: i,j,k, inucl, kk, jj
  double precision :: r, r_inv, r_vec(3)

  do jj=1,3
    do j=1,point_num

     ! Initial contribution of the MO
     mo_hessian(:,:,j,jj) = 0.d0
     do k=1,ao_num
        if (ao_hessian(k,1,j,jj) /= 0.d0) then
          inucl = ao_nucl(k)+1
          if ( (en_distance(inucl,j) > r_cusp(inucl)) .or. (ao_ang_mom(k) > 0) ) then
            do kk=1,3
              do i=1,mo_num
                  mo_hessian(i,kk,j,jj) = mo_hessian(i,kk,j,jj) + coefficient_t(i,k) * ao_hessian(k,kk,j,jj)
              end do
            end do
          end if
       end if
     end do

      ! Cusp adjustment
      do inucl=1,nucl_num
        r = en_distance(inucl,j)
        if (r < r_cusp(inucl)) then

            r_vec(1:3) = point_coord(j,1:3) - nucl_coord(inucl,1:3)
            r_inv = 1.d0/r
            do kk=1,3
              do i=1,mo_num

                mo_hessian(i,kk,j,jj) = mo_hessian(i,kk,j,jj) - &
                    cusp_param(i,2,inucl) * r_vec(kk) * r_vec(jj) * r_inv**3 + &
                    cusp_param(i,4,inucl) * 3 * r_vec(kk) * r_vec(jj) * r_inv

                if (kk == jj) then
                  mo_hessian(i,kk,j,jj) = mo_hessian(i,kk,j,jj) + &
                      cusp_param(i,2,inucl) * r_inv + &
                      cusp_param(i,3,inucl) * 2 + &
                      cusp_param(i,4,inucl) * 3 * r
                end if
              end do
            end do
        end if
      end do ! inucl
    end do
  end do ! jj
  info = QMCKL_SUCCESS

end function qmckl_compute_mo_basis_mo_hessian_cusp_doc_f

qmckl_exit_code qmckl_compute_mo_basis_mo_hessian_cusp_doc_f (
      const qmckl_context context,
      const int64_t nucl_num,
      const int64_t ao_num,
      const int64_t mo_num,
      const int64_t point_num,
      const int64_t* ao_nucl,
      const int32_t* ao_ang_mom,
      const double* en_distance,
      const double* nucl_coord,
      const double* point_coord,
      const double* r_cusp,
      const double* cusp_param,
      const double* coefficient_t,
      const double* ao_hessian,
      double* const mo_hessian );

qmckl_exit_code qmckl_compute_mo_basis_mo_hessian_cusp (
      const qmckl_context context,
      const int64_t nucl_num,
      const int64_t ao_num,
      const int64_t mo_num,
      const int64_t point_num,
      const int64_t* ao_nucl,
      const int32_t* ao_ang_mom,
      const double* en_distance,
      const qmckl_matrix nucl_coord,
      const qmckl_matrix point_coord,
      const double* r_cusp,
      const qmckl_tensor cusp_param,
      const double* coefficient_t,
      const double* ao_hessian,
      double* const mo_hessian );

qmckl_exit_code
qmckl_compute_mo_basis_mo_hessian_cusp (const qmckl_context context,
                                    const int64_t nucl_num,
                                    const int64_t ao_num,
                                    const int64_t mo_num,
                                    const int64_t point_num,
                                    const int64_t* ao_nucl,
                                    const int32_t* ao_ang_mom,
                                    const double* en_distance,
                                    const qmckl_matrix nucl_coord_matrix,
                                    const qmckl_matrix point_coord_matrix,
                                    const double* r_cusp,
                                    const qmckl_tensor cusp_param_tensor,
                                    const double* coefficient_t,
                                    const double* ao_hessian,
                                    double* const mo_hessian )
{
  qmckl_exit_code rc;


  double * nucl_coord  = qmckl_alloc_double_of_matrix(context, nucl_coord_matrix);
  double * point_coord = qmckl_alloc_double_of_matrix(context, point_coord_matrix);
  double * cusp_param  = qmckl_alloc_double_of_tensor(context, cusp_param_tensor);

  rc = qmckl_compute_mo_basis_mo_hessian_cusp_doc_f (context, nucl_num, ao_num, mo_num, point_num,
                                               ao_nucl, ao_ang_mom, en_distance, nucl_coord,
                                               point_coord, r_cusp, cusp_param, coefficient_t,
                                               ao_hessian, mo_hessian );

  qmckl_free(context, nucl_coord);
  qmckl_free(context, point_coord);
  qmckl_free(context, cusp_param);

  return rc;
}

2.5. Computation of cusp-corrected MOs: values only

2.5.1. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`nucl_num`	`int64_t`	in	Number of nuclei
`ao_num`	`int64_t`	in	Number of AOs
`mo_num`	`int64_t`	in	Number of MOs
`point_num`	`int64_t`	in	Number of points
`ao_nucl`	`int64_t[ao_num]`	in	Nucleus on which the AO is centered
`ao_ang_mom`	`int32_t[ao_num]`	in	Angular momentum of the shell
`en_distance`	`double[point_num][nucl_num]`	in	Electron-nucleus distances
`r_cusp`	`double[nucl_num]`	in	Cusp-adjustment radius
`cusp_param`	`double[nucl_num][4][mo_num]`	in	Cusp-adjustment parameters
`coefficient_t`	`double[ao_num][mo_num]`	in	Transpose of the AO to MO transformation matrix
`ao_value`	`double[point_num][ao_num]`	in	Value of the AOs
`mo_value`	`double[point_num][mo_num]`	out	Cusp correction for the values of the MOs

function qmckl_compute_mo_basis_mo_value_cusp_doc(context, &
     nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
     r_cusp, cusp_param, coefficient_t, ao_value, mo_value) &
     result(info) bind(C)
  use qmckl_constants
  implicit none
  integer(qmckl_context), intent(in), value  :: context
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: ao_num
  integer (c_int64_t) , intent(in)  , value :: mo_num
  integer (c_int64_t) , intent(in)  , value :: point_num
  integer (c_int64_t) , intent(in)          :: ao_nucl(ao_num)
  integer (c_int32_t) , intent(in)          :: ao_ang_mom(ao_num)
  real    (c_double ) , intent(in)          :: en_distance(nucl_num,point_num)
  real    (c_double ) , intent(in)          :: r_cusp(nucl_num)
  real    (c_double ) , intent(in)          :: cusp_param(mo_num,4,nucl_num)
  real    (c_double ) , intent(in)          :: coefficient_t(mo_num,ao_num)
  real    (c_double ) , intent(in)          :: ao_value(ao_num,point_num)
  real    (c_double ) , intent(out)         :: mo_value(mo_num,point_num)

  integer(qmckl_exit_code) :: info

  integer*8 :: i, j, k, inucl
  double precision :: r

  info = QMCKL_SUCCESS

  do i=1,point_num
     mo_value(:,i) = 0.d0
     do k=1,ao_num
        if (ao_value(k,i) == 0.d0) cycle
        inucl = ao_nucl(k)+1
        if ( (en_distance(inucl,i) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
        mo_value(:,i) = mo_value(:,i) + coefficient_t(:,k) * ao_value(k,i)
     end do ! k

     do inucl=1,nucl_num
        r = en_distance(inucl,i)
        if (r > r_cusp(inucl)) cycle

        do j=1,mo_num
           mo_value(j,i) = mo_value(j,i) + &
                cusp_param(j,1,inucl) + r*(cusp_param(j,2,inucl) + r*(  &
                cusp_param(j,3,inucl) + r* cusp_param(j,4,inucl)     ))
        enddo
     enddo ! inucl
  enddo ! i

end function qmckl_compute_mo_basis_mo_value_cusp_doc

2.6. Computation of cusp-corrected MOs: values, gradient, Laplacian

2.6.1. Compute

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`nucl_num`	`int64_t`	in	Number of nuclei
`ao_num`	`int64_t`	in	Number of AOs
`mo_num`	`int64_t`	in	Number of MOs
`point_num`	`int64_t`	in	Number of points
`ao_nucl`	`int64_t[ao_num]`	in	Nucleus on which the AO is centered
`ao_ang_mom`	`int32_t[ao_num]`	in	Angular momentum of the shell
`en_distance`	`double[point_num][nucl_num]`	in	Electron-nucleus distances
`nucl_coord`	`double[3][nucl_num]`	in	Nuclear coordinates
`point_coord`	`double[3][point_num]`	in	Electron coordinates
`r_cusp`	`double[nucl_num]`	in	Cusp-adjustment radius
`cusp_param`	`double[nucl_num][4][mo_num]`	in	Cusp-adjustment parameters
`coefficient_t`	`double[ao_num][mo_num]`	in	Transpose of the AO to MO transformation matrix
`ao_vgl`	`double[point_num][5][ao_num]`	in	Value, gradients and Laplacian of the AOs
`mo_vgl`	`double[point_num][5][mo_num]`	out	Value, gradients and Laplacian of the MOs

function qmckl_compute_mo_basis_mo_vgl_cusp_doc(context, &
     nucl_num, ao_num, mo_num, point_num, ao_nucl, ao_ang_mom, en_distance, &
     nucl_coord, point_coord, r_cusp, cusp_param, coefficient_t, ao_vgl, mo_vgl) &
     result(info) bind(C)
  use qmckl_constants
  implicit none
  integer(qmckl_context), intent(in), value  :: context
  integer (c_int64_t) , intent(in)  , value :: nucl_num
  integer (c_int64_t) , intent(in)  , value :: ao_num
  integer (c_int64_t) , intent(in)  , value :: mo_num
  integer (c_int64_t) , intent(in)  , value :: point_num
  integer (c_int64_t) , intent(in)          :: ao_nucl(ao_num)
  integer (c_int32_t) , intent(in)          :: ao_ang_mom(ao_num)
  real    (c_double ) , intent(in)          :: en_distance(nucl_num,point_num)
  real    (c_double ) , intent(in)          :: nucl_coord(nucl_num,3)
  real    (c_double ) , intent(in)          :: point_coord(point_num,3)
  real    (c_double ) , intent(in)          :: r_cusp(nucl_num)
  real    (c_double ) , intent(in)          :: cusp_param(mo_num,4,nucl_num)
  real    (c_double ) , intent(in)          :: coefficient_t(mo_num,ao_num)
  real    (c_double ) , intent(in)          :: ao_vgl(ao_num,5,point_num)
  real    (c_double ) , intent(out)         :: mo_vgl(mo_num,5,point_num)

  integer(qmckl_exit_code) :: info

  integer*8 :: i,j,k, inucl
  double precision :: c1, c2, c3, c4, c5
  double precision :: r, r_inv, r_vec(3)


  do j=1,point_num

     ! Initial contribution of the MO
     mo_vgl(:,:,j) = 0.d0
     do k=1,ao_num
        if (ao_vgl(k,1,j) /= 0.d0) then
          inucl = ao_nucl(k)+1
          if ( (en_distance(inucl,j) > r_cusp(inucl)) .or. (ao_ang_mom(k) > 0) ) then
             c1 = ao_vgl(k,1,j)
             c2 = ao_vgl(k,2,j)
             c3 = ao_vgl(k,3,j)
             c4 = ao_vgl(k,4,j)
             c5 = ao_vgl(k,5,j)
             do i=1,mo_num
                mo_vgl(i,1,j) = mo_vgl(i,1,j) + coefficient_t(i,k) * c1
                mo_vgl(i,2,j) = mo_vgl(i,2,j) + coefficient_t(i,k) * c2
                mo_vgl(i,3,j) = mo_vgl(i,3,j) + coefficient_t(i,k) * c3
                mo_vgl(i,4,j) = mo_vgl(i,4,j) + coefficient_t(i,k) * c4
                mo_vgl(i,5,j) = mo_vgl(i,5,j) + coefficient_t(i,k) * c5
             end do
          end if
       end if
     end do

    ! Cusp adjustment
    do inucl=1,nucl_num
       r = en_distance(inucl,j)
       if (r < r_cusp(inucl)) then

          r_vec(1:3) = point_coord(j,1:3) - nucl_coord(inucl,1:3)
          r_inv = 1.d0/r

          do i=1,mo_num
             mo_vgl(i,1,j) = mo_vgl(i,1,j) +  &
                  cusp_param(i,1,inucl) + r*( &
                  cusp_param(i,2,inucl) + r*( &
                  cusp_param(i,3,inucl) + r*( &
                  cusp_param(i,4,inucl)    )))

             c1 = r_inv * cusp_param(i,2,inucl) + 2.d0*cusp_param(i,3,inucl) +  &
                  r * 3.d0 * cusp_param(i,4,inucl)

             mo_vgl(i,2,j) = mo_vgl(i,2,j) + r_vec(1) * c1
             mo_vgl(i,3,j) = mo_vgl(i,3,j) + r_vec(2) * c1
             mo_vgl(i,4,j) = mo_vgl(i,4,j) + r_vec(3) * c1

             mo_vgl(i,5,j) = mo_vgl(i,5,j) +         &
                  2.d0*cusp_param(i,2,inucl)*r_inv + &
                  6.d0*cusp_param(i,3,inucl) +       &
                  12.d0*cusp_param(i,4,inucl)*r

          end do
       end if
    end do ! inucl
  end do
  info = QMCKL_SUCCESS

end function qmckl_compute_mo_basis_mo_vgl_cusp_doc

2.7. Rescaling of MO coefficients

When evaluating Slater determinants, the value of the determinants may get out of the range of double precision. A simple fix is to rescale the MO coefficients to put back the determinants in the correct range.

qmckl_exit_code
qmckl_mo_basis_rescale(qmckl_context context,
                          const double scaling_factor);