TREX/qmckl

mirror of https://github.com/TREX-CoE/qmckl.git synced 2024-11-19 20:42:50 +01:00

Anthony Scemama 2ecfc55dbc Added eN cusp fitting in MOs

2023-03-17 14:54:58 +01:00

108 KiB

Raw Blame History

Molecular Orbitals

Context
Computation

The molecular orbitals (MOs) are defined in the basis of AOs along with a AO to MO coefficient matrix $C$. Using these coefficients (e.g. from Hartree Fock method) the MOs are defined as follows:

\[ \phi_i(\mathbf{r}) = C_i * \chi_i (\mathbf{r}) \]

By default, all the MOs are computed. A subset of MOs can be selected for computation, for example to remove computation of the useless virtual MOs present in a MO coefficient matrix.

For each nucleus, a radius $r_{\text{cusp}}$ can be given such that if the electron-nucleus distance is smaller than $r_{\text{cusp}}$, all the MOs are replaced by functions with the correct cusp condition and such that the values and the gradients of the MOs are continuous at $r_{\text{cusp}}$.

Molecular orbitals (MOs) are defined in the basis of atomic orbitals (AOs) using a coefficient matrix $C$, which determines how the AOs are combined to form the MOs.

The equation for the MOs $\phi_i(\mathbf{r})$ is given by:

$$\phi_i(\mathbf{r}) = \sum_{\mu} C_{\mu i} \chi_{\mu}(\mathbf{r})$$

where $i$ labels the MO, $\mu$ labels the AO, $C_{\mu i}$ is the coefficient of AO $\mu$ in MO $i$, and $\chi_{\mu}(\mathbf{r})$ is the AO itself.

In some cases, it may be desirable to only compute a subset of the MOs, for example to exclude virtual MOs that do not contribute to the wave function. This can be achieved by selecting a subset of columns from the coefficient matrix $C$.

The exact wave function must have a cusp at the positions of the nuclei to ensure that the kinetic energy diverges and cancels out the divergence of the potential, resulting in a finite total energy. To ensure that the cusp condition is satisfied, a modification can be made to the molecular orbitals (MOs) when the electron-nucleus distance is smaller than a certain radius $r_{\text{cusp}}$. At distances smaller than $r_{\text{cusp}}$, the MOs are replaced by functions that have the correct electron-nucleus cusp condition and that ensure that the values and the gradients of the MOs are continuous at $r_{\text{cusp}}$.

A radius $r_{\text{cusp}\, A}$ is given for each nucleus $A$, default is zero. If an electron is closer to the nucleus $A$ than $r_{\text{cusp}\, A}$, the MOs are locally replaced as \[ \phi_{\text{cusp}\, i}(\mathbf{r}) = \phi_i(\mathbf{r}) - \phi_{s_A i}(\mathbf{r}) + \sum_{k=0}^{3} f_k\, |\mathbf{r}-\mathbf{R}_A|^k \] where $\phi_{s_A i}$ are the contributions of the $s$ AOs centered at $A$ to MO $i$. The coefficients $f_k$ are such that the value and gradient of the MO are continuous at $r_{\text{cusp}}$, and the electron-nucleus cusp is exact.

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

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;
}

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

#+NAME:pre

#+NAME:post

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

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.

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

Fortran interfaces

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

Fortran interface

Computation

Parameters of the cusp-correction functions

Computation of MOs: values only

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

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;
}

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
}

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 );
#endif

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

#ifdef HAVE_OPENMP
#pragma omp parallel 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]);

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

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];

#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;
  }
}

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

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

Computation of MOs: values, gradient, Laplacian

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

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;
}

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

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, &
 coefficient_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)  :: coefficient_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

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

! info = qmckl_dgemm(context,'N', 'N', mo_num, point_num, ao_num, 1.d0, &
!      coefficient_t, int(size(coefficient_t,1),8),      &
!      ao_vgl, int(size(ao_vgl,1),8), 0.d0,                  &
!      mo_vgl, int(size(mo_vgl,1),8))

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* coefficient_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* coefficient_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* coefficient_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, coefficient_t, ao_vgl, mo_vgl);
#else
return qmckl_compute_mo_basis_mo_vgl_doc (context, ao_num, mo_num, point_num, coefficient_t, ao_vgl, mo_vgl);
#endif
}

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* coefficient_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 coefficient_t,
                               const double* restrict ao_vgl,
                               double* restrict const mo_vgl )
{
assert (context != QMCKL_NULL_CONTEXT);

#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) {
  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=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];

  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;
  }
}

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

#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

