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

typedef struct qmckl_mo_basis_struct {
  int64_t  mo_num;
  double * restrict coefficient;
  double * restrict coefficient_t;
  double * restrict r_cusp;

  double * restrict mo_vgl;
  double * restrict mo_value;
  qmckl_tensor cusp_param;

  uint64_t  mo_vgl_date;
  uint64_t  mo_value_date;

  int32_t   uninitialized;
  bool      provided;
} qmckl_mo_basis_struct;

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.

qmckl_exit_code qmckl_init_mo_basis(qmckl_context context);

qmckl_exit_code qmckl_init_mo_basis(qmckl_context context) {

  if (qmckl_context_check(context) == QMCKL_NULL_CONTEXT) {
    return false;
  }

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

  ctx->mo_basis.r_cusp = NULL;

  ctx->mo_basis.uninitialized = (1 << 2) - 1;

  return QMCKL_SUCCESS;
}

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.

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.

integer function qmckl_compute_mo_basis_mo_value_doc_f(context, &
     ao_num, mo_num, point_num, &
     coefficient_t, ao_value, mo_value) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context), intent(in)  :: context
  integer*8             , intent(in)  :: ao_num, mo_num
  integer*8             , intent(in)  :: point_num
  double precision      , intent(in)  :: ao_value(ao_num,point_num)
  double precision      , intent(in)  :: coefficient_t(mo_num,ao_num)
  double precision      , intent(out) :: mo_value(mo_num,point_num)
  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_f

qmckl_exit_code qmckl_compute_mo_basis_mo_value (
      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_value,
      double* const mo_value );

qmckl_exit_code qmckl_compute_mo_basis_mo_value_doc (
      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_value,
      double* const mo_value );

qmckl_exit_code
qmckl_compute_mo_basis_mo_value (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_value,
                                 double* const mo_value )
{
#ifdef HAVE_HPC
  return qmckl_compute_mo_basis_mo_value_hpc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#else
  return qmckl_compute_mo_basis_mo_value_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_value, mo_value);
#endif
}

2.2.4 HPC version

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (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_value,
                                     double* const mo_value );

qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc_sp (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_value,
                                        double* const mo_value );
#endif

2.2.4.1 Single-precision

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc_sp (const qmckl_context context,
                                        const int64_t ao_num,
                                        const int64_t mo_num,
                                        const int64_t point_num,
                                        const double* restrict coefficient_t,
                                        const double* restrict ao_value,
                                        double* restrict const mo_value )
{
  assert (context != QMCKL_NULL_CONTEXT);

  float* __attribute__((aligned(64))) coefficient_t_sp = calloc(ao_num*mo_num, sizeof(float));
  if (coefficient_t_sp == NULL) {
    return qmckl_failwith( context,
                           QMCKL_ALLOCATION_FAILED,
                           "qmckl_compute_mo_basis_mo_value_hpc_sp",
                           "coefficient_t_sp");
  };

  IVDEP
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
    for (int64_t i=0 ; i<mo_num*ao_num ; ++i) {
      coefficient_t_sp[i] = (float) coefficient_t[i];
    }

#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
    {
      int64_t* __attribute__((aligned(64))) idx = calloc((size_t) ao_num, sizeof(int64_t));
      float*   __attribute__((aligned(64))) av1 = calloc((size_t) ao_num, sizeof(float));
      assert (idx != NULL);
      assert (av1 != NULL);

#ifdef HAVE_OPENMP
#pragma omp for
#endif
    for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
      double* restrict const vgl1  = &(mo_value[ipoint*mo_num]);
      const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);

      float*   __attribute__((aligned(64))) vgl_sp = calloc((size_t) mo_num, sizeof(float));
      assert (vgl_sp != NULL);

      int64_t nidx=0;
      for (int64_t k=0 ; k<ao_num ; ++k) {
        if (avgl1[k] != 0.) {
          idx[nidx] = k;
          av1[nidx] = (float) avgl1[k];
          ++nidx;
        }
      }

      int64_t n=0;

      for (n=0 ; n < nidx-4 ; n+=4) {
        const float* restrict ck1 = coefficient_t_sp + idx[n  ]*mo_num;
        const float* restrict ck2 = coefficient_t_sp + idx[n+1]*mo_num;
        const float* restrict ck3 = coefficient_t_sp + idx[n+2]*mo_num;
        const float* restrict ck4 = coefficient_t_sp + idx[n+3]*mo_num;

        const float a11 = av1[n  ];
        const float a21 = av1[n+1];
        const float a31 = av1[n+2];
        const float a41 = av1[n+3];

        IVDEP
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
          for (int64_t i=0 ; i<mo_num ; ++i) {
            vgl_sp[i] += ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
          }
      }

      for (int64_t m=n ; m < nidx ; m+=1) {
        const float* restrict ck = coefficient_t_sp + idx[m]*mo_num;
        const float a1 = av1[m];

        IVDEP
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
          for (int64_t i=0 ; i<mo_num ; ++i) {
            vgl_sp[i] += ck[i] * a1;
          }
      }

      IVDEP
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
      for (int64_t i=0 ; i<mo_num ; ++i) {
        vgl1[i] = (double) vgl_sp[i];
      }
      free(vgl_sp);
    }
    free(av1);
    free(idx);
  }
  free(coefficient_t_sp);
  return QMCKL_SUCCESS;
}
#endif

