UP | HOME

Molecular Orbitals

Table of Contents

1 Context

The following arrays are stored in the context:

mo_num   Number of MOs
coefficient [mo_num][ao_num] Orbital coefficients
coefficient_t [ao_num][mo_num] Transposed of the Orbital coefficients

Computed data:

mo_vgl [point_num][5][mo_num] Value, gradients, Laplacian of the MOs at point positions
mo_vgl_date uint64_t Late modification date of 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 mo_vgl;
  uint64_t  mo_vgl_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* const) context;
  assert (ctx != NULL);

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

  return QMCKL_SUCCESS;
}

1.2 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.2.1 Fortran interfaces

interface
  integer(c_int32_t) function qmckl_get_mo_basis_mo_num (context, &
       mo_num) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    integer (c_int64_t) , intent(out)         :: mo_num
  end function qmckl_get_mo_basis_mo_num
end interface

interface
  integer(c_int32_t) function qmckl_get_mo_basis_coefficient(context, &
       coefficient, size_max) bind(C)
    use, intrinsic :: iso_c_binding
    import
    implicit none
    integer (c_int64_t) , intent(in)  , value :: context
    double precision, intent(out)             :: coefficient(*)
    integer (c_int64_t) , intent(in)  , value :: size_max
  end function qmckl_get_mo_basis_coefficient
end interface

1.3 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);

When the basis set is completely entered, other data structures are computed to accelerate the calculations.

2 Computation

2.1 Computation of MOs

2.1.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.1.2 Provide

2.1.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
coef_normalized_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

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_vgl_doc_f(context, &
     ao_num, mo_num, point_num, &
     coef_normalized_t, ao_vgl, mo_vgl) &
     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_vgl(ao_num,5,point_num)
  double precision      , intent(in)  :: coef_normalized_t(mo_num,ao_num)
  double precision      , intent(out) :: mo_vgl(mo_num,5,point_num)
  integer*8 :: i,j,k
  double precision :: c1, c2, c3, c4, c5

  integer*8 :: LDA, LDB, LDC

  info = QMCKL_SUCCESS
  if (.True.)  then    ! fast algorithm
     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) + coef_normalized_t(i,k) * c1
                 mo_vgl(i,2,j) = mo_vgl(i,2,j) + coef_normalized_t(i,k) * c2
                 mo_vgl(i,3,j) = mo_vgl(i,3,j) + coef_normalized_t(i,k) * c3
                 mo_vgl(i,4,j) = mo_vgl(i,4,j) + coef_normalized_t(i,k) * c4
                 mo_vgl(i,5,j) = mo_vgl(i,5,j) + coef_normalized_t(i,k) * c5
              end do
           end if
        end do
     end do

  else ! dgemm

    LDA = size(coef_normalized_t,1)
    LDB = size(ao_vgl,1) 
    LDC = size(mo_vgl,1)

    info = qmckl_dgemm(context,'N', 'N', mo_num, point_num*5_8, ao_num*1_8, 1.d0,     &
                                    coef_normalized_t, LDA, ao_vgl, LDB, &
                                    0.d0, mo_vgl, LDC)

  end if

end function qmckl_compute_mo_basis_mo_vgl_doc_f

qmckl_exit_code qmckl_compute_mo_basis_mo_vgl (
      const qmckl_context context,
      const int64_t ao_num,
      const int64_t mo_num,
      const int64_t point_num,
      const double* coef_normalized_t,
      const double* ao_vgl,
      double* const mo_vgl ); 

qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_doc (
      const qmckl_context context,
      const int64_t ao_num,
      const int64_t mo_num,
      const int64_t point_num,
      const double* coef_normalized_t,
      const double* ao_vgl,
      double* const mo_vgl ); 

qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl (const qmckl_context context,
                            const int64_t ao_num,
                            const int64_t mo_num,
                            const int64_t point_num,
                            const double* coef_normalized_t,
                            const double* ao_vgl,
                            double* const mo_vgl )
{
#ifdef HAVE_HPC
  return qmckl_compute_mo_basis_mo_vgl_hpc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl);
#else
  return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coef_normalized_t, ao_vgl, mo_vgl);
#endif
}

2.1.4 HPC version

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
                                   const int64_t ao_num,
                                   const int64_t mo_num,
                                   const int64_t point_num,
                                   const double* coef_normalized_t,
                                   const double* ao_vgl,
                                   double* const mo_vgl );
#endif

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_hpc (const qmckl_context context,
                                   const int64_t ao_num,
                                   const int64_t mo_num,
                                   const int64_t point_num,
                                   const double* restrict coef_normalized_t,
                                   const double* restrict ao_vgl,
                                   double* restrict const mo_vgl )
{
#ifdef HAVE_OPENMP
  #pragma omp parallel for
#endif
  for (int64_t ipoint=0 ; ipoint < point_num ; ++ipoint) {
    double* restrict const vgl1 = &(mo_vgl[ipoint*5*mo_num]);
    double* restrict const vgl2 =  vgl1 + mo_num;
    double* restrict const vgl3 =  vgl1 + (mo_num << 1);
    double* restrict const vgl4 =  vgl1 + (mo_num << 1) + mo_num;
    double* restrict const vgl5 =  vgl1 + (mo_num << 2);

    const double* restrict avgl1 = &(ao_vgl[ipoint*5*ao_num]);
    const double* restrict avgl2 = avgl1 + ao_num;
    const double* restrict avgl3 = avgl1 + (ao_num << 1);
    const double* restrict avgl4 = avgl1 + (ao_num << 1) + ao_num;
    const double* restrict avgl5 = avgl1 + (ao_num << 2);

    for (int64_t i=0 ; i<mo_num ; ++i) {
      vgl1[i] = 0.;
      vgl2[i] = 0.;
      vgl3[i] = 0.;
      vgl4[i] = 0.;
      vgl5[i] = 0.;
    }

    int64_t nidx=0;
    int64_t idx[ao_num];
    double  av1[ao_num];
    double  av2[ao_num];
    double  av3[ao_num];
    double  av4[ao_num];
    double  av5[ao_num];
    for (int64_t k=0 ; k<ao_num ; ++k) {
      const double* restrict ck1 = coef_normalized_t + k*mo_num;
      if (avgl1[k] != 0.) {
        idx[nidx] = k;
        av1[nidx] = avgl1[k];
        av2[nidx] = avgl2[k];
        av3[nidx] = avgl3[k];
        av4[nidx] = avgl4[k];
        av5[nidx] = avgl5[k];
        ++nidx;
      }
    }

    int64_t n;

    for (n=0 ; n < nidx-4 ; n+=4) {
      int64_t k = idx[n];
      const double* restrict ck1 = coef_normalized_t + idx[n  ]*mo_num;
      const double* restrict ck2 = coef_normalized_t + idx[n+1]*mo_num;
      const double* restrict ck3 = coef_normalized_t + idx[n+2]*mo_num;
      const double* restrict ck4 = coef_normalized_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];

      const double a12 = av2[n  ];
      const double a22 = av2[n+1];
      const double a32 = av2[n+2];
      const double a42 = av2[n+3];

      const double a13 = av3[n  ];
      const double a23 = av3[n+1];
      const double a33 = av3[n+2];
      const double a43 = av3[n+3];

      const double a14 = av4[n  ];
      const double a24 = av4[n+1];
      const double a34 = av4[n+2];
      const double a44 = av4[n+3];

      const double a15 = av5[n  ];
      const double a25 = av5[n+1];
      const double a35 = av5[n+2];
      const double a45 = av5[n+3];

#ifdef HAVE_OPENMP
#pragma omp simd
#endif
      for (int64_t i=0 ; i<mo_num ; ++i) {
        vgl1[i] = vgl1[i] + ck1[i] * a11 + ck2[i] * a21 + ck3[i] * a31 + ck4[i] * a41;
        vgl2[i] = vgl2[i] + ck1[i] * a12 + ck2[i] * a22 + ck3[i] * a32 + ck4[i] * a42;
        vgl3[i] = vgl3[i] + ck1[i] * a13 + ck2[i] * a23 + ck3[i] * a33 + ck4[i] * a43;
        vgl4[i] = vgl4[i] + ck1[i] * a14 + ck2[i] * a24 + ck3[i] * a34 + ck4[i] * a44;
        vgl5[i] = vgl5[i] + ck1[i] * a15 + ck2[i] * a25 + ck3[i] * a35 + ck4[i] * a45;
      }
    }

    const int64_t n0 = n < 0 ? 0 : n;
    for (int64_t n=n0 ; n < nidx ; n+=1) {
      const double* restrict ck = coef_normalized_t + idx[n]*mo_num;
      const double a1 = av1[n];
      const double a2 = av2[n];
      const double a3 = av3[n];
      const double a4 = av4[n];
      const double a5 = av5[n];

#ifdef HAVE_OPENMP
  #pragma omp simd
#endif
      for (int64_t i=0 ; i<mo_num ; ++i) {
        vgl1[i] += ck[i] * a1;
        vgl2[i] += ck[i] * a2;
        vgl3[i] += ck[i] * a3;
        vgl4[i] += ck[i] * a4;
        vgl5[i] += ck[i] * a5;
      }
    }
  }
  return QMCKL_SUCCESS;
}
#endif

2.1.5 Test

Author: TREX CoE

Created: 2022-05-05 Thu 18:50

Validate