Computation of cusp-corrected MOs: values only

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[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	Cusp correction for the values of the MOs

integer function qmckl_compute_mo_basis_mo_value_cusp_doc_f(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)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: nucl_num, ao_num, mo_num, point_num
integer*8             , intent(in)  :: ao_nucl(ao_num)
integer*4             , intent(in)  :: ao_ang_mom(ao_num)
double precision      , intent(in)  :: en_distance(nucl_num, point_num)
double precision      , intent(in)  :: r_cusp(nucl_num)
double precision      , intent(in)  :: cusp_param(mo_num, 4, nucl_num)
double precision      , intent(in)  :: coefficient_t(mo_num, ao_num)
double precision      , intent(in)  :: ao_value(ao_num, point_num)
double precision      , intent(out) :: mo_value(mo_num, point_num)

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_f

qmckl_exit_code qmckl_compute_mo_basis_mo_value_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 double* r_cusp,
      const qmckl_tensor cusp_param,
      const double* coefficient_t,
      const double* ao_value,
      double* const mo_value );

qmckl_exit_code qmckl_compute_mo_basis_mo_value_cusp_doc (
      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* r_cusp,
      const double* cusp_param,
      const double* coefficient_t,
      const double* ao_value,
      double* const mo_value );

qmckl_exit_code
qmckl_compute_mo_basis_mo_value_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 double* r_cusp,
                                  const qmckl_tensor cusp_param_tensor,
                                  const double* coefficient_t,
                                  const double* ao_value,
                                  double* const mo_value )
{
qmckl_exit_code rc;

#ifdef HAVE_HPC
rc = qmckl_compute_mo_basis_mo_value_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num, 
                                             ao_nucl, ao_ang_mom, en_distance, r_cusp,
                                             cusp_param_tensor, coefficient_t, ao_value, mo_value ); 
#else
double * cusp_param  = qmckl_alloc_double_of_tensor(context, cusp_param_tensor);

rc = 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 ); 

qmckl_free(context, cusp_param);
#endif
return rc;
}

HPC version

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp_hpc (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* r_cusp,
                                      const qmckl_tensor cusp_param,
                                      const double* coefficient_t,
                                      const double* ao_value,
                                      double* const mo_value ); 
#endif

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_value_cusp_hpc (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* r_cusp,
                                      const qmckl_tensor cusp_param,
                                      const double* coefficient_t,
                                      const double* ao_value,
                                      double* const mo_value)
{
assert (context != QMCKL_NULL_CONTEXT);

#ifdef HAVE_OPENMP
#pragma omp parallel 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]);
const double* restrict ria   = &(en_distance[ipoint*nucl_num]);

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

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.) {
    const int64_t inucl = ao_nucl[k];
    if (ria[inucl] > r_cusp[inucl] || ao_ang_mom[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];

#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;
  }
}

for (int64_t m=n ; m < nidx ; m+=1) {
  const double* restrict ck = coefficient_t + idx[m]*mo_num;
  const double a1 = av1[m];
  
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
  for (int64_t i=0 ; i<mo_num ; ++i) {
    vgl1[i] += ck[i] * a1;
  }
}

for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
  if (ria[inucl] < r_cusp[inucl]) {
    const double r = ria[inucl];
#ifdef HAVE_OPENMP
#pragma omp simd
#endif
    for (int64_t i=0 ; i<mo_num ; ++i) {
      vgl1[i] += qmckl_ten3(cusp_param,i,0,inucl) + r*(
                 qmckl_ten3(cusp_param,i,1,inucl) + r*(
                 qmckl_ten3(cusp_param,i,2,inucl) + r*(
                 qmckl_ten3(cusp_param,i,3,inucl) )));
    }
  }
}

}
return QMCKL_SUCCESS;
}
#endif

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

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[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

integer function qmckl_compute_mo_basis_mo_vgl_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_vgl, mo_vgl) &
 result(info)
use qmckl
implicit none
integer(qmckl_context), intent(in)  :: context
integer*8             , intent(in)  :: nucl_num, ao_num, mo_num, point_num
integer*8             , intent(in)  :: ao_nucl(ao_num)
integer*4             , intent(in)  :: ao_ang_mom(ao_num)
double precision      , intent(in)  :: en_distance(nucl_num, point_num)
double precision      , intent(in)  :: nucl_coord(nucl_num,3)
double precision      , intent(in)  :: point_coord(point_num,3)
double precision      , intent(in)  :: r_cusp(nucl_num)
double precision      , intent(in)  :: cusp_param(mo_num,4,nucl_num)
double precision      , intent(in)  :: coefficient_t(mo_num,ao_num)
double precision      , intent(in)  :: ao_vgl(ao_num,5,point_num)
double precision      , intent(out) :: mo_vgl(mo_num,5,point_num)
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) cycle
    inucl = ao_nucl(k)+1
    if ( (en_distance(inucl,j) < r_cusp(inucl)) .and. (ao_ang_mom(k) == 0) ) cycle
    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 do
 
 ! Cusp adjustment
 do inucl=1,nucl_num
    r = en_distance(inucl,j)
    if (r > r_cusp(inucl)) cycle

    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

    enddo
 enddo ! inucl
end do
info = QMCKL_SUCCESS

end function qmckl_compute_mo_basis_mo_vgl_cusp_doc_f

qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_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_vgl,
      double* const mo_vgl );

qmckl_exit_code qmckl_compute_mo_basis_mo_vgl_cusp_doc (
      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_vgl,
      double* const mo_vgl );

qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_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_vgl,
                                double* const mo_vgl )
{
qmckl_exit_code rc;

#ifdef HAVE_HPC
rc = qmckl_compute_mo_basis_mo_vgl_cusp_hpc (context, nucl_num, ao_num, mo_num, point_num,
                                           ao_nucl, ao_ang_mom, en_distance, nucl_coord_matrix,
                                           point_coord_matrix, r_cusp, cusp_param_tensor,
                                           coefficient_t, ao_vgl, mo_vgl );
#else
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_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 );

qmckl_free(context, nucl_coord);
qmckl_free(context, point_coord);
qmckl_free(context, cusp_param);
#endif
return rc;
}

HPC version

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_cusp_hpc (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_vgl,
                                    double* const mo_vgl ); 
#endif

#ifdef HAVE_HPC
qmckl_exit_code
qmckl_compute_mo_basis_mo_vgl_cusp_hpc (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_vgl,
                                    double* const mo_vgl )
{
assert (context != QMCKL_NULL_CONTEXT);

#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.;
}

const double* restrict ria   = &(en_distance[ipoint*nucl_num]);

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) {
  if (avgl1[k] != 0.) {
    const int64_t inucl = ao_nucl[k];
    if (ria[inucl] > r_cusp[inucl] || ao_ang_mom[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=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];

  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;
  }
}

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

#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;
  }
}

// TODO
for (int64_t inucl=0 ; inucl<nucl_num ; ++inucl) {
  if (ria[inucl] < r_cusp[inucl]) {
    const double r = ria[inucl];
    const double r_vec[3] = {
      qmckl_mat(point_coord,ipoint,0) - qmckl_mat(nucl_coord,inucl,0),
      qmckl_mat(point_coord,ipoint,1) - qmckl_mat(nucl_coord,inucl,1),
      qmckl_mat(point_coord,ipoint,2) - qmckl_mat(nucl_coord,inucl,2) };
    const double r_inv = 1./r;

#ifdef HAVE_OPENMP
#pragma omp simd
#endif
    for (int64_t i=0 ; i<mo_num ; ++i) {
      vgl1[i] += qmckl_ten3(cusp_param,i,0,inucl) + r*(
                 qmckl_ten3(cusp_param,i,1,inucl) + r*(
                 qmckl_ten3(cusp_param,i,2,inucl) + r*(
                 qmckl_ten3(cusp_param,i,3,inucl) )));

      const double c1 = r_inv * qmckl_ten3(cusp_param,i,1,inucl) +
        2.0*qmckl_ten3(cusp_param,i,2,inucl) +
        r * 3.0 * qmckl_ten3(cusp_param,i,3,inucl);

       vgl2[i] += r_vec[0] * c1;
       vgl3[i] += r_vec[1] * c1;
       vgl4[i] += r_vec[2] * c1;

       vgl5[i] +=  2.0*qmckl_ten3(cusp_param,i,1,inucl)*r_inv +
                   6.0*qmckl_ten3(cusp_param,i,2,inucl) +
                  12.0*qmckl_ten3(cusp_param,i,3,inucl)*r;
    }
  }
}
}
return QMCKL_SUCCESS;
}
#endif

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

108 KiB Raw Blame History

Molecular Orbitals

Context

Data structure

Initialization functions

Cusp adjsutment functions

Access functions

Fortran interfaces

Update

Fortran interface

Computation

Parameters of the cusp-correction functions

Computation of MOs: values only

Get

Provide

Compute

HPC version

Computation of MOs: values, gradient, Laplacian

Get

Provide

Compute

HPC version

Computation of cusp-corrected MOs: values only

Compute

HPC version

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

Compute

HPC version

Rescaling of MO coefficients

Fortran interface

Test

108 KiB

Raw Blame History