2.2.4.2 Double-precision

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_hpc (const qmckl_context context,
                                     const int64_t ao_num,
                                     const int64_t mo_num,
                                     const int64_t point_num,
                                     const double* restrict coefficient_t,
                                     const double* restrict ao_value,
                                     double* restrict const mo_value )
{
  assert (context != QMCKL_NULL_CONTEXT);
  qmckl_context_struct* const ctx = (qmckl_context_struct*) context;
  assert (ctx != NULL);

  /* Don't compute polynomials when the radial part is zero. */
  const int precision = ctx->numprec.precision;
  const bool single_precision = precision <= 24;

  if (single_precision) {
    return qmckl_compute_mo_basis_mo_value_hpc_sp (context,
                                                   ao_num,
                                                   mo_num,
                                                   point_num,
                                                   coefficient_t,
                                                   ao_value,
                                                   mo_value );
  }

#ifdef HAVE_OPENMP
#pragma omp parallel
#endif
  {
    int64_t* __attribute__((aligned(64))) idx = calloc(ao_num, sizeof(int64_t));
    double*  __attribute__((aligned(64))) av1 = calloc(ao_num, sizeof(double));
    assert (idx != NULL);
    assert (av1 != NULL);

#ifdef HAVE_OPENMP
#pragma omp for
#endif
    for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
      double* restrict const vgl1  = &(mo_value[ipoint*mo_num]);
      const double* restrict avgl1 = &(ao_value[ipoint*ao_num]);

      memset(vgl1, 0, mo_num*sizeof(double));

      int64_t nidx=0;
      int64_t idx[ao_num];
      double  av1[ao_num];

      for (int64_t k=0 ; k<ao_num ; ++k) {
        if (avgl1[k] != 0.) {
          idx[nidx] = k;
          av1[nidx] = avgl1[k];
          ++nidx;
        }
      }

      int64_t n=0;

      for (n=0 ; n < nidx-4 ; n+=4) {
        const double* restrict ck1 = coefficient_t + idx[n  ]*mo_num;
        const double* restrict ck2 = coefficient_t + idx[n+1]*mo_num;
        const double* restrict ck3 = coefficient_t + idx[n+2]*mo_num;
        const double* restrict ck4 = coefficient_t + idx[n+3]*mo_num;

        const double a11 = av1[n  ];
        const double a21 = av1[n+1];
        const double a31 = av1[n+2];
        const double a41 = av1[n+3];

        IVDEP
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
          for (int64_t i=0 ; i<mo_num ; ++i) {
            vgl1[i] += ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
          }
      }

      for (int64_t m=n ; m < nidx ; m+=1) {
        const double* restrict ck = coefficient_t + idx[m]*mo_num;
        const double a1 = av1[m];

        IVDEP
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
          for (int64_t i=0 ; i<mo_num ; ++i) {
            vgl1[i] += ck[i] * a1;
          }
      }
    }
    free(av1);
    free(idx);
  }
  return QMCKL_SUCCESS;
}
#endif

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[mo_num][ao